Chapter 5

Special Functions

This chapter describes the special functions provided by the PLT Scheme Science Collection.

The functions described in this chapter are defined in the special-functions sub-collection of the science collection. The entire special-functions sub-collection can be made available using the form:

(require (planet "special-functions.ss"
                 ("williams" "science.plt")))

The individual modules in the special-functions sub-collection can also be made available as described in the sections below.

5.1 Error Functions

The error function is described in Abramowitz and Stegun, Chapter 7. The functions are defined in the error.ss file in the special-functions sub-collection of the science collection and are made available using the form:

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))

5.1.1 Error Function

erf

Function:

(erf x)

Contract:

(-> real? (real-in -1.0 1.0))

This function computes the error function:

$\mathnormal{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt.$

Example: Plot of (erf x) on the interval [ - 4, 4].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line erf)
      (x-min -4) (x-max 4)
      (y-min -1) (y-max 1)
      (title "Error Function, erf(x)"))

Figure 2 shows the resulting plot.

Figure 2: Plot of Error Function on [ - 4, 4]

5.1.2 Complementary Error Function

erfc

Function:

(erfc x)

Contract:

(-> real? (real-in 0.0 2.0))

This function computes the complementary error function:

$\mathnormal{erfc}(x) = 1 - \mathnormal{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^{2}} dt.$

Example: Plot of (erfc x) on the interval [ - 4, 4].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line erfc)
      (x-min -4) (x-max 4)
      (y-min 0) (y-max 2)
      (title "Complementary Error Function, erfc(x)"))

Figure 3 shows the resulting plot.

Figure 3: Plot of Complementary Error Function on [ - 4, 4]

5.1.3 Hazard Function

The hazard function for the normal distribution, also known as the inverse Mill's ratio, is the ratio of the probability function, P(x), to the survival function, S(x), and is defined as:

$h(x) = \frac{P(x)}{S(x)} = \sqrt{\frac{2 \pi e^{\frac{x^{2}}{2}}}{erfc({x^{\sqrt{2}}})}}.$

It decreases rapidly as x approaches - infty and asymptopes to h(x) x as x approaches + infty.

hazard

Function:

(hazard x)

Contract:

(-> real? (>=/c 0.0))

This function computes the hazard function for the normal distribution.

Example: Plot of (hazard x) on the interval [ - 5, 10].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line hazard)
      (x-min -5) (x-max 10)
      (y-min 0) (y-max 10)
      (title "Hazard Function, hazard(x)"))

Figure 4 shows the resulting plot.

Figure 4: Plot of Hazard Function on [ - 5, 10]

5.2 Exponential Integral Functions

Information on the exponential integrals can be found in Abramowitz and Stegun, Chapter 5. The functions are defined in the exponential-integral.ss file in the special-functions sub-collection of the science collection and are made available using the form:

(require (planet "exponential-integral.ss"
                 ("williams" "science.plt")
                 "special-functions"))

5.2.1 Exponential Integral

expint-E1 expint-E1-scaled

Function:

(expint-E1 x)

Function:

(expint-E1-scaled x)

Contract:

(-> real? real?)

These functions compute the exponential integral E₁(x)

$E_1(x) = \int_1^\infty \frac{e^{-tx}}{t} dt.$

Example: Plot of (expint-E1 x) on the interval [ - 4, 4].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line expint_E1)
      (x-min -4) (x-max 4)
      (y-min -10) (y-max 10)
      (title "Exponential Integral E_1"))

Figure 5 shows the resulting plot.

Figure 5: Plot of Error Function on [ - 4, 4]

expint-E2 expint-E2-scaled

Function:

(expint-E2 x)

Function:

(expint-E2-scaled x)

Contract:

(-> real? real?)

These functions compute the second-order exponential integral E₂(x)

$E_2(x) = \int_1^\infty \frac{e^{-xt}}{t^2} dt.$

Example: Plot of (expint-E2 x) on the interval [ - 4, 4].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line expint_E2)
      (x-min -4) (x-max 4)
      (y-min -10) (y-max 10)
      (title "Exponential Integral E_2"))

Figure 6 shows the resulting plot.

Figure 6: Plot of Error Function on [ - 4, 4]

5.2.2 Ei(x)

expint-Ei expint-Ei-scaled

Function:

(expint-Ei-scaled x)

Function:

(expint-Ei x)

