\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra liouv.spad}
\author{Manuel Bronstein}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package LF LiouvillianFunction}
<<package LF LiouvillianFunction>>=
)abbrev package LF LiouvillianFunction
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 10 August 1994
++ Keywords: liouvillian, function, primitive, exponential.
++ Examples:  )r LF INPUT
++ Description:
++   This package provides liouvillian functions over an integral domain.
LiouvillianFunction(R, F): Exports == Implementation where
  R:Join(OrderedSet, IntegralDomain)
  F:Join(FunctionSpace R,RadicalCategory,TranscendentalFunctionCategory)

  OP  ==> BasicOperator
  PR  ==> Polynomial R
  K   ==> Kernel F
  SE  ==> Symbol
  O   ==> OutputForm
  INP ==> InputForm
  INV ==> error "Invalid argument"

  SPECIALDIFF ==> "%specialDiff"
  SPECIALDISP ==> "%specialDisp"
  SPECIALINPUT ==> "%specialInput"
  SPECIALEQUAL ==> "%specialEqual"

  Exports ==> with
    belong? : OP -> Boolean
      ++ belong?(op) checks if op is Liouvillian
    operator: OP -> OP
      ++ operator(op) returns the Liouvillian operator based on op
    Ei      : F  -> F
      ++ Ei(f) denotes the exponential integral
    Si      : F  -> F
      ++ Si(f) denotes the sine integral
    Ci      : F  -> F
      ++ Ci(f) denotes the cosine integral
    li      : F  -> F
      ++ li(f) denotes the logarithmic integral
    erf     : F  -> F
      ++ erf(f) denotes the error function
    dilog   : F  -> F
      ++ dilog(f) denotes the dilogarithm
    integral: (F, SE) -> F
      ++ integral(f,x) indefinite integral of f with respect to x.
    integral: (F, SegmentBinding F) -> F
      ++ integral(f,x = a..b) denotes the definite integral of f with
      ++ respect to x from \spad{a} to b.

  Implementation ==> add
    iei        : F  -> F
    isi        : F  -> F
    ici        : F  -> F
    ierf       : F  -> F
    ili        : F  -> F
    ili2       : F  -> F
    iint       : List F -> F
    eqint      : (K,K) -> Boolean
    dvint      : (List F, SE) -> F
    dvdint     : (List F, SE) -> F
    ddint      : List F -> O
    integrand  : List F -> F

    dummy := new()$SE :: F

    opint  := operator("integral"::Symbol)$CommonOperators
    opdint := operator("%defint"::Symbol)$CommonOperators
    opei   := operator("Ei"::Symbol)$CommonOperators
    opli   := operator("li"::Symbol)$CommonOperators
    opsi   := operator("Si"::Symbol)$CommonOperators
    opci   := operator("Ci"::Symbol)$CommonOperators
    opli2  := operator("dilog"::Symbol)$CommonOperators
    operf  := operator("erf"::Symbol)$CommonOperators

    Si x                == opsi x
    Ci x                == opci x
    Ei x                == opei x
    erf x               == operf x
    li  x               == opli x
    dilog x             == opli2 x

    belong? op     == has?(op, "prim")
    isi x          == kernel(opsi, x)
    ici x          == kernel(opci, x)
    ierf x         == (zero? x => 0; kernel(operf, x))
