Initial population.

1 files changed, 456 insertions, 0 deletions

diff --git a/src/algebra/special.spad.pamphlet b/src/algebra/special.spad.pamphlet
new file mode 100644
index 00000000..1765eebd
--- /dev/null
+++ b/src/algebra/special.spad.pamphlet
@@ -0,0 +1,456 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra special.spad}
+\author{Bruce W. Char, Stephen M. Watt}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package DFSFUN DoubleFloatSpecialFunctions}
+<<package DFSFUN DoubleFloatSpecialFunctions>>=
+)abbrev package DFSFUN DoubleFloatSpecialFunctions
+++ Author: Bruce W. Char, Stephen M. Watt
+++ Date Created:  1990
+++ Date Last Updated: June 25, 1991
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:
+++   This package provides special functions for double precision
+++   real and complex floating point.
+
+DoubleFloatSpecialFunctions(): Exports == Impl where
+    NNI ==> NonNegativeInteger
+    R   ==> DoubleFloat
+    C   ==> Complex DoubleFloat
+
+    Exports ==> with
+	Gamma: R -> R
+	    ++ Gamma(x) is the Euler gamma function, \spad{Gamma(x)}, defined by
+	    ++   \spad{Gamma(x) = integrate(t^(x-1)*exp(-t), t=0..%infinity)}.
+	Gamma: C -> C
+	    ++ Gamma(x) is the Euler gamma function, \spad{Gamma(x)}, defined by
+	    ++   \spad{Gamma(x) = integrate(t^(x-1)*exp(-t), t=0..%infinity)}.
+
+	Beta: (R, R) -> R
+	    ++ Beta(x, y) is the Euler beta function, \spad{B(x,y)}, defined by
+            ++   \spad{Beta(x,y) = integrate(t^(x-1)*(1-t)^(y-1), t=0..1)}.
+	    ++ This is related to \spad{Gamma(x)} by
+	    ++   \spad{Beta(x,y) = Gamma(x)*Gamma(y) / Gamma(x + y)}.
+	Beta: (C, C) -> C
+	    ++ Beta(x, y) is the Euler beta function, \spad{B(x,y)}, defined by
+            ++   \spad{Beta(x,y) = integrate(t^(x-1)*(1-t)^(y-1), t=0..1)}.
+	    ++ This is related to \spad{Gamma(x)} by
+	    ++   \spad{Beta(x,y) = Gamma(x)*Gamma(y) / Gamma(x + y)}.
+
+	logGamma: R -> R
+	    ++ logGamma(x) is the natural log of \spad{Gamma(x)}.
+	    ++ This can often be computed even if \spad{Gamma(x)} cannot.
+	logGamma: C -> C
+	    ++ logGamma(x) is the natural log of \spad{Gamma(x)}.
+	    ++ This can often be computed even if \spad{Gamma(x)} cannot.
+
+	digamma: R -> R
+	    ++ digamma(x) is the function, \spad{psi(x)}, defined by
+	    ++   \spad{psi(x) = Gamma'(x)/Gamma(x)}.
+	digamma: C -> C
+	    ++ digamma(x) is the function, \spad{psi(x)}, defined by
+	    ++   \spad{psi(x) = Gamma'(x)/Gamma(x)}.
+
+	polygamma: (NNI, R) -> R
+	    ++ polygamma(n, x) is the n-th derivative of \spad{digamma(x)}.
+	polygamma: (NNI, C) -> C
+	    ++ polygamma(n, x) is the n-th derivative of \spad{digamma(x)}.
+
+
+	besselJ: (R,R) -> R
+	    ++ besselJ(v,x) is the Bessel function of the first kind,
+            ++ \spad{J(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
+	besselJ: (C,C) -> C
+	    ++ besselJ(v,x) is the Bessel function of the first kind,
+            ++ \spad{J(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
+
+	besselY: (R, R) -> R
+	    ++ besselY(v,x) is the Bessel function of the second kind,
+            ++ \spad{Y(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
+            ++ Note: The default implmentation uses the relation
+            ++   \spad{Y(v,x) = (J(v,x) cos(v*%pi) - J(-v,x))/sin(v*%pi)}
+            ++ so is not valid for integer values of v.
+	besselY: (C, C) -> C
+	    ++ besselY(v,x) is the Bessel function of the second kind,
+            ++ \spad{Y(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
+            ++ Note: The default implmentation uses the relation
+            ++   \spad{Y(v,x) = (J(v,x) cos(v*%pi) - J(-v,x))/sin(v*%pi)}
+            ++ so is not valid for integer values of v.
+
+	besselI: (R,R) -> R
+	    ++ besselI(v,x) is the modified Bessel function of the first kind,
+            ++ \spad{I(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
+	besselI: (C,C) -> C
+	    ++ besselI(v,x) is the modified Bessel function of the first kind,
+            ++ \spad{I(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
+
+	besselK: (R, R) -> R
+	    ++ besselK(v,x) is the modified Bessel function of the first kind,
+            ++ \spad{K(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
+            ++ Note: The default implmentation uses the relation
+            ++   \spad{K(v,x) = %pi/2*(I(-v,x) - I(v,x))/sin(v*%pi)}.
+            ++ so is not valid for integer values of v.
+	besselK: (C, C) -> C
+	    ++ besselK(v,x) is the modified Bessel function of the first kind,
+            ++ \spad{K(v,x)}.
+	    ++ This function satisfies the differential equation:
+	    ++   \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
+            ++ Note: The default implmentation uses the relation
+            ++   \spad{K(v,x) = %pi/2*(I(-v,x) - I(v,x))/sin(v*%pi)}
+            ++ so is not valid for integer values of v.
+
+        airyAi:   C -> C
+            ++ airyAi(x) is the Airy function \spad{Ai(x)}.
+            ++ This function satisfies the differential equation:
+            ++   \spad{Ai''(x) - x * Ai(x) = 0}.
+        airyAi:   R -> R
+            ++ airyAi(x) is the Airy function \spad{Ai(x)}.
+            ++ This function satisfies the differential equation:
+            ++   \spad{Ai''(x) - x * Ai(x) = 0}.
+
+        airyBi:   R -> R
+            ++ airyBi(x) is the Airy function \spad{Bi(x)}.
+            ++ This function satisfies the differential equation:
+            ++   \spad{Bi''(x) - x * Bi(x) = 0}.
+        airyBi:   C -> C
+            ++ airyBi(x) is the Airy function \spad{Bi(x)}.
+            ++ This function satisfies the differential equation:
+            ++   \spad{Bi''(x) - x * Bi(x) = 0}.
+
+	hypergeometric0F1: (R, R) -> R
+	    ++ hypergeometric0F1(c,z) is the hypergeometric function
+            ++ \spad{0F1(; c; z)}.
+	hypergeometric0F1: (C, C) -> C
+	    ++ hypergeometric0F1(c,z) is the hypergeometric function
+            ++ \spad{0F1(; c; z)}.
+
+
+    Impl ==> add
+	a, v, w, z: C
+	n, x, y: R
+
+        -- These are hooks to Bruce's boot code.
+	Gamma z         == CGAMMA(z)$Lisp
+	Gamma x         == RGAMMA(x)$Lisp
+
+	polygamma(k,z)  == CPSI(k, z)$Lisp
+	polygamma(k,x)  == RPSI(k, x)$Lisp
+
+	logGamma z      == CLNGAMMA(z)$Lisp
+	logGamma x      == RLNGAMMA(x)$Lisp
+
+	besselJ(v,z)    == CBESSELJ(v,z)$Lisp
+	besselJ(n,x)    == RBESSELJ(n,x)$Lisp
+
+	besselI(v,z)    == CBESSELI(v,z)$Lisp
+	besselI(n,x)    == RBESSELI(n,x)$Lisp
+
+	hypergeometric0F1(a,z) == CHYPER0F1(a, z)$Lisp
+	hypergeometric0F1(n,x) == retract hypergeometric0F1(n::C, x::C)
+
+
+        -- All others are defined in terms of these.
+	digamma x == polygamma(0, x)
+	digamma z == polygamma(0, z)
+
+	Beta(x,y) == Gamma(x)*Gamma(y)/Gamma(x+y)
+	Beta(w,z) == Gamma(w)*Gamma(z)/Gamma(w+z)
+
+        fuzz := (10::R)**(-7)
+
+        import IntegerRetractions(R)
+        import IntegerRetractions(C)
+
+        besselY(n,x) ==
+            if integer? n then n := n + fuzz
+            vp := n * pi()$R
+            (cos(vp) * besselJ(n,x) - besselJ(-n,x) )/sin(vp)
+	besselY(v,z) ==
+            if integer? v then v := v + fuzz::C
+            vp := v * pi()$C
+            (cos(vp) * besselJ(v,z) - besselJ(-v,z) )/sin(vp)
+
+        besselK(n,x) ==
+            if integer? n then n := n + fuzz
+            p    := pi()$R
+            vp   := n*p
+            ahalf:= 1/(2::R)
+            p * ahalf * ( besselI(-n,x) - besselI(n,x) )/sin(vp)
+        besselK(v,z) ==
+            if integer? v then v := v + fuzz::C
+            p    := pi()$C
+            vp   := v*p
+            ahalf:= 1/(2::C)
+            p * ahalf * ( besselI(-v,z) - besselI(v,z) )/sin(vp)
+
+        airyAi x ==
+            ahalf  := recip(2::R)::R
+            athird := recip(3::R)::R
+            eta := 2 * athird * (-x) ** (3*ahalf)
+            (-x)**ahalf * athird * (besselJ(-athird,eta) + besselJ(athird,eta))
+        airyAi z ==
+            ahalf  := recip(2::C)::C
+            athird := recip(3::C)::C
+            eta := 2 * athird * (-z) ** (3*ahalf)
+            (-z)**ahalf * athird * (besselJ(-athird,eta) + besselJ(athird,eta))
+
+        airyBi x ==
+            ahalf  := recip(2::R)::R
+            athird := recip(3::R)::R
+            eta := 2 * athird * (-x) ** (3*ahalf)
+            (-x*athird)**ahalf * ( besselJ(-athird,eta) - besselJ(athird,eta) )
+
+        airyBi z ==
+            ahalf  := recip(2::C)::C
+            athird := recip(3::C)::C
+            eta := 2 * athird * (-z) ** (3*ahalf)
+            (-z*athird)**ahalf * ( besselJ(-athird,eta) - besselJ(athird,eta) )
+
+@
+\section{package ORTHPOL OrthogonalPolynomialFunctions}
+<<package ORTHPOL OrthogonalPolynomialFunctions>>=
+)abbrev package ORTHPOL OrthogonalPolynomialFunctions
+++ Author: Stephen M. Watt
+++ Date Created:  1990
+++ Date Last Updated: June 25, 1991
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:
+++   This package provides orthogonal polynomials as functions on a ring.
+
+OrthogonalPolynomialFunctions(R: CommutativeRing): Exports == Impl where
+    NNI ==> NonNegativeInteger
+    RN  ==> Fraction Integer
+
+    Exports ==> with
+
+        chebyshevT: (NNI, R) -> R
+           ++ chebyshevT(n,x) is the n-th Chebyshev polynomial of the first
+           ++ kind, \spad{T[n](x)}.  These are defined by
+           ++ \spad{(1-t*x)/(1-2*t*x+t**2) = sum(T[n](x) *t**n, n = 0..)}.
+
+        chebyshevU: (NNI, R) -> R
+           ++ chebyshevU(n,x) is the n-th Chebyshev polynomial of the second
+           ++ kind, \spad{U[n](x)}. These are defined by
+           ++ \spad{1/(1-2*t*x+t**2) = sum(T[n](x) *t**n, n = 0..)}.
+
+        hermiteH:   (NNI, R) -> R
+           ++ hermiteH(n,x) is the n-th Hermite polynomial, \spad{H[n](x)}.
+           ++ These are defined by
+           ++ \spad{exp(2*t*x-t**2) = sum(H[n](x)*t**n/n!, n = 0..)}.
+
+        laguerreL:  (NNI, R) -> R
+           ++ laguerreL(n,x) is the n-th Laguerre polynomial, \spad{L[n](x)}.
+           ++ These are defined by
+           ++ \spad{exp(-t*x/(1-t))/(1-t) = sum(L[n](x)*t**n/n!, n = 0..)}.
+
+        laguerreL:  (NNI, NNI, R) -> R
+           ++ laguerreL(m,n,x) is the associated Laguerre polynomial,
+           ++ \spad{L<m>[n](x)}.  This is the m-th derivative of \spad{L[n](x)}.
+
+        if R has Algebra RN then
+            legendreP:  (NNI, R) -> R
+                ++ legendreP(n,x) is the n-th Legendre polynomial,
+                ++ \spad{P[n](x)}.  These are defined by
+                ++ \spad{1/sqrt(1-2*x*t+t**2) = sum(P[n](x)*t**n, n = 0..)}.
+
+    Impl ==> add
+        p0, p1: R
+        cx:     Integer
+
+        import IntegerCombinatoricFunctions()
+
+        laguerreL(n, x) ==
+            n = 0 => 1
+            (p1, p0) := (-x + 1, 1)
+            for i in 1..n-1 repeat
+                (p1, p0) := ((2*i::R + 1 - x)*p1 - i**2*p0, p1)
+            p1
+        laguerreL(m, n, x) ==
+            ni := n::Integer
+            mi := m::Integer
+            cx := (-1)**m * binomial(ni,ni-mi) * factorial(ni)
+            p0 := 1
+            p1 := cx::R
+            for j in 1..ni-mi repeat
+                cx := -cx*(ni-mi-j+1)
+                cx := (cx exquo ((mi+j)*j))::Integer
+                p0 := p0 * x
+                p1 := p1 + cx*p0
+            p1
+        chebyshevT(n, x) ==
+            n = 0 => 1
+            (p1, p0) := (x, 1)
+            for i in 1..n-1 repeat
+                (p1, p0) := (2*x*p1 - p0, p1)
+            p1
+        chebyshevU(n, x) ==
+            n = 0 => 1
+            (p1, p0) := (2*x, 1)
+            for i in 1..n-1 repeat
+                (p1, p0) := (2*x*p1 - p0, p1)
+            p1
+        hermiteH(n, x) ==
+            n = 0 => 1
+            (p1, p0) := (2*x, 1)
+            for i in 1..n-1 repeat
+                (p1, p0) := (2*x*p1 - 2*i*p0, p1)
+            p1
+        if R has Algebra RN then
+            legendreP(n, x) ==
+                n = 0 => 1
+                p0 := 1
+                p1 := x
+                for i in 1..n-1 repeat
+                    c: RN := 1/(i+1)
+                    (p1, p0) := (c*((2*i+1)*x*p1 - i*p0), p1)
+                p1
+
+@
+\section{package NTPOLFN NumberTheoreticPolynomialFunctions}
+<<package NTPOLFN NumberTheoreticPolynomialFunctions>>=
+)abbrev package NTPOLFN NumberTheoreticPolynomialFunctions
+++ Author: Stephen M. Watt
+++ Date Created:  1990
+++ Date Last Updated: June 25, 1991
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:
+++   This package provides polynomials as functions on a ring.
+
+NumberTheoreticPolynomialFunctions(R: CommutativeRing): Exports == Impl where
+    NNI ==> NonNegativeInteger
+    RN  ==> Fraction Integer
+
+    Exports ==> with
+
+        cyclotomic: (NNI, R) -> R
+		++ cyclotomic(n,r) \undocumented
+
+        if R has Algebra RN then
+            bernoulliB: (NNI, R) -> R
+		++ bernoulliB(n,r) \undocumented
+            eulerE:     (NNI, R) -> R
+		++ eulerE(n,r) \undocumented
+
+    Impl ==> add
+
+        import PolynomialNumberTheoryFunctions()
+
+        I   ==> Integer
+        SUP ==> SparseUnivariatePolynomial
+
+        -- This is the wrong way to evaluate the polynomial.
+        cyclotomic(k, x) ==
+            p: SUP(I) := cyclotomic(k)
+            r: R      := 0
+            while p ^= 0 repeat
+                d := degree p
+                c := leadingCoefficient p
+                p := reductum p
+                r := c*x**d + r
+            r
+
+        if R has Algebra RN then
+            eulerE(k, x) ==
+                p: SUP(RN) := euler(k)
+                r: R       := 0
+                while p ^= 0 repeat
+                    d := degree p
+                    c := leadingCoefficient p
+                    p := reductum p
+                    r := c*x**d + r
+                r
+            bernoulliB(k, x) ==
+                p: SUP(RN) := bernoulli(k)
+                r: R       := 0
+                while p ^= 0 repeat
+                    d := degree p
+                    c := leadingCoefficient p
+                    p := reductum p
+                    r := c*x**d + r
+                r
+
+@
+\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>>
+
+<<package DFSFUN DoubleFloatSpecialFunctions>>
+<<package ORTHPOL OrthogonalPolynomialFunctions>>
+<<package NTPOLFN NumberTheoreticPolynomialFunctions>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}