Contract:

(-> real? real?)

These functions compute the exponential integral Ei(x).

$E_i(x) = - PV \int_{-x}^\infty \frac{e^{-t}}{t} dt.$

Example: Plot of (expint-Ei x) on the interval [ - 4, 4].

(require (planet "error.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line expint_Ei)
      (x-min -4) (x-max 4)
      (y-min -10) (y-max 10)
      (title "Exponential Integral E_1"))

Figure 7 shows the resulting plot.

Figure 7: Plot of Error Function on [ - 4, 4]

5.3 Gamma and Beta Functions

The gamma functions are defined in the gamma.ss file in the special-functions sub-collection of the science collection and are made available using the form:

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))

Note that the gamma functions (Section 5.2), psi functions (Section 5.3), and the zeta functions (Section 5.4) are defined in the same module, gamma.ss. This is because their definitions are interdependent and PLT Scheme does not allow circular module dependencies.

5.3.1 Gamma Function

The gamma function is defined by the integral:

$\Gamma(x) \equiv \int_{0}^{\infty} t^{x-1} e^{-t} dt.$

It is related to the factorial function by Gamma(n) = (n - 1)! for positive integer n. Further information on the gamma function can be found in Abramowitz & Stegun, Chapter 6.

gamma

Function:

(gamma x)

Contract:

(-> real? real?)

This function computes the gamma function, Gamma(x), subject to x not being a negative integer. The function is computed using the real Lanczos method. The maximum value of x such that Gamma(x) is not considered an overflow is given by the constant gamma-xmax and is 171.0.

Example: Plot of (gamma x) on the interval (0, 6].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line gamma)
      (x-min 0.001) (x-max 6)
      (y-min 0) (y-max 120)
      (title "Gamma Function, gamma(x)"))

Figure 8 shows the resulting plot.

Figure 8: Plot of Gamma Function on (0, 6]

Example: Plot of (gamma x) on the interval ( - 1, 0).

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line gamma)
      (x-min -0.999) (x-max -0.001)
      (y-min -120) (y-max 0)
      (title "Gamma Function, gamma(x)"))

Figure 9 shows the resulting plot.

Figure 9: Plot of Gamma Function on ( - 1, 0)

lngamma

Function:

(lngamma x)

Contract:

(-> real? real?)

This function computes the logarithm of the gamma function, log(Gamma(x)), subject to x not being a negative integer. For x < 0, the real part of log(Gamma(x)) is returned, which is equivalent to log(|Gamma(x)|). The function is computed using the real Lanczos method.

lngamma-sgn

Function:

(lngamma-sgn x)

Contract:

(-> real? (values real? (integer-in -1 1))

This function computes the logarithm of the magnitude of the gamma function and its sign, subject to 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 Gamma(x) = sgn * exp(resultlg), where resultlg and sgn are the returned values.

Example: Plot of (lngamma x) on the interval (0, 6].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line lngamma)
      (x-min 0.001) (x-max 6)
      (y-min 0) (y-max 120)
      (title "Log Gamma Function, lngamma(x)"))

Figure 10 shows the resulting plot.

Figure 10: Plot of Log Gamma Function on (0, 6]

gamma-inv

Function:

(gamma-inv x)

Contract:

(-> real? real?)

This function computes the reciprocal of the gamma function, (1/Gamma(x)), using the real Lanczos method.

5.3.2 Regulated Gamma Function

The regulated gamma function is given by

$\Gamma^*(x) = \frac{\Gamma(x)}{\sqrt{2 \pi}}x^{x-\frac{1}{2}}e^{-x}.$

gammastar

Function:

(gammastar x)

Contract:

(-> (>/c 0.0) real?)

This function computes the regulated gamma function, Gamma^*(x), for x > 0.

Example: Plot of (erf x) on the interval (0, 4].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line gammastar)
      (x-min 0.001) (x-max 4)
      (y-min 0) (y-max 10)
      (title "Regulated Gamma Function, gamma*(x)"))

Figure 11 shows the resulting plot.

Figure 11: Plot of Regulated Gamma Function on (0, 4]

5.3.3 Incomplete Gamma Function

gamma-inc-Q

Function:

(gamma-inc-Q a x)

Contract:

(-> (>/c 0.0) (>=/c 0.0) real?)

This function computes the normalized incomplete gamma function,