--    ili2 x         == (one? x => INV; kernel(opli2, x))
    ili2 x         == ((x = 1) => INV; kernel(opli2, x))
    integrand l    == eval(first l, retract(second l)@K, third l)
    integral(f:F, x:SE) == opint [eval(f, k:=kernel(x)$K, dummy), dummy, k::F]

    iint l ==
      zero? first l => 0
      kernel(opint, l)

    ddint l ==
      int(integrand(l)::O * hconcat("d"::SE::O, third(l)::O),
                                    third(rest l)::O, third(rest rest l)::O)

    eqint(k1,k2) == 
      a1:=argument k1
      a2:=argument k2
      res:=operator k1 = operator k2
      if not res then return res
      res:= a1 = a2
      if res then return res
      res:= (a1.3 = a2.3) and (subst(a1.1,[retract(a1.2)@K],[a2.2]) = a2.1)

    dvint(l, x) ==
      k  := retract(second l)@K
      differentiate(third l, x) * integrand l
          + opint [differentiate(first l, x), second l, third l]


    dvdint(l, x) ==
      x = retract(y := third l)@SE => 0
      k := retract(d := second l)@K
      differentiate(h := third rest rest l,x) * eval(f := first l, k, h)
        - differentiate(g := third rest l, x) * eval(f, k, g)
             + opdint [differentiate(f, x), d, y, g, h]

    integral(f:F, s: SegmentBinding F) ==
      x := kernel(variable s)$K
      opdint [eval(f,x,dummy), dummy, x::F, lo segment s, hi segment s]

    ili x ==
      x = 1 => INV
      is?(x, "exp"::Symbol) => Ei first argument(retract(x)@K)
      kernel(opli, x)

    iei x ==
      x = 0 => INV
      is?(x, "log"::Symbol) => li first argument(retract(x)@K)
      kernel(opei, x)

    operator op ==
      is?(op, "integral"::Symbol)   => opint
      is?(op, "%defint"::Symbol)    => opdint
      is?(op, "Ei"::Symbol)         => opei
      is?(op, "Si"::Symbol)         => opsi
      is?(op, "Ci"::Symbol)         => opci
      is?(op, "li"::Symbol)         => opli
      is?(op, "erf"::Symbol)        => operf
      is?(op, "dilog"::Symbol)      => opli2
      error "Not a Liouvillian operator"

    evaluate(opei,   iei)$BasicOperatorFunctions1(F)
    evaluate(opli,   ili)
    evaluate(opsi,   isi)
    evaluate(opci,   ici)
    evaluate(operf,  ierf)
    evaluate(opli2,  ili2)
    evaluate(opint,  iint)
    derivative(opsi, sin(#1) / #1)
    derivative(opci, cos(#1) / #1)
    derivative(opei, exp(#1) / #1)
    derivative(opli, inv log(#1))
    derivative(operf, 2 * exp(-(#1**2)) / sqrt(pi()))
    derivative(opli2, log(#1) / (1 - #1))
    setProperty(opint,SPECIALEQUAL,eqint@((K,K) -> Boolean) pretend None)
    setProperty(opint,SPECIALDIFF,dvint@((List F,SE) -> F) pretend None)
    setProperty(opdint,SPECIALDIFF,dvdint@((List F,SE)->F) pretend None)
    setProperty(opdint, SPECIALDISP, ddint@(List F -> O) pretend None)

    if R has ConvertibleTo INP then
      inint : List F -> INP
      indint: List F -> INP
      pint  : List INP -> INP


      pint l  == convert concat(convert("integral"::SE)@INP, l)
      inint l == 
        r2:= convert([convert("::"::SE)@INP,convert(third l)@INP,convert("Symbol"::SE)@INP]@List INP)@INP
        pint [convert(integrand l)@INP, r2]

      indint l ==
        pint [convert(integrand l)@INP,
              convert concat(convert("="::SE)@INP,
                            [convert(third l)@INP,
                             convert concat(convert("SEGMENT"::SE)@INP,
                                            [convert(third rest l)@INP,
                                             convert(third rest rest l)@INP])])]

      setProperty(opint, SPECIALINPUT, inint@(List F -> INP) pretend None)
      setProperty(opdint, SPECIALINPUT, indint@(List F -> INP) pretend None)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

-- SPAD files for the functional world should be compiled in the
-- following order:
--
--   op  kl  fspace  algfunc elemntry LIOUV expr

<<package LF LiouvillianFunction>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}