Breakthroughs in Special Functions and Transforms#

mathematicskit.special_functions hub · History · Examples · API reference

“I have had my results for a long time, but I do not yet know how I am to arrive at them.” – attributed to Carl Friedrich Gauss

The “special functions” – gamma, Bessel, the orthogonal polynomial families – are special because they keep recurring, unbidden, as exact solutions of a huge range of unrelated differential equations, from vibrating drumheads to quantum wavefunctions. This chronology traces the ideas behind mathematicskit.special_functions, from Leonhard Euler’s extension of the factorial to the fast Fourier transform, an algorithm rediscovered from a much older idea that made digital signal processing computationally feasible, and on to the transforms, filters, and wavelets that signal processing is built from.

1729 – Euler’s Gamma Function #

Christian Goldbach asked how the factorial could be extended to non-integer values. Leonhard Euler answered in an October 1729 letter with an infinite-product formula for such an interpolation, and gave the integral form early the next year. Adrien-Marie Legendre later named the function \(\Gamma\). It satisfies \(\Gamma(n) = (n-1)!\) at positive integers and extends smoothly everywhere else, with poles at zero and the negative integers. The closely related beta integral \(B(a,b)\), which Euler also studied and Jacques Binet named in 1839, can be written directly in terms of \(\Gamma\).

\[\Gamma(x) = \int_0^\infty t^{x-1}e^{-t}\,dt, \qquad B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\]

Implementation: mathematicskit.special_functions.systems.gamma_beta.gamma_function() and beta_function() wrap scipy.special.gamma()/beta() directly. log_gamma_function() wraps the numerically stable scipy.special.gammaln() for arguments large enough that \(\Gamma\) itself would overflow.

References: L. Euler, letter to C. Goldbach, 13 October 1729, in P.-H. Fuss, ed., Correspondance mathématique et physique de quelques célèbres géomètres du XVIIIème siècle, vol. 1 (St. Petersburg, 1843), 56-60.

The gamma function extends the factorial

1730 – Stirling’s Formula #

Abraham de Moivre, studying the binomial distribution, found in 1730 that \(n!\) grows like \(C\,n^{n+1/2}e^{-n}\) for some constant \(C\). In the same year James Stirling’s Methodus Differentialis identified the constant as \(\sqrt{2\pi}\):

\[n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n .\]

The approximation is already good to 1% at \(n = 10\), and its relative error shrinks like \(1/(12n)\). Adding correction terms built from the Bernoulli numbers gives Stirling’s series for \(\ln\Gamma(x)\). That series never converges. For each fixed \(x\) its terms first shrink and then grow without bound. It was one of the first asymptotic series: truncated at the right point, it is extremely accurate, and it is still how \(\ln\Gamma\) is computed for large arguments.

Implementation: mathematicskit.special_functions.systems.gamma_beta.stirling_factorial() evaluates the leading formula, and stirling_log_gamma() sums the series with a chosen number of Bernoulli terms from scipy.special.bernoulli(). The tests check the \(1/(12n)\) error and show that at small \(x\) the error eventually grows again as terms are added.

References: J. Stirling, Methodus Differentialis: sive Tractatus de Summatione et Interpolatione Serierum Infinitarum (London, 1730), Proposition 28; A. de Moivre, Miscellanea Analytica de Seriebus et Quadraturis (London, 1730); NIST Digital Library of Mathematical Functions, Sec. 5.11.

Stirling’s formula and its asymptotic series

Stirling's formula and its asymptotic series

1758-1996 – The Lambert W Function #

In 1758 Johann Heinrich Lambert gave a series solution of the trinomial equation \(x = q + x^m\), and in 1783 Leonhard Euler recast it in a more symmetric form. Hidden inside their work is the function \(W(z)\) that solves

\[W(z)\,e^{W(z)} = z .\]

For two centuries it had no name, even though it kept turning up in delay equations, enzyme kinetics, and combinatorial counting of trees. In 1996 Robert Corless, Gaston Gonnet, David Hare, David Jeffrey, and Donald Knuth collected these scattered appearances into one paper, fixed the notation \(W\), and defined its infinitely many complex branches. For real \(-1/e \le z < 0\) there are two real solutions, the principal branch \(W_0\) and the lower branch \(W_{-1}\), which meet at \(z = -1/e\). Many equations that mix exponentials and powers can be solved with \(W\). For example, \(x^x = y\) has the solution \(x = \ln y / W(\ln y)\).

Implementation: mathematicskit.special_functions.systems.lambert_w.lambert_w() wraps scipy.special.lambertw() and returns a real result on the real branches. The tests check \(we^w = z\), the omega constant \(W(1) = e^{-W(1)}\), and both branches at \(z = -\tfrac12\ln 2\).

References: J. H. Lambert, “Observationes variae in mathesin puram,” Acta Helvetica 3 (1758), 128-168; L. Euler, “De serie Lambertina plurimisque eius insignibus proprietatibus,” Acta Academiae Scientiarum Imperialis Petropolitanae 1779:2 (1783), 29-51; R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, “On the Lambert W function,” Advances in Computational Mathematics 5 (1996), 329-359.

The two real branches of the Lambert W function

1782-1854 – Legendre and Chebyshev Polynomials #

Adrien-Marie Legendre’s 1782 work on the gravitational attraction of spheroids introduced the polynomial family that now bears his name. The polynomials solve Laplace’s equation in spherical coordinates and are orthogonal on \([-1,1]\) with the simplest possible weight, 1. Pafnuty Chebyshev’s 1854 memoir on least-deviation approximation (see the numerical-analysis chronology) introduced a second family, orthogonal with respect to the weight \(1/\sqrt{1-x^2}\) and given in closed form by \(T_n(x) = \cos(n\arccos x)\). Its minimal-oscillation property is exactly what makes Chebyshev nodes the cure for Runge’s phenomenon.