$Q(a,x) = 1/\Gamma(a)\int_a^\infty t^{a-1} e^{-t} dt,$

for a > 0 and x > = 0.

gamma-inc-P

Function:

(gamma-inc-Q a x)

Contract:

(-> (>/c 0.0) (>=/c 0.0) real?)

This function computes the complementary normalized incomplete gamma function,

$P(a,x) = 1/\Gamma(a)\int_0^x t^{a-1} e^{-t} dt,$

for a > 0 and x > = 0.

gamma-inc

Function:

(gamma-inc a x)

Contract:

(-> real? (>=/c 0.0) real?)

This function computes the incomplete gamma function,

$\Gamma(a,x) = \int_x^\infty t^{a-1} e^{-t} dt,$

for a real and x > = 0.

5.3.4 Factorial Function

The factorial of a positive integer n, n!, is defined as n! = n × (n - 1) × ... × 2 × 1. By definition, 0! = 1. The factorial function is related to the gamma function by n! = Gamma(n + 1).

factorial

Function:

(fact n)

Contract:

(-> natural-number? (>/c 1.0))

This function computes the factorial of n, n!.

lnfact

Function:

(lnfact n)

Contract:

(-> natural-number? (>=/c 0.0))

This function computes the logarithm of the factorial of n, log(n!). The algorithm is faster than computing ln gamma(n + 1) via lngamma for n < 170, but defers for larger n.

5.3.5 Double Factorial Function

The double factorial of n, n!!, is defined as n!! = n × (n - 2) × (n - 4) × .... By definition, - 1!! = 0!! = 1.

double-fact

Function:

(double-fact n)

Contract:

(-> natural-number? (>/c 1.0))

This function computes the double factorial of n, n!!.

lndouble-fact

Function:

(lndouble-fact n)

Contract:

(-> natural-number? (>=/c 0.0))

This function computes the logarithm of the double factorial of n, log(n!!).

5.3.6 Binomial Coefficient Function

The binomial coefficient, b choose m, is defined as:

choose

Function:

(choose n m)

Contract:

(-> natural-number? natural-number? (>/c 1.0))

This function computes the binomial coefficient for n and m, n choose m.

lnchoose

Function:

(lnchoose n m)

Contract:

(-> natural-number? natural-number? (>/c 0.0))

This function computes the logarithm of the binomial coefficient for n and m,
log (n choose m).

5.3.7 Beta Functions

beta

Function:

(beta a b)

Contract:

(-> (>/c 0.0) (>/c 0.0) real?)

This function computes the Beta function,

$B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b),$

for a > 0 and b > 0.

lnbeta

Function:

(beta a b)

Contract:

(-> (>/c 0.0) (>/c 0.0) real?)

This function computes the logorithm of the Beta function, log(B(a,b)), for a > 0 and b > 0.

5.4 Psi Functions

The psi functions are defined in gamma.ss in the special-functions sub-collection of the science collection and are made available using the form:

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))

Note that the gamma functions (Section 5.2), psi functions (Section 5.3), and zeta functions (Section 5.4) are defined in the same module, gamma.ss. This is because their definitions are interdependent and PLT Scheme does not allow circular module dependencies.

5.4.1 Psi (Digamma) Functions

psi-int

Function:

(psi-int n)

Contract:

(-> (integer-in 1 +inf.0) real?)

This function computes the digamma function, psi(n), for positive integer n. The digamma function is also called the Psi function.

psi

Function:

(psi x)

Contract:

(-> real? real?)

This function returns the digamma function, psi(x), for general x, x != 0.

psi-1piy

Function:

(psi-1piy y)

Contract:

(-> real? real?)

This function computes the real part of the digamma function on the line 1 + i y, Re[psi(1 + i y)].

Example: Plot of (psi x) on the interval (0, 5].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line psi)
      (x-min 0.001) (x-max 5)
      (y-min -5) (y-max 5)
      (title "Psi (Digamma) Function, psi(x)"))

Figure 12 shows the resulting plot.

Figure 12: Plot of Psi (Digamma) Function on (0, 5]

5.4.2 Psi-1 (Trigamma) Functions

psi-1-int

Function:

(psi-1-int n)

Contract:

(-> (integer-in 1 +inf.0) real?)

This function computes the trigamma function psi'(n) for positive integer n.

psi-1

Function:

(psi-1 x)

Contract:

(-> real? real?)

