#lang scribble/doc @(require scribble/manual scribble/struct scribblings/icons (for-label scheme/base "../science-with-graphics.ss")) @title[#:tag "special-functions"]{Special Functions} @local-table-of-contents[] This chapter describes the special functions provided by the PLT Scheme Science Collection. The functions described in this chapter are defined in the @schemefont{special-functions} sub-collection of the science collection. The entire @schemefont{special-functions} sub-collection can be made available using the form: @defmodule[(planet williams/science/special-functions)] The individual modules in the @schemefont{special-functions} sub-collection can also be made available as describes in the sections below. @section{Error Functions} The @tech{error function} is described in Abramowitz and Stegun @cite["Abramowitz64"], Chapter 7. The functions are defined in the @filepath{error.ss} file in the @schemefont{special-functions} sub-collection of the science collection and are made available using the form: @defmodule[(planet williams/science/special-functions/error)] @subsection{Error Function} @(margin-note magnify @link["http://mathworld.wolfram.com/Erf.html"]{Erf} " from Wolfram MathWorld.") @defproc[(erf (x real?)) (real-in -1.0 1.0)]{ Computes the error function: @image["scribblings/images/erf-equation.gif"]{erf equation}.} @bold{Example:} Plot of @scheme[(erf x)] on the interval [-4, 4]. @schememod[ scheme (require (planet williams/science/special-functions/error)) (require plot/plot) (plot (line erf) (x-min -4.0) (x-max 4.0) (y-min -1.0) (y-max 1.0) (title "Error Function, erf(x)"))] The following figure shows the resulting plot: @image["scribblings/images/erf.bmp"]{erf plot} @subsection{Complementary Error Function} @(margin-note magnify @link["http://mathworld.wolfram.com/Erfc.html"]{Erfc} " from Wolfram MathWorld.") @defproc[(erfc (x real?)) (real-in 0.0 2.0)]{ Computes the complementary error function: @image["scribblings/images/erfc-equation.gif"]{erfc equation}.} @bold{Example:} Plot of @scheme[(erfc x)] on the interval [-4, 4]. @schememod[ scheme (require (planet williams/science/special-functions/error)) (require plot/plot) (plot (line erfc) (x-min -4.0) (x-max 4.0) (y-min 0.0) (y-max 2.0) (title "Complementary Error Function, erfc(x)"))] The following figure shows the resulting plot: @image["scribblings/images/erfc.bmp"]{erfc plot} @subsection{Hazard Function} @(margin-note magnify @link["http://mathworld.wolfram.com/HazardFunction.html"]{Hazard Function} " from Wolfram MathWorld.") The @tech{hazard function} for the normal distribution, also known as the inverse Mill's ratio, is the ratio of the probability function, @math{P(x)}, to the survival function, @math{S(x)}, and is defined as: @image["scribblings/images/hazard-equation.gif"]{hazard equation} @defproc[(hazard (x real?)) (>=/c 0.0)]{ Computes the hazard function for the normal distribution.} @bold{Example:} Plot of @scheme[(hazard x)] on the interval [-5, 10]. @schememod[ scheme (require (planet williams/science/special-functions/error)) (require plot/plot) (plot (line hazard) (x-min -5.0) (x-max 10.0) (y-min 0.0) (y-max 10.0) (title "Hazard Function, hazard(x)"))] The following figure shows the resulting plot: @image["scribblings/images/hazard.bmp"]{hazard plot} @section{Exponential Integral Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/ExponentialIntegral.html"]{Exponential Intgral} " from Wolfram MathWorld.") Information on the exponential integral functions can be found in Abramowitz and Stegun @cite["Abramowitz64"], Chapter 5. The functions are defined in the @filepath{exponential-integral.ss} file in the @schemefont{special-functions} sub-collection of the science collection are made available using the form: @defmodule[(planet williams/science/special-functions/exponential-integral)] @subsection{First-Order Exponential Integral} @defproc*[(((expint-E1 (x real?)) real?) ((expint-E1-scaled (x real?)) real?))]{ Computes the exponential integral @math{E_1(x)}: @image["scribblings/images/expint-E1-equation.gif"]} @bold{Example:} Plot of @scheme[(expint-E1 x)] on the interval [-4, 4]. @schememod[ scheme (require (planet williams/science/special-functions/exponential-integral)) (require plot/plot) (plot (line expint-E1) (x-min -4.0) (x-max 4.0) (y-min -10.0) (y-max 10.0) (title "Exponential Integral, E1(x)"))] The following figure shows the resulting plot: @image["scribblings/images/expint-E1.gif"] @subsection{Second-Order Exponential Integral} @defproc*[(((expint-E2 (x real?)) real?) ((expint-E2-scaled (x real?)) real?))]{ Computes the second-order exponential integral @math{E_2(x)}: @image["scribblings/images/expint-E2-equation.gif"]} @bold{Example:} Plot of @scheme[(expint-E2 x)] on the interval [-4, 4]. @schememod[ scheme (require (planet williams/science/special-functions/exponential-integral)) (require plot/plot) (plot (line expint-E2) (x-min -4.0) (x-max 4.0) (y-min -10.0) (y-max 10.0) (title "Exponential Integral, E2(x)"))] The following figure shows the resulting plot: @image["scribblings/images/expint-E2.gif"] @subsection{General Exponential Integral} @defproc*[(((expint-Ei (x real?)) real?) ((expint-Ei-scaled (x real?)) real?))]{ Computes the exponential integral @math{E_i(x)}: @image["scribblings/images/expint-Ei-equation.gif"]} @bold{Example:} Plot of @scheme[(expint-Ei x)] on the interval [-4, 4]. @schememod[ scheme (require (planet williams/science/special-functions/exponential-integral)) (require plot/plot) (plot (line expint-Ei) (x-min -4.0) (x-max 4.0) (y-min -10.0) (y-max 10.0) (title "Exponential Integral, Ei(x)"))] The following figure shows the resulting plot: @image["scribblings/images/expint-Ei.gif"] @section{Gamma and Beta Functions} The gamma functions are defined in the @filepath{gamma.ss} file in the @schemefont{special-functions} sub-collection of the science collection and are made available using the form: @defmodule[(planet williams/science/special-functions/gamma)] Note that the gamma functions (Section 5.3), psi functions (Section 5.4), and the zeta functions (Section 5.5) are defined in the same module, @filepath{gamma.ss}. This is because their definitions are interdependent and PLT Scheme does not allow circular module dependencies. @subsection{Gamma Function} @(margin-note magnify @link["http://mathworld.wolfram.com/GammaFunction.html"]{Gamma Function} " from Wolfram MathWorld.") The @tech{gamma function} is defined by the integral: @image["scribblings/images/gamma-equation.gif"] It is related to the factorial function by @math{Γ(n) = (n - 1)!} for positive integer @math{n}. Further information on the gamma function can be found in Abramowitz and Stegun @cite["Abramowitz64"], Chapter 6. @defproc[(gamma (x real?)) real?]{ Computes the gamma function, @math{Γ(x)}, subject to @scheme[x] not being a negative integer. This function is computed using the real Lanczos method. The maximum value of @scheme[x] such that @math{Γ(x)} is not considered an overflow is given by the constant @scheme[gamma-xmax] and is 171.0.} @bold{Example:} Plot of @scheme[(gamma x)] on the interval (0, 6]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line gamma) (x-min 0.001) (x-max 6.0) (y-min 0.0) (y-max 120.0) (title "Gamma Function, Gamma(x)"))] The following figure shows the resulting plot: @image["scribblings/images/gamma.bmp"] @bold{Example:} Plot of @scheme[(gamma x)] on the interval (-1, 0). @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line gamma) (x-min -0.999) (x-max -0.001) (y-min -120.0) (y-max 0.0) (title "Gamma Function, Gamma(x)"))] The following figure shows the resulting plot: @image["scribblings/images/gamma-1.bmp"] @(margin-note magnify @link["http://mathworld.wolfram.com/LogGammaFunction.html"]{Log Gamma Function} " from Wolfram MathWorld.") @defproc[(lngamma (x real?)) real?]{ Computes the logarithm of the gamma function, @math{log Γ(x)}, subject to @scheme[x] not being a negative integer. For @math{x < 0}, the real part of @math{log Γ(x)} is returned, which is equivalent to @math{log |Γ(x)|}. The function is computed using the real Lanczos method.} @defproc[(lngamma-sgn (x real?)) (values real? (integer-in -1 1))]{ Computes the logarithm of the magnitude of the gamma function and its sign, subject to @scheme[x] not being a negative integer, and returns them as multiple values. The function is computed using the real Lanczos method. The value of the gamma function can be reconstructed using the relation @math{Γ(x) = sgn × exp(resultlg)}, where @math{resultlg} and @math{sgn} are the returned values.} @bold{Example:} Plot of @scheme[(lngamma x)] on the interval (0, 6]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line lngamma) (x-min 0.001) (x-max 6.0) (y-min 0.0) (y-max 5.0) (title "Log Gamma Function, log Gamma(x)"))] The following figure shows the resulting plot: @image["scribblings/images/lngamma.bmp"] @defproc[(gamma-inv (x real?)) real?]{ Computes the reciprocal of the gamma function, @math{1/Γ(x)}, using the real Lanczos method.} @subsection{Regulated Gamma Function} The @tech{regulated gamma function} is given by @image["scribblings/images/gammastar-equation.gif"] @defproc[(gammastar (x (>/c 0.0))) real?]{ Computes the regulated gamma function, @math{Γ*(x)}, for @math{x > 0}.} @bold{Example:} Plot of @scheme[(gammastar x)] on the interval (0, 4]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line gammastar) (x-min 0.001) (x-max 4.0) (y-min 0.0) (y-max 10.0) (title "Regulated Gamma Function, Gamma*(x)"))] The following figure shows the resulting plot: @image["scribblings/images/gammastar.gif"] @subsection{Incomplete Gamma Function} @(margin-note magnify @link["http://mathworld.wolfram.com/IncompleteGammaFunction.html"]{Incomplete Gamma Function} " from Wolfram MathWorld.") @defproc[(gamma-inc-Q (a (>/c 0.0)) (x (>=/c 0.0))) real?]{ Computes the normalized incomplete gamma function, @image["scribblings/images/gamma-inc-Q-equation.gif"] for @math{a > 0} and @math{x ≥ 0}.} @defproc[(gamma-inc-P (a (>/c 0.0)) (x (>=/c 0.0))) real?]{ Computes the complementary normalized incomplete gamma function, @image["scribblings/images/gamma-inc-P-equation.gif"] for @math{a > 0} and @math{x ≥ 0}.} @defproc[(gamma-inc (a real?) (x (>=/c 0.0))) real?]{ Computes the unnormalized incomplete gamma function, @image["scribblings/images/gamma-inc-equation.gif"] for @scheme[a] real and @math{x ≥ 0}.} @subsection{Factorial Function} The factorial of a positive integer @math{n}, @math{n!}, is defined as @math{n! = n × (n - 1) × ... × 2 × 1}. By definition, @math{0! = 1}. The related function is related to the gamma function by @math{Γ(n) = (n - 1)!}. @defproc[(fact (n natural-number/c)) (>=/c 1.0)]{ Computes the factorial of @scheme[n], @math{n!}.} @defproc[(lnfact (n natural-number/c)) (>=/c 0.0)]{ Computes the logarithm of the factorial of @scheme[n], @math{log n!}. The algorithm is faster than computing @math{ln Γ(n + 1)} via @scheme[lngamma] for @math{n < 170}, but defers for larger @math{n}.} @subsection{Double Factorial Function} The double factorial of @scheme[n], @math{n!!}, is defined as @math{n! = n × (n - 2) × (n - 4) × ...}. By definition, @math{-1!! = 0!! = 1}. @defproc[(double-fact (n natural-number/c)) (>=/c 1.0)]{ Computes the double factorial of @scheme[n], @math{n!!}.} @defproc[(lndouble-fact (n natural-number/c)) (>=/c 0.0)]{ Computes the logarithm of the double factorial of @scheme[n], @math{log n!!}.} @subsection{Binomial Coefficient Function} The binomial coefficient, @math{n choose m}, is defined as: @image["scribblings/images/choose-equation.gif"] @defproc[(choose (n natural-number/c) (m natural-number/c)) (>=/c 1.0)]{ Computes the binomial coefficient for @scheme[n] and @scheme[m], @math{n choose m}.} @defproc[(lnchoose (n natural-number/c) (m natural-number/c)) (>=/c 0.0)]{ Computes the logarithm of the binomial coefficient for @scheme[n] and @scheme[m], @math{log (n choose m)}.} @subsection{Beta Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/BetaFunction.html"]{Beta Function} " from Wolfram MathWorld.") @defproc[(beta (a (>/c 0.0)) (b (>/c 0.0))) real?]{ Computes the beta function, @image["scribblings/images/beta-equation.gif"] for @math{a > 0} and @math{b > 0}.} @defproc[(lnbeta (a (>/c 0.0)) (b (>/c 0.0))) real?]{ Computes the logarithm of the beta function, @math{log Β(a, b)}, for @math{a > 0} and @math{b > 0}.} @section{Psi Functions} The psi functions are defined in the @filepath{gamma.ss} file in the @schemefont{special-functions} sub-collection of the science collection and are made available using the form: @;@defmodule[(planet williams/science/special-functions/gamma)] @scheme[(require (planet williams/science/special-functions/gamma))] Note that the gamma functions (Section 5.3), psi functions (Section 5.4), and the zeta functions (Section 5.5) are defined in the same module, @filepath{gamma.ss}. This is because their definitions are interdependent and PLT Scheme does not allow circular module dependencies. @subsection{Psi (Digamma) Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/DigammaFunction.html"]{Digamma Function} " from Wolfram MathWorld.") @defproc[(psi-int (n (integer-in 1 +inf.0))) real?]{ Computes the digamma function, @math{Ψ(n)}, for positive integer @scheme[n]. The digamma function is also called the Psi function.} @defproc[(psi (x real?)) real?]{ Returns the digamma function, @math{Ψ(x)}, for general @scheme[x], @math{x ≠ 0}.} @defproc[(psi-1piy (y real?)) real?]{ Computes the real part of the digamma function on the line @math{1 + iy}, @math{Re Ψ(1 + iy)}.} @bold{Example:} Plot of @scheme[(psi x)] on the interval (0, 5]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line psi) (x-min 0.001) (x-max 5.0) (y-min -5.0) (y-max 5.0) (title "Psi (Digamma) Function, Psi(x)"))] The following figure shows the resulting plot: @image["scribblings/images/psi.bmp"] @subsection{Psi-1 (Trigamma) Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/TrigammaFunction.html"]{Trigamma Function} " from Wolfram MathWorld.") @defproc[(psi-1-int (n (integer-in 1 +inf.0))) real?]{ Computes the trigamma function, @math{Ψ'(n)}, for positive integer @scheme[n].} @defproc[(psi-1 (x real?)) real?]{ Computes the trigamma function, @math{Ψ'(x)}, for general @scheme[x].} @bold{Example:} Plot of @scheme[(psi-1 x)] on the interval (0, 5]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line psi-1) (x-min 0.001) (x-max 5.0) (y-min 0.0) (y-max 5.0) (title "Psi-1 (Trigamma) Function, Psi-1(x)"))] The following figure shows the resulting plot: @image["scribblings/images/psi-1.bmp"] @subsection{Psi-n (Polygamma) Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/PolygammaFunction.html"]{Polygamma Function} " from Wolfram MathWorld.") @defproc[(psi-n (n natural-number/c) (x (>/c 0.0))) real?]{ Computes the polygamma function, @math{Ψ^m(x)}, for @math{m ≥ 0}, @math{x > 0}.} @bold{Example:} Plot of @scheme[(psi-n n x)] for @math{n = 3} on the interval (0, 5]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line (lambda (x) (psi-n 3 x))) (x-min 0.001) (x-max 5.0) (y-min 0.0) (y-max 10.0) (title "Psi-n (Polygamma) Function, Psi-n(3, x)"))] The following figure shows the resulting plot: @image["scribblings/images/psi-n.bmp"] @section{Zeta Functions} The Riemann zeta function is defined in Abramowitz and Stegun @cite["Abramowitz64"], Section 23.3. The zeta functions are defined in the @filepath{gamma.ss} file in the @schemefont{special-functions} subcollection of the science collection are are made available using the form: @;@defmodule[(planet williams/science/special-functions/gamma)] @scheme[(require (planet williams/science/special-functions/gamma))] Note that the gamma functions (Section 5.3), psi functions (Section 5.4), and the zeta functions (Section 5.5) are defined in the same module, @filepath{gamma.ss}. This is because their definitions are interdependent and PLT Scheme does not allow circular module dependencies. @subsection{Riemann Zeta Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/RiemannZetaFunction.html"]{Riemann Zeta Function} " from Wolfram MathWorld.") The Riemann zeta function is defined by the infinite sum: @image["scribblings/images/zeta-equation.gif"] @defproc[(zeta-int (n integer?)) real?]{ Computes the Reimann zeta function, @math{ζ(n)}, for integer @scheme[n], @math{n ≠ 1}.} @defproc[(zeta (s real?)) real?]{ Computes the Riemann zeta function, @math{ζ(s)}, for arbitrary @scheme[s], @math{s ≠ 1}.} @bold{Example:} Plot of @scheme[(zeta x)] on the interval [-5, 5]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line zeta) (x-min -5.0) (x-max 5.0) (y-min -5.0) (y-max 5.0) (title "Riemann Zeta Function, zeta(x)"))] The following figure shows the resulting plot: @image["scribblings/images/zeta.bmp"] @subsection{Riemann Zeta Functions Minus One} For large positive arguments, the Riemann zeta function approached one. In this region the fractional part is interesting and, therefore, we need a function to evaluate it explicitly. @defproc[(zetam1-int (n integer?)) real?]{ Computes @math{ζ(n) - 1} for integer @scheme[n], @math{n ≠ 1}.} @defproc[(zetam1 (s real?)) real?]{ Computes @math{ζ(s) - 1} for argitrary @scheme[s], @math{s ≠ 1}.} @subsection{Hutwitz Zeta Function} @(margin-note magnify @link["http://mathworld.wolfram.com/HurwitzZetaFunction.html"]{Hurwitz Zeta Function} " from Wolfram MathWorld.") The Hurwitz zeta function is defined by: @image["scribblings/images/hzeta-equation.gif"] @defproc[(hzeta (s (>/c 1.0)) (q (>/c 0.0))) real?]{ Computes the Hurwitz zeta function, @math{ζ(s, q)}, for @math{s > 1}, @math{q > 0}.} @bold{Example:} Plot of @scheme[(hzeta x 2.0)] on the interval (1, 5]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line (lambda (x) (hzeta x 2.0))) (x-min 1.001) (x-max 5.0) (y-min 0.0) (y-max 5.0) (title "Hutwitz Zeta Function, hzeta(x, 2.0)"))] The following figure shows the resulting plot: @image["scribblings/images/hzeta.bmp"] @subsection{Eta Functions} @(margin-note magnify @link["http://mathworld.wolfram.com/DirichletEtaFunction.html"]{Dirichlet Eta Function} " from Wolfram MathWorld.") The eta function is defined by: @image["scribblings/images/eta-equation.gif"] @defproc[(eta-int (n integer?)) real?]{ Computes the eta function, @math{η(n)}, for integer @scheme[n].} @defproc[(eta (s real?)) real?]{ Computes the eta function, @math{η(s)}, for arbitrary @scheme[s].} @bold{Example:} Plot of @scheme[(eta x)] on the interval [-10, 10]. @schememod[ scheme (require (planet williams/science/special-functions/gamma)) (require plot/plot) (plot (line eta) (x-min -10.0) (x-max 10.0) (y-min -5.0) (y-max 5.0) (title "Eta Function, eta(x)"))] The following figure shows the resulting plot: @image["scribblings/images/eta.bmp"]