\documentclass{article} \usepackage{open-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} <>= )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} <>= )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[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 p0: R := 1 p1: R := -x + 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 p0: R := 1 p1: R := x for i in 1..n-1 repeat (p1, p0) := (2*x*p1 - p0, p1) p1 chebyshevU(n, x) == n = 0 => 1 p0: R := 1 p1: R := 2*x for i in 1..n-1 repeat (p1, p0) := (2*x*p1 - p0, p1) p1 hermiteH(n, x) == n = 0 => 1 p0: R := 1 p1: R := 2*x 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} <>= )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} <>= --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. @ <<*>>= <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}