Implementation: mathematicskit.special_functions.systems.orthogonal_polynomials.legendre_polynomial() and chebyshev_polynomial() wrap numpy.polynomial.legendre.legval()/chebval() directly, and their orthogonality is verified numerically with mathematicskit.special_functions.utils.orthogonality.inner_product().

References: A.-M. Legendre, “Recherches sur l’attraction des sphéroïdes homogènes,” Mémoires de Mathématique et de Physique, présentés à l’Académie Royale des Sciences 10 (1785), 411-434 (presented 1782).

Legendre and Chebyshev polynomials on [-1, 1]

1799 – Gauss’s Arithmetic-Geometric Mean #

Start with two positive numbers and repeatedly replace them by their arithmetic and geometric means:

\[a_{n+1} = \frac{a_n + b_n}{2}, \qquad b_{n+1} = \sqrt{a_n b_n}.\]

The two sequences converge to a common limit, the arithmetic-geometric mean \(\mathrm{AGM}(a, b)\), and the number of correct digits roughly doubles at every step. On 30 May 1799 the 21-year-old Carl Friedrich Gauss noted in his diary that \(1/\mathrm{AGM}(1,\sqrt2)\) agreed to eleven decimals with a constant from the arc length of the lemniscate, and predicted that this would open “an entirely new field of analysis.” The explanation is an identity for the complete elliptic integral of the first kind:

\[K(m) = \int_0^{\pi/2}\frac{d\theta}{\sqrt{1 - m\sin^2\theta}} = \frac{\pi}{2\,\mathrm{AGM}\bigl(1, \sqrt{1-m}\bigr)} .\]

Elliptic integrals had been studied since the early eighteenth century because they give the arc length of an ellipse and the period of a pendulum, but no one could evaluate them quickly. Gauss’s identity turns them into a few square roots. The same quadratic convergence drives the Brent-Salamin algorithm of 1976, which has been used to compute \(\pi\) to billions of digits.

Implementation: mathematicskit.special_functions.systems.elliptic.arithmetic_geometric_mean(), complete_elliptic_integral_first_kind(), and complete_elliptic_integral_second_kind() wrap scipy.special.agm(), ellipk(), and ellipe(). The tests check Gauss’s identity, compare both integrals with direct quadrature, and verify Legendre’s relation \(EK' + E'K - KK' = \pi/2\).

References: C. F. Gauss, Werke, vol. 3 (Göttingen, 1866), 361-403; D. A. Cox, “The Arithmetic-Geometric Mean of Gauss,” L’Enseignement Mathématique 30 (1984), 275-330; J. M. Borwein and P. B. Borwein, Pi and the AGM (New York: Wiley, 1987).

Gauss’s arithmetic-geometric mean and the elliptic integral K

Gauss's arithmetic-geometric mean and the elliptic integral K

1810-1879 – Hermite and Laguerre Polynomials #

Pierre-Simon Laplace in 1810, and Pafnuty Chebyshev more systematically in 1859, studied polynomials orthogonal with respect to the Gaussian weight \(e^{-x^2}\). Charles Hermite’s 1864 memoir gave them their definitive treatment and the name they carry today. Edmond Laguerre’s 1879 paper introduced the complementary family, orthogonal on the half-line \([0,\infty)\) with weight \(e^{-x}\). Nearly half a century later, both families turned out to be exactly the eigenfunctions quantum mechanics needed: Hermite polynomials for the quantum harmonic oscillator, and Laguerre polynomials for the radial part of the hydrogen atom’s wavefunctions.

Implementation: mathematicskit.special_functions.systems.orthogonal_polynomials.hermite_polynomial() and laguerre_polynomial() wrap numpy.polynomial.hermite.hermval()/lagval(). Their orthogonality against each family’s weight is verified numerically with scipy.integrate.quad() over an infinite domain.

References: C. Hermite, “Sur un nouveau développement en série des fonctions,” Comptes Rendus de l’Académie des Sciences 58 (1864), 93-100, 266-273; E. Laguerre, “Sur l’intégrale \(\int_x^\infty e^{-x}x^{-1}dx\),” Bulletin de la Société Mathématique de France 7 (1879), 72-81.

Hermite and Laguerre polynomials on infinite domains

1812 – Gauss’s Hypergeometric Series #

Leonhard Euler had already studied the series

\[{}_2F_1(a, b; c; z) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n}\,\frac{z^n}{n!}, \qquad (q)_n = q(q+1)\cdots(q+n-1),\]

and the second-order differential equation it satisfies. Carl Friedrich Gauss’s memoir, presented to the Göttingen Society in January 1812, made it a unified object of study. He showed that the one family contains the logarithm, the inverse trigonometric functions, and many other functions of analysis as special cases, for example \(\ln(1+z) = z\,{}_2F_1(1,1;2;-z)\). He gave one of the first rigorous convergence proofs for an infinite series, and he found the value at \(z = 1\):

\[{}_2F_1(a, b; c; 1) = \frac{\Gamma(c)\,\Gamma(c-a-b)}{\Gamma(c-a)\,\Gamma(c-b)}, \qquad \operatorname{Re}(c - a - b) > 0 .\]

When \(a\) is a negative integer the series stops after finitely many terms, and every classical orthogonal polynomial family on a finite interval is a special case. Ernst Kummer and Bernhard Riemann later studied the function through its differential equation, which has exactly three singular points, at \(0\), \(1\), and \(\infty\).

