path: root/src/algebra/special.spad.pamphlet

\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}