This function computes the trigamma function psi'(x) for general x.

Example: Plot of (psi-1 x) on the interval (0, 5].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line psi-1)
      (x-min 0.005) (x-max 5)
      (y-min 0) (y-max 5)
      (title "Psi-1 (Trigamma) Function, psi-1(x)"))

Figure 13 shows the resulting plot.

Figure 13: Plot of Psi-1 (Trigamma) Function on (0, 5]

5.4.3 Psi-n (Polygamma) Function

psi-n

Function:

(psi-n n x)

Contract:

(-> natural-number? (>/c 0.0) real?)

This function computes the polygamma function psi^m(x) for m > = 0, x > 0.

Example: Plot of (psi-n n x) on the interval (0, 5].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line (lambda (x) (psi-n 3 x)))
      (x-min 0.05) (x-max 5)
      (y-min 0) (y-max 10)
      (title "Psi-n (Polygamma) Function, psi-n(3, x)"))

Figure 14 shows the resulting plot.

Figure 14: Plot of Psi-n (Polygamma), n = 3, Function on (0, 5]

5.5 Zeta Functions

The Riemann zeta function is defined in Abramowitz and Stegun, Section 23.2. The zeta functions are defined in the gamma.ss file in the special-functions sub-collection of the science collection and are made available using the form:

(require (lib "gamma.ss" ("williams" "science.plt")
              "special-functions"))

5.5.1 Riemann Zeta Functions

The Riemann zeta function is defined by the infinite sum

$\zeta(s) = \sum_{k=1}^\infty k^{-s}.$

zeta-int

Function:

(zeta-int n)

Contract:

(-> integer? real?)

This function computes the Reimann zeta function zeta(n) for integer n, n != 1.

zeta

Function:

(zeta s)

Contract:

(-> real? real?)

This function computes the Reimann zera function zeta(s) for arbitraty s, s != 1.

Example: Plot of (zeta x) on the interval [ - 5, 5].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line zeta)
      (x-min -5) (x-max 5)
      (y-min -5) (y-max 5)
      (title "Reimann Zeta Function, zeta(x)"))

Figure 15 shows the resulting plot.

Figure 15: Plot of Reimann Zeta Function on [ - 5, 5]

5.5.2 Riemann Zeta Functions Minus One

For large positive arguments, the Reimann zeta function approaches one. In this region the fractional part is interesting amd, therefore, we need a function to evaluate it explicitly.

zetam1-int

Function:

(zetam1-int n)

Contract:

(-> integer? real?)

This function computes zeta(n) - 1 for integer n, n != 1.

zetam1

Function:

(zetam1 s)

Contract:

(-> real? real?)

This function computes zeta(n) - 1 for arbitrary s, s != 1.

5.5.3 Hurwitz Zeta Function

The Hurwitx zeta function is defined by

$\zeta(s,q) = \sum_{k=0}^\infty (k+q)^{-s}.$

hzeta

Function:

(hzeta s q)

Contract:

(-> (>/c 1.0) (>/c 0.0) real?)

This function computes the Hurwitz zeta function zeta(s,q) for s > 1, q > 0.

Example: Plot of (hzeta x) on the interval (1, 5].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line (lambda (x) (hzeta x 2.0)))
      (x-min 1.001) (x-max 5)
      (y-min 0) (y-max 5)
      (title "Hurwitz Zeta Function, hzeta(x, 2.0)"))

Figure 16 shows the resulting plot.

Figure 16: Plot of Hurwitz Zeta Function, q = 2.0, on (1, 5]

5.5.4 Eta Functions

The eta function is defined by

$\eta(s) = (1-2^{1-s})\zeta(s).$

eta-int

Function:

(eta-int n)

Contract:

(-> integer? real?)

This function computes the eta function eta(n) for integer n.

eta

Function:

(eta s)

Contract:

(-> real? real?)

This function computes the eta function eta(s) for arbitrary s.

Example: Plot of (eta x) on the interval [ - 10, 10].

(require (planet "gamma.ss" ("williams" "science.plt")
                 "special-functions"))
(require (lib "plot.ss" "plot"))

(plot (line eta)
      (x-min -10) (x-max 10)
      (y-min -5) (y-max 5)
      (title "Eta Function, eta(x)"))

Figure 17 shows the resulting plot.

Figure 17: Plot of Eta Function on [ - 10, 10]