Implementation: mathematicskit.special_functions.systems.hypergeometric.hypergeometric_2f1() wraps scipy.special.hyp2f1(). The tests recover \(\ln(1+z)\), \((1+z)^{2.5}\), \(\arcsin z\), and the Chebyshev polynomial \(T_4\) from it, and check Gauss’s summation theorem against scipy.special.gamma().

References: C. F. Gauss, “Disquisitiones generales circa seriem infinitam \(1 + \frac{\alpha\beta}{1\cdot\gamma}x + \cdots\),” Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores 2 (1813); NIST Digital Library of Mathematical Functions, Ch. 15.

Gauss’s hypergeometric function contains the elementary functions

Gauss's hypergeometric function contains the elementary functions

1812 – Laplace’s Transform #

Leonhard Euler had already used integrals of the form \(\int e^{-st}f(t)\,dt\) in the 1740s and 1760s to solve differential equations. Pierre-Simon Laplace, in the 1780s and then in his 1812 Théorie analytique des probabilités, used them systematically as generating functions that turn a problem about \(f(t)\) into one about a new function \(F(s)\):

\[F(s) = \int_0^\infty f(t)\,e^{-st}\,dt, \qquad \mathcal{L}\{f'\}(s) = sF(s) - f(0).\]

The second identity is the reason the transform matters: derivatives become multiplication by \(s\), so a linear differential equation with constant coefficients becomes an algebraic equation. Oliver Heaviside’s operational calculus of the 1890s exploited exactly this for electrical circuits, without proofs. Thomas Bromwich’s 1916 inversion integral along a vertical line in the complex \(s\)-plane, and Gustav Doetsch’s textbooks of the 1930s, gave the method a rigorous footing, and in that form it became the working language of control engineering.

Implementation: mathematicskit.special_functions.systems.laplace_transform.laplace_transform() evaluates the defining integral with scipy.integrate.quad(), for real or complex \(s\). The tests check it against the table pairs for \(e^{-at}\), \(t^2\), \(te^{-t}\), \(1-e^{-t}\), and \(\sin t\).

References: P.-S. Laplace, Théorie analytique des probabilités (Paris: Courcier, 1812); T. J. I’A. Bromwich, “Normal coordinates in dynamical systems,” Proceedings of the London Mathematical Society (2) 15 (1916), 401-448; G. Doetsch, Theorie und Anwendung der Laplace-Transformation (Berlin: Springer, 1937); M. A. B. Deakin, “The development of the Laplace transform, 1737-1937,” Archive for History of Exact Sciences 25 (1981), 343-390.

Laplace’s transform: from functions of t to functions of s

Laplace's transform: from functions of t to functions of s

1818 – Fresnel Integrals and the Diffraction of Light #

Augustin-Jean Fresnel’s 1818 prize memoir for the French Academy explained diffraction by treating light as a wave. The light reaching a point past an obstacle is the sum of wavelets from every point of the unobstructed wavefront. For a straight edge that sum reduces to the integrals

\[C(t) = \int_0^t \cos\!\left(\frac{\pi\tau^2}{2}\right)d\tau, \qquad S(t) = \int_0^t \sin\!\left(\frac{\pi\tau^2}{2}\right)d\tau .\]

Fresnel tabulated them by hand, and his predicted fringe patterns matched experiment. The judge Siméon Denis Poisson objected that the theory implied a bright spot at the centre of a circular shadow. François Arago looked for the spot and found it, and the prize went to Fresnel. In 1874 Alfred Cornu plotted \((C(t), S(t))\) as a curve, a spiral that winds into the points \(\pm(\tfrac12, \tfrac12)\). Its curvature grows in proportion to arc length. For that reason the same curve, often called the Euler or clothoid spiral, is used as the transition curve in road and railway design.

Implementation: mathematicskit.special_functions.systems.error_functions.fresnel_integrals() wraps scipy.special.fresnel() and returns a FresnelResult. The tests check the limits \(\pm\tfrac12\), the derivatives, and that the spiral has unit speed and curvature \(\pi t\).

References: A. Fresnel, “Mémoire sur la diffraction de la lumière,” Mémoires de l’Académie Royale des Sciences de l’Institut de France 5 (1826), 339-475; A. Cornu, “Méthode nouvelle pour la discussion des problèmes de diffraction dans le cas d’une onde cylindrique,” Journal de Physique Théorique et Appliquée 3 (1874), 5-15, 44-52.

Fresnel integrals and the Cornu spiral

1822-1965 – Fourier’s Theorem and the Fast Fourier Transform #

Joseph Fourier’s 1822 treatise on heat conduction claimed, controversially at the time, that essentially any periodic function can be written as a sum of sines and cosines. The claim, first presented to the Paris Academy in 1807, was made rigorous only gradually over the following century, but it became the conceptual foundation of the discrete Fourier transform once digital computers made the transform something to actually compute. James Cooley and John Tukey’s 1965 paper showed how to compute a length-\(N\) discrete Fourier transform recursively from two length-\(N/2\) transforms, cutting the cost from \(O(N^2)\) to \(O(N\log N)\). They had largely independently rediscovered an algorithm that Carl Friedrich Gauss sketched in an unpublished 1805 manuscript on interpolating asteroid orbits. For the signal lengths that radar and telecommunications needed, the speedup was so large that the algorithm is routinely credited with making digital signal processing practical.

\[X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}\]

Implementation: mathematicskit.special_functions.systems.fourier_transform.fft_numpy() wraps numpy.fft.fft() as the primary API. dft_naive() and fft_radix2() hand-roll the naive \(O(N^2)\) sum and the Cooley-Tukey radix-2 recursion. compare_fft_methods() times all three against each other to make the \(O(N^2)\) versus \(O(N\log N)\) gap visible. This is the one submodule in this domain where hand-rolling, rather than a direct library call, is the point.

References: J. Fourier, Théorie analytique de la chaleur (Paris: Firmin Didot, 1822); J. W. Cooley and J. W. Tukey, “An Algorithm for the Machine Calculation of Complex Fourier Series,” Mathematics of Computation 19(90) (1965), 297-301; M. T. Heideman, D. H. Johnson, and C. S. Burrus, “Gauss and the History of the Fast Fourier Transform,” Archive for History of Exact Sciences 34(3) (1985), 265-277 (documenting Gauss’s unpublished 1805 anticipation).

O(n^2) vs. O(n log n): naive DFT vs. radix-2 FFT vs. numpy.fft

1824 – Bessel Functions #

In the 1730s, Daniel Bernoulli’s study of the oscillations of a hanging chain had already produced solutions of what is now recognized as a special case of Bessel’s equation. Friedrich Wilhelm Bessel’s 1824 memoir on planetary perturbations gave the family its first systematic treatment and tabulation. His motivation was the periodic problem of Kepler’s equation for planetary orbits. \(J_\nu\), regular at the origin, and \(Y_\nu\), singular there and introduced later by Carl Neumann and Heinrich Weber, are two linearly independent solutions of the same second-order equation. They appear throughout physics and engineering wherever a problem has cylindrical symmetry: vibrating drumheads, waveguides, heat conduction in a cylinder.

\[x^2 y'' + xy' + (x^2 - \nu^2)y = 0\]

Implementation: mathematicskit.special_functions.systems.bessel.bessel_first_kind() and bessel_second_kind() wrap scipy.special.jv()/yv() directly. This domain’s tests check them against the defining differential equation itself, using finite differences.

References: F. W. Bessel, “Untersuchung des Theils der planetarischen Störungen, welcher aus der Bewegung der Sonne entsteht,” Abhandlungen der Königlichen Akademie der Wissenschaften zu Berlin (1824), 1-52.

Bessel functions of the first and second kind

1829 – Jacobi’s Elliptic Functions #

Niels Henrik Abel and Carl Gustav Jacobi, working independently in 1827-1829, had the idea of inverting the elliptic integral. The sine can be defined as the inverse of \(\int_0^x dt/\sqrt{1-t^2}\). In the same way, if

\[u = \int_0^{\varphi}\frac{d\theta}{\sqrt{1 - m\sin^2\theta}},\]

then Jacobi defined \(\operatorname{sn} u = \sin\varphi\), \(\operatorname{cn} u = \cos\varphi\), and \(\operatorname{dn} u = \sqrt{1 - m\sin^2\varphi}\). The new functions are doubly periodic: they repeat along two independent directions in the complex plane. On the real line \(\operatorname{sn}\) and \(\operatorname{cn}\) have period \(4K(m)\). At \(m = 0\) they reduce to \(\sin\) and \(\cos\), and at \(m = 1\) to hyperbolic functions. Jacobi’s 1829 Fundamenta nova developed the theory in full. The functions give exact solutions of the large-amplitude pendulum, of the rigid body spinning freely, and of the soliton waves of the Korteweg-de Vries equation.

Implementation: mathematicskit.special_functions.systems.elliptic.jacobi_elliptic_functions() wraps scipy.special.ellipj() and returns a JacobiEllipticResult. The tests check \(\operatorname{sn}^2 + \operatorname{cn}^2 = 1\), \(\operatorname{dn}^2 + m\operatorname{sn}^2 = 1\), the trigonometric and hyperbolic limits, and the period \(4K(m)\).

References: N. H. Abel, “Recherches sur les fonctions elliptiques,” Journal für die reine und angewandte Mathematik 2 (1827), 101-181; C. G. J. Jacobi, Fundamenta nova theoriae functionum ellipticarum (Königsberg: Bornträger, 1829).

Jacobi elliptic functions and the pendulum

1837 – Kummer’s Confluent Hypergeometric Function #

Ernst Kummer’s work of 1836-1837 studied Gauss’s hypergeometric series through its differential equation. In the same work he introduced the simpler series

\[M(a, b, z) = {}_1F_1(a; b; z) = \sum_{n=0}^{\infty}\frac{(a)_n}{(b)_n}\frac{z^n}{n!},\]

which solves \(zw'' + (b - z)w' - aw = 0\). It arises from \({}_2F_1\) when two of the three singular points are merged, or made confluent: \({}_1F_1(a; b; z) = \lim_{c\to\infty}{}_2F_1(a, c; b; z/c)\). Kummer proved the transformation \(M(a, b, z) = e^z M(b - a, b, -z)\). The confluent function contains the exponential, the error function, the incomplete gamma function, and the Laguerre and Hermite polynomials. It later gave the exact solution for the hydrogen atom and for Coulomb scattering. Edmund Whittaker’s 1904 normalization, the Whittaker functions, is the standard form in quantum mechanics.

Implementation: mathematicskit.special_functions.systems.hypergeometric.confluent_hypergeometric_1f1() wraps scipy.special.hyp1f1(). The tests check \({}_1F_1(a; a; z) = e^z\), Kummer’s transformation, the error function as a special case, the differential equation, and the \(O(1/c)\) rate of the confluent limit.

References: E. E. Kummer, “Über die hypergeometrische Reihe,” Journal für die reine und angewandte Mathematik 15 (1836), 39-83, 127-172; E. E. Kummer, “De integralibus quibusdam definitis et seriebus infinitis,” Journal für die reine und angewandte Mathematik 17 (1837), 228-242; E. T. Whittaker, “An expression of certain known functions as generalized hypergeometric functions,” Bulletin of the American Mathematical Society 10 (1904), 125-134.

Kummer’s confluent hypergeometric function

Kummer's confluent hypergeometric function

1838 – The Airy Function #

The Astronomer Royal George Biddell Airy wanted to explain the supernumerary bows, the faint extra bands just inside a rainbow. Near a caustic, where light rays bunch together, simple ray optics predicts infinite brightness. Airy’s 1838 wave calculation replaced that infinity with the integral

\[\operatorname{Ai}(x) = \frac{1}{\pi}\int_0^\infty \cos\!\left(\frac{t^3}{3} + xt\right)dt,\]

which solves the simplest equation whose solutions change character, \(y'' = xy\). For \(x < 0\) it oscillates, and for \(x > 0\) it decays exponentially. Airy computed the integral numerically by hand. George Gabriel Stokes’s 1850 asymptotic expansions made it practical to evaluate, and in studying them Stokes found the Stokes phenomenon. A second solution \(\operatorname{Bi}\), which grows for \(x > 0\), completes the pair. Every smooth turning point looks locally like Airy’s equation: in quantum mechanics this is the point where a particle’s energy equals the potential. For this reason the Airy functions are the connecting pieces of the WKB approximation.

Implementation: mathematicskit.special_functions.systems.airy.airy_functions() wraps scipy.special.airy() and returns an AiryResult. The tests check the exact values at \(x = 0\) in terms of \(\Gamma(1/3)\) and \(\Gamma(2/3)\), the differential equation, the Wronskian \(\operatorname{Ai}\operatorname{Bi}' - \operatorname{Ai}'\operatorname{Bi} = 1/\pi\), and the decay rate for large \(x\).

References: G. B. Airy, “On the intensity of light in the neighbourhood of a caustic,” Transactions of the Cambridge Philosophical Society 6 (1838), 379-402; G. G. Stokes, “On the numerical calculation of a class of definite integrals and infinite series,” Transactions of the Cambridge Philosophical Society 9 (1850), 166-187.

The Airy functions Ai and Bi

1859 – Riemann’s Zeta Function #

In 1735 Leonhard Euler solved the Basel problem, showing that \(\sum 1/n^2 = \pi^2/6\). Two years later he found the link between such sums and the primes:

\[\zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^s} = \prod_{p\ \text{prime}}\frac{1}{1 - p^{-s}}, \qquad s > 1,\]

a restatement of unique prime factorization. Bernhard Riemann’s only paper on number theory, in 1859, treated \(\zeta\) as a function of a complex variable. He extended it to the whole plane except \(s = 1\) and proved the functional equation relating \(\zeta(s)\) to \(\zeta(1-s)\). He then gave an explicit formula for the number of primes below \(x\) in terms of the zeros of \(\zeta\). The extension gives finite values to divergent sums, such as \(\zeta(-1) = -1/12\), and has “trivial” zeros at \(s = -2, -4, \ldots\). Riemann conjectured that all the other zeros lie on the line \(\operatorname{Re} s = \tfrac12\). That conjecture, the Riemann hypothesis, is still unproved.

Implementation: mathematicskit.special_functions.systems.zeta.riemann_zeta() wraps scipy.special.zeta() on the real line, including its analytic continuation to \(s < 1\). euler_product() hand-rolls the truncated product over primes from a sieve of Eratosthenes. The tests check \(\zeta(2)\), \(\zeta(4)\), and \(\zeta(6)\), the values \(\zeta(0) = -\tfrac12\) and \(\zeta(-1) = -\tfrac1{12}\), the trivial zeros, the functional equation, and the convergence of the product.

References: L. Euler, “De summis serierum reciprocarum,” Commentarii Academiae Scientiarum Petropolitanae 7 (1740), 123-134; L. Euler, “Variae observationes circa series infinitas,” Commentarii Academiae Scientiarum Petropolitanae 9 (1744), 160-188; B. Riemann, “Ueber die Anzahl der Primzahlen unter einer gegebenen Grösse,” Monatsberichte der Königlichen Preußischen Akademie der Wissenschaften zu Berlin (1859), 671-680.

The Riemann zeta function and Euler’s product

The Riemann zeta function and Euler's product

1868 – Mathieu Functions #

Émile Mathieu’s 1868 study of a vibrating drumhead of elliptical shape separated the wave equation in elliptic coordinates. It produced an equation with a periodic coefficient:

\[y'' + (a - 2q\cos 2x)\,y = 0 .\]

Unlike Bessel’s equation, this one has periodic solutions only for special characteristic values \(a = a_m(q)\) (even solutions \(\operatorname{ce}_m\)) and \(a = b_m(q)\) (odd solutions \(\operatorname{se}_m\)). As \(q \to 0\) these reduce to \(m^2\), \(\cos mx\), and \(\sin mx\). Gaston Floquet’s 1883 theory of equations with periodic coefficients and George William Hill’s 1886 work on the motion of the Moon placed Mathieu’s equation in a larger setting. Its stability chart, which shows the values of \((a, q)\) for which solutions stay bounded, explains parametric resonance: a swing pumped at twice its natural frequency. The same chart is the design principle of the Paul ion trap, for which Wolfgang Paul shared the 1989 Nobel Prize in Physics.

Implementation: mathematicskit.special_functions.systems.mathieu.mathieu_characteristic_a(), mathieu_characteristic_b(), mathieu_even(), and mathieu_odd() wrap scipy.special.mathieu_a(), mathieu_b(), mathieu_cem(), and mathieu_sem(), taking angles in radians rather than scipy’s degrees. The tests check the \(q = 0\) limit, the small-\(q\) expansion of \(a_1\), the interlacing of characteristic values, and the differential equation.

References: É. Mathieu, “Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique,” Journal de Mathématiques Pures et Appliquées (2) 13 (1868), 137-203; G. Floquet, “Sur les équations différentielles linéaires à coefficients périodiques,” Annales Scientifiques de l’École Normale Supérieure (2) 12 (1883), 47-88.

Mathieu functions and their characteristic values

1871 – The Error Function #

The integral of the Gaussian \(e^{-t^2}\) appears wherever the normal distribution does: in Abraham de Moivre’s 1733 approximation of the binomial, in Pierre-Simon Laplace’s probability theory, and in Carl Friedrich Gauss’s 1809 theory of observational errors. It has no elementary antiderivative, so it had to be tabulated. Christian Kramp did so in 1799 for work on atmospheric refraction. In 1871 James Whitbread Lee Glaisher named the function, coining the notation “erf” (for “error function”) and “erfc” for its complement. In today’s normalization

\[\operatorname{erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt, \qquad \operatorname{erfc}(x) = 1 - \operatorname{erf}(x),\]

so that the chance that a normal variable lies within \(k\) standard deviations of its mean is \(\operatorname{erf}(k/\sqrt2)\). That formula gives the familiar 68-95-99.7 rule. The complementary function is not just a convenience. For \(x\) beyond about 6, \(1 - \operatorname{erf}(x)\) rounds to exactly zero in double precision, while \(\operatorname{erfc}(x)\) is still accurate.

Implementation: mathematicskit.special_functions.systems.error_functions.error_function() and complementary_error_function() wrap scipy.special.erf()/erfc(). The tests check oddness, the defining integral, the 68-95-99.7 rule, agreement with scipy.stats.norm() cumulative probabilities, and the asymptotic tail of \(\operatorname{erfc}\).

References: J. W. L. Glaisher, “On a class of definite integrals,” Philosophical Magazine (4) 42 (1871), 294-302; C. Kramp, Analyse des réfractions astronomiques et terrestres (Strasbourg, 1799); C. F. Gauss, Theoria motus corporum coelestium in sectionibus conicis solem ambientium (Hamburg, 1809).

The error function and the normal distribution

1910 – Haar’s Orthonormal System #

David Hilbert asked whether every continuous function could be expanded in some orthonormal system whose series always converged uniformly, something Fourier series famously fail to do. His student Alfréd Haar answered in his 1909 dissertation, published in 1910, with a system of step functions: a constant, then functions equal to \(+1\) on one half of a dyadic interval and \(-1\) on the other, rescaled and shifted to every scale.

\[\begin{split}\psi(x) = \begin{cases} 1, & 0 \le x < \tfrac12 \\ -1, & \tfrac12 \le x < 1 \end{cases}, \qquad \psi_{j,k}(x) = 2^{j/2}\psi(2^j x - k).\end{split}\]

Its coefficients are local averages and differences, so a jump in the signal affects only one coefficient per scale. Seventy years later the Haar system was recognized as the first and simplest wavelet basis.

Implementation: mathematicskit.special_functions.systems.wavelets.discrete_wavelet_transform() with wavelet="haar" computes pairwise sums and differences over \(\sqrt2\) level by level, and inverse_discrete_wavelet_transform() inverts it. Neither numpy nor scipy (since scipy 1.15) has wavelet routines, so both are hand-rolled. The tests check the pairwise formula, energy preservation, and exact reconstruction.

References: A. Haar, “Zur Theorie der orthogonalen Funktionensysteme,” Mathematische Annalen 69 (1910), 331-371.

Haar’s wavelet: averages and differences at every scale

Haar's wavelet: averages and differences at every scale

1930 – Butterworth’s Maximally Flat Filter #

Radio engineers of the 1920s built filters by trial and error, and the passbands rippled. Stephen Butterworth, a British physicist, asked for the filter whose gain is as flat as possible at zero frequency, and found that the squared magnitude

\[|H(i\Omega)|^2 = \frac{1}{1 + (\Omega/\Omega_c)^{2N}}\]

has its first \(2N - 1\) derivatives zero at \(\Omega = 0\). Its poles are equally spaced on a circle of radius \(\Omega_c\) in the left half-plane, it passes through half power at the cutoff for every order, and it falls off at \(20N\) dB per decade. Butterworth also showed how to build it from identical amplifier stages. The Chebyshev and elliptic (Cauer) filters of the following decades trade some ripple for a sharper cutoff. For digital filters, the bilinear transform \(s = 2(1 - z^{-1})/(1 + z^{-1})\) maps the analog design to discrete time.

Implementation: mathematicskit.special_functions.systems.filters.butterworth_filter() and chebyshev1_filter() wrap scipy.signal.butter()/cheby1(); frequency_response() wraps scipy.signal.freqz() and apply_filter() wraps scipy.signal.lfilter()/filtfilt(). The tests check the Butterworth magnitude against the bilinear-transformed closed form for several orders, and the Chebyshev passband ripple bounds.

References: S. Butterworth, “On the theory of filter amplifiers,” Experimental Wireless and the Wireless Engineer 7 (1930), 536-541.

Butterworth’s maximally flat filter

Butterworth's maximally flat filter

1947-1952 – The Z-Transform #

Radar fire-control systems of the Second World War sampled their inputs at discrete times, and the Laplace transform no longer fit. Witold Hurewicz, working on servomechanisms at the MIT Radiation Laboratory, described in 1947 how to treat sequences with a generating function in a complex variable. John Ragazzini and Lotfi Zadeh at Columbia named it the z-transform in 1952:

\[X(z) = \sum_{n=0}^{\infty} x[n]\,z^{-n}, \qquad \mathcal{Z}\{x[n-1]\} = z^{-1}X(z).\]

A delay becomes multiplication by \(z^{-1}\), so a linear difference equation becomes a rational function \(H(z) = B(z)/A(z)\). Each pole \(p\) contributes a mode \(p^n\) to the impulse response, so the system is stable exactly when every pole lies inside the unit circle. On the circle itself, \(z = e^{i\omega}\), the transform is the discrete-time Fourier transform. The idea is the discrete counterpart of Laplace’s, and ultimately of de Moivre’s generating functions.

Implementation: mathematicskit.special_functions.systems.z_transform.z_transform() and transfer_function() evaluate the defining sums by Horner’s rule; poles_zeros() wraps scipy.signal.tf2zpk() and reports stability; and inverse_z_transform() expands \(H(z)\) in partial fractions with scipy.signal.residuez() and inverts each term, including repeated poles. The tests check the geometric-series pair, the DFT on the unit circle, the double-pole pair \((n+1)p^n\), and the damped cosine from a complex pole pair.

References: W. Hurewicz, “Filters and servo systems with pulsed data,” in H. M. James, N. B. Nichols, and R. S. Phillips, eds., Theory of Servomechanisms, MIT Radiation Laboratory Series 25 (New York: McGraw-Hill, 1947), Ch. 5; J. R. Ragazzini and L. A. Zadeh, “The analysis of sampled-data systems,” Transactions of the AIEE 71, Part II (1952), 225-234; E. I. Jury, Theory and Application of the z-Transform Method (New York: Wiley, 1964).

The Z-transform: poles, stability, and partial fractions

1966 – Fast Convolution #

Convolution, \((x * h)[n] = \sum_k x[k]\,h[n-k]\), is how every linear time-invariant system acts on its input, and evaluated directly it costs \(O(NM)\) operations. The convolution theorem, \(\mathcal{F}\{x * h\} = \mathcal{F}\{x\}\,\mathcal{F}\{h\}\), had been known since the nineteenth century, but computing Fourier transforms was itself too slow for it to help. A year after Cooley and Tukey’s FFT, Thomas Stockham at MIT showed that zero-padding both sequences to length \(N + M - 1\), multiplying their FFTs, and transforming back gives the linear convolution in \(O((N+M)\log(N+M))\) time. Without the padding, the product gives the circular convolution instead. Stockham went on to use FFT-based deconvolution in 1975 to restore Enrico Caruso’s acoustic recordings, and FFT convolution is now how long filters, reverberation, and cross-correlation are computed.

Implementation: mathematicskit.special_functions.systems.convolution.convolve_direct() and convolve_fft() wrap numpy.convolve() and scipy.signal.fftconvolve(); circular_convolve() applies the DFT convolution theorem directly; cross_correlate() wraps scipy.signal.correlate(); and compare_convolution_methods() times the direct and FFT methods against each other. The tests check polynomial multiplication, agreement of the two methods, circular convolution as wrapped linear convolution, and recovery of a known shift from a correlation peak.

References: T. G. Stockham Jr., “High-speed convolution and correlation,” AFIPS Conference Proceedings 28, Spring Joint Computer Conference (1966), 229-233; A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed. (Upper Saddle River: Pearson, 2010), Sec. 8.7.

Fast convolution: the convolution theorem and the FFT

1970-1979 – Numerical Inversion of the Laplace Transform #

When \(F(s)\) is known only numerically, or has no table entry, Bromwich’s inversion integral must be computed, and along a vertical line it oscillates without decaying. Harald Stehfest’s 1970 algorithm, refining Donald Gaver’s 1966 probabilistic method, avoids complex arithmetic entirely:

\[f(t) \approx \frac{\ln 2}{t}\sum_{k=1}^{n} V_k\,F\!\left(\frac{k\ln 2}{t}\right),\]

with alternating integer weights \(V_k\) that grow like \(10^{n/2}\), so each extra term costs digits to cancellation and oscillating \(f\) cannot be resolved. Alan Talbot’s 1979 method instead deforms the Bromwich line into a contour that wraps around the negative real axis, where \(e^{st}\) decays rapidly, so the trapezoidal rule on it converges geometrically. Joseph Abate and Peter Valkó’s 2004 fixed Talbot choice of contour needs a single parameter, the number of nodes \(m\), and reaches close to machine precision.

Implementation: mathematicskit.special_functions.systems.laplace_transform.inverse_laplace_stehfest() (with stehfest_coefficients(), computed in exact rational arithmetic) and inverse_laplace_talbot() are hand-rolled, since neither numpy nor scipy offers a numerical inverse Laplace transform. The tests invert closed-form pairs, the oscillating \(\sin t\), and the non-rational \(e^{-\sqrt s}/s\), check the Stehfest weights’ identities, and show that Talbot’s contour must enclose oscillating poles.

References: D. P. Gaver Jr., “Observing stochastic processes, and approximate transform inversion,” Operations Research 14 (1966), 444-459; H. Stehfest, “Algorithm 368: Numerical inversion of Laplace transforms,” Communications of the ACM 13 (1970), 47-49; A. Talbot, “The accurate numerical inversion of Laplace transforms,” IMA Journal of Applied Mathematics 23 (1979), 97-120; J. Abate and P. P. Valkó, “Multi-precision Laplace transform inversion,” International Journal for Numerical Methods in Engineering 60 (2004), 979-993.

Numerical Laplace inversion: Stehfest and Talbot

1974 – Kaiser’s Window and FIR Filter Design #

The ideal lowpass filter has a sinc-shaped impulse response that never ends. Truncating it to \(N\) taps gives a finite impulse response (FIR) filter with exactly linear phase, but the abrupt truncation leaves Gibbs ripple in the stopband of about \(-21\) dB, however long the filter. Tapering the truncation with a smooth window reduces the ripple at the cost of a wider transition band. James Kaiser at Bell Labs looked for the best window of a given length and found a near-optimal family built from the modified Bessel function \(I_0\):

\[w[n] = \frac{I_0\!\left(\beta\sqrt{1 - \left(\frac{2n}{N-1} - 1\right)^2}\right)}{I_0(\beta)}, \qquad 0 \le n \le N - 1.\]

The single parameter \(\beta\) trades main-lobe width against stopband attenuation, and Kaiser gave empirical formulas for the \(\beta\) and \(N\) that meet a given specification.

Implementation: mathematicskit.special_functions.systems.filters.fir_window_filter() wraps scipy.signal.firwin() with any window, including ("kaiser", beta). The tests check that the filter is symmetric (linear phase), has unit DC gain, and has the phase of a pure delay of \((N-1)/2\) samples.

References: J. F. Kaiser, “Nonrecursive digital filter design using the \(I_0\)-sinh window function,” Proceedings of the IEEE International Symposium on Circuits and Systems (1974), 20-23.

Kaiser’s window and FIR filter design

Kaiser's window and FIR filter design

1982-1984 – Morlet’s Continuous Wavelet Transform #

Jean Morlet, a geophysicist at Elf Aquitaine, analysed seismic echoes whose high frequencies last only briefly while low frequencies persist. Dennis Gabor’s 1946 windowed Fourier transform uses one window width for every frequency and so cannot follow both. Morlet instead correlated the signal with shifted and dilated copies of one Gaussian-windowed oscillation, so the window shrinks as the frequency rises:

\[W(s, \tau) = \frac{1}{\sqrt{s}}\int x(t)\,\overline{\psi\!\left(\frac{t-\tau}{s}\right)}\,dt, \qquad \psi(\eta) = \pi^{-1/4}e^{i\omega_0\eta}e^{-\eta^2/2}.\]

His 1982 papers with colleagues introduced the method; his 1984 paper with the physicist Alex Grossmann proved that the signal can be reconstructed from \(W\), and gave the name wavelet its mathematical meaning. The squared modulus, the scalogram, shows which frequencies are present when.

Implementation: mathematicskit.special_functions.systems.wavelets.morlet_cwt() computes the transform scale by scale as a product of FFTs, following Torrence and Compo’s normalization, and returns a ScalogramResult with the Fourier frequency of each scale. It is hand-rolled because scipy 1.15 removed scipy.signal.cwt. The tests check that a pure tone peaks at its own frequency and that a frequency switch is located in time.

References: J. Morlet, G. Arens, E. Fourgeau, and D. Giard, “Wave propagation and sampling theory,” Geophysics 47 (1982), 203-236; A. Grossmann and J. Morlet, “Decomposition of Hardy functions into square integrable wavelets of constant shape,” SIAM Journal on Mathematical Analysis 15 (1984), 723-736; C. Torrence and G. P. Compo, “A practical guide to wavelet analysis,” Bulletin of the American Meteorological Society 79 (1998), 61-78.

Grossmann and Morlet’s continuous wavelet transform

Grossmann and Morlet's continuous wavelet transform

1988 – Daubechies Wavelets #

Haar’s wavelet is compactly supported but discontinuous, and the smooth orthonormal wavelets Yves Meyer built in 1985-1986 had infinite support. Ingrid Daubechies found wavelets that have both properties. The wavelet with \(p\) vanishing moments, orthogonal to every polynomial of degree below \(p\), comes from a lowpass filter \(h\) of only \(2p\) taps, obtained by factorizing a polynomial \(P\) that she wrote down explicitly:

\[|H(\omega)|^2 = 2\cos^{2p}(\omega/2)\,P\!\left(\sin^2\tfrac{\omega}{2}\right), \qquad P(y) = \sum_{k=0}^{p-1}\binom{p-1+k}{k}y^k .\]

For \(p = 2\) the filter is \((1+\sqrt3,\,3+\sqrt3,\,3-\sqrt3,\,1-\sqrt3)/(4\sqrt2)\). Stéphane Mallat’s 1989 multiresolution analysis showed that the whole transform is a cascade of these filters followed by downsampling, costing only \(O(N)\). Because smooth pieces of a signal produce near-zero detail coefficients, discarding small coefficients compresses it; the JPEG 2000 image standard is built on the same idea.

Implementation: mathematicskit.special_functions.systems.wavelets.daubechies_filter() computes the filter by spectral factorization of \(P\), and discrete_wavelet_transform() implements Mallat’s pyramid with periodic boundaries. Both are hand-rolled, since scipy 1.15 removed scipy.signal.daub. The tests check the closed-form \(p = 2\) filter, orthonormality and vanishing moments up to \(p = 8\), and vanishing detail coefficients on a linear signal.

References: I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Communications on Pure and Applied Mathematics 41 (1988), 909-996; S. G. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence 11 (1989), 674-693; I. Daubechies, Ten Lectures on Wavelets (Philadelphia: SIAM, 1992).

Daubechies wavelets: vanishing moments and compression