Gauss’s constant

I hadn’t heard of Gauss’s constant until recently. I imagine I’d seen it before and paid no attention. But once I paid attention, I started seeing it more often. There’s a psychology term for this—reticular activation?—like when you buy a green Toyota and suddenly you see green Toyotas everywhere.

Our starting point is the arithmetic-geometric mean or AGM. It takes two numbers, takes the arithmetic (ordinary) mean and the geometric mean, then repeats the process over and over. The limit of this process is the AGM of the two numbers.

Gauss’s constant is the reciprocal of the AGM of 1 and √2.

g \equiv \frac{1}{\text{AGM}(1, \sqrt{2})} = 0.834686\ldots

Gauss’s constant can be expressed in terms of the Gamma function:

g = (2\pi)^{-3/2} \, \Gamma\left(\frac{1}{4}\right)^2

Exponential sums

Last week I wrote about the sum

\sum_{n=-\infty}^\infty e^{-\pi n^2} = \frac{\sqrt[4] \pi}{\Gamma\left(\frac 3 4\right)}

which we will see shortly is related to Gauss’s constant.

We can the definition of g in terms of Γ(1/4) and the reflection identify for the gamma function to show that the sum above can be written in terms of Gauss’s constant:

\sum_{n=-\infty}^\infty e^{-\pi n^2} = 2^{1/4} \sqrt{g}

The alternating version of the sum is also related to g:

\sum_{n=-\infty}^\infty (-1)^n e^{-\pi n^2} = g^2

Integrals

Another place where g comes up is in integrals of hyperbolic functions. For example

\begin{align*} \int_0^\infty \sqrt{\text{sech}(x)} \, dx &= \pi g \\ \int_0^\infty \sqrt{\text{csch}(x)} \, dx &= \sqrt{2} \pi g \end{align*}

Here

\begin{align*} \text{sech}(x) &= \frac{1}{\cosh(x)} = \frac{2}{e^x + e^{-x}} \\ \text{csch}(x) &= \frac{1}{\sinh(x)} = \frac{2}{e^x - e^{-x}} \\ \end{align*}

Beta probabilities

Another place Gauss’s constant comes up is in special values of beta distribution probabilities.

Define the incomplete beta function by

B(a, b, x) = \int_0^x t^{a-1} (1-t)^{b-1}\, dt

If X is a random variable with a beta(a, b) distribution, then the CDF of X is the normalized incomplete beta function, i.e.

\text{Prob}(X \leq x) = \frac{B(a, b, x)}{B(a, b, 1)}

A couple special values of B involve Gauss’s constant, namely

B\left(\frac{1}{4}, \frac{1}{4}, \frac{1}{2}\right) = \sqrt{2} \pi g

and

B\left(\frac{3}{4}, \frac{3}{4}, \frac{1}{2}\right) = \frac{1}{\sqrt{2} g}

Source:  An Atlas of Functions

Ubiquitous constant

The constant M defined by 1/M = √2 g is the so-called “ubiquitous constant,” though this seems like an unjustified name for such an arcane constant. Perhaps there is some context in which the constant is ubiquitous in the colloquial sense.

Gamma of integer plus one over integer

The gamma function satisfies

Γ(x+1) = x Γ(x)

and so in principle you could calculate the gamma function for any positive real number if you can calculate it on the interval (0, 1). For example,

\begin{align*} \Gamma(\pi) &= (\pi - 1)\Gamma(\pi -1) \\ &= (\pi - 1)(\pi - 2) \Gamma(\pi - 2) \\ &= (\pi - 1)(\pi - 2) (\pi - 3)\Gamma(\pi - 3) \end{align*}

So if you’re able to compute Γ(π-3) then you could compute Γ(π).

If n is a positive integer ε is some number between 0 and 1, is there a more direct way to express Γ(n + ε) in terms of Γ(ε)?

There is when ε is the reciprocal of an integer. For example,

\begin{align*} \Gamma\left(n + \frac{1}{2}\right) &= \frac{(2n-1)!!}{2^n} \Gamma\left(\frac{1}{2}\right) \\ \Gamma\left(n + \frac{1}{3}\right) &= \frac{(3n-2)!!!}{3^n} \Gamma\left(\frac{1}{3}\right) \\ \Gamma\left(n + \frac{1}{4}\right) &= \frac{(4n-3)!!!!}{4^n} \Gamma\left(\frac{1}{4}\right) \\ \end{align*}

The multiple exclamation marks may look strange if you haven’t seen this notation before. See the previous post for an explanation.

The general equation for integer k is

\Gamma\left(n + \frac{1}{k}\right) = \frac{\Pi_k(kn - k + 1)}{k^n} \Gamma\left(\frac{1}{k}\right)

where Πk is a notation for k-factorial I suggested in the previous post. A more common notation would be to put !(k) to the right of the argument rather than Πk on the left.

Here’s Python code demonstrating that the equation above holds for randomly selected n and k.

from operator import mul
from functools import reduce
from scipy.special import gamma
import numpy as np

def multifactorial(n, k):
    return reduce(mul, range(n, 0, -k), 1)

def f1(n, k):
    return gamma(n + 1/k)

def f2(n, k):
    return multifactorial(k*n - k + 1, k)*gamma(1/k)/k**n

np.random.seed(20211014)
for _ in range(1000):
    n = np.random.randint(1, 50)
    k = np.random.randint(1, 50)    
    a, b = f1(n,k), f2(n, k)
    assert( abs(a/b - 1) <= 1e-10 )

Multifactorial

The factorial of an integer n is the product of the positive integers up to n.

The double factorial of an even (odd) integer n is the product of the positive even (odd) integers up to n. For example,

8!! = 8 × 6 × 4 × 2

and

9!! = 9 × 7 × 5 × 3 × 1.

Double factorials come up fairly often, and sometimes triple or higher multifactorials do too.

In general, the k-factorial of n is the product of the positive integers up to n that are congruent to n mod k. Here’s how you might implement k-factorials in Python.

    from operator import mul
    from functools import reduce

    def multifactorial(n, k):
        return reduce(mul, [i for i in range(1, n+1) if (n-i)%k == 0], 1)

Update: As pointed out in the comments, multifactorial could be implemented more sucinctly as

    def multifactorial(n, k):
        return reduce(mul, range(n, 0, -k), 1)

Factorial notation is a little awkward, and multifactorial notation is even more awkward. Someone seeing n!! for the first time might reasonably assume it means (n!)!, but it does not.

Adding exclamation points works, albeit awkwardly, for specific values of k, but not for variable numbers of k. One way I’ve seen to express k-factorials for variable k is to add a subscript (k) to the factorial:

n!_{(k)} \equiv n\underbrace{! \cdots !}_\text{$k$ times}

I’ve also seen the subscript on the exclamation point written as a superscript instead.

I’d like to suggest a different notation: Πk. by analogy with the Π function.

\Pi_k(n) = \prod_{\stackrel{1 \,\leq \,i \leq \,n}{i \,\equiv\, n \pmod{k}}} i

When k ≡ 1 the condition in (mod k) always holds, and we have factorial, i.e.

Π1(n) = Π(n) = n!.

One downside of this notation is that while the function Π is defined for all complex numbers (aside from singularities at negative integers) the function Πk is only defined for positive integer arguments. Still, in my opinion, this is less awkward than the alternatives I’ve seen.

By the way, the term “double factorial” seems backward. Maybe it should have been “half factorial” because you’re multiplying half as many numbers together, not twice as many. “Multifactorial” in general seems like an unfortunate name. Subfactorial might have been better, but unfortunately that name means something else.

I understand why someone may have come up with the idea of using two exclamation points for what we call double factorial; it would be hard to draw half an exclamation point. An advantage of the notation suggested here is that it suggests that there’s a variation on factorial somehow associated with k, but there’s no implication of multiplying or dividing by k.

Adding “mod k” as a subscript would be even more explicit than a subscript k. Someone who hadn’t seen

Πmod k (n)

before might be able to figure out what it means; I think they would at least easily remember what it means once told. And the presence of “mod” might suggest to the reader that the argument needs to be an integer.

Discrete sum analog of Gaussian integral

A comment on a recent post lead me to a page of series on Wikipedia. The last series on that page caught my eye:

\sum_{n=-\infty}^\infty e^{-\pi n^2} = \frac{\sqrt[4] \pi}{\Gamma\left(\frac 3 4\right)}

It’s a lot more common to see exp(-πx²) inside an integral than inside a sum. If the summation symbol were replaced with an integration sign, the integral would be 1. You could derive this from scratch using the polar coordinate trick, of your could look at it and say that the integrand is the PDF of a Gaussian random variable with mean 0 and variance 1/2π.

Last week I wrote about the Pi function Π(z) = Γ(z + 1) and how some equations are simpler or easier to remember when expressed in terms of Π rather than Γ, and the equation above is an example. It can be rewritten as

\sum_{n=-\infty}^\infty e^{-\pi n^2} = \frac{1}{\pi^{-\frac{1}{4}}\,\Pi(-\frac{1}{4})}

It’s curious that there’s such a similarity between the constant pi and the function Pi in the denominator. After all, what’s the connection between the two? Their names come from “perimeter” and “product”, both which start with “p.” That explains a linguistic connection but not a mathematical connection. And yet π and Π often appear together. For example, the volume of an n-dimensional sphere is

\frac{\pi^{\frac{n}{2}}}{\Pi(\frac{n}{2})}

We could rewrite this changing πn/2 in the numerator to πn/2 in the denominator, making it look more like the equation for the sum above, except the sign of the exponent on π would be the opposite of the sign on the argument to Π.

Incidentally, the summation result above says that we could define a probability distribution on the integers with probability mass function

p(n) = \frac{1}{\pi^{1/4} \Pi(1/4)} \exp(-\pi n^2)

which would be a sort of discrete analog of the normal distribution. I imagine someone has done this before and given it a name.

The Calculus of Finite Differences

The Calculus of Finite Differences by L. M. Milne-Thompson is a classic. It covers a great deal of elegant and useful [1] mathematics that isn’t widely known, at least not any more. For a taste of the subject matter of the book, see this post.

The book is now in the public domain—at least in the United States—and so you can find scans of it online. The copy I downloaded and a hard copy that I purchased both have some problems. I’m writing up how to fix the problems in case anyone else runs into the same difficulty.

I bought a hard copy of the book because I prefer reading from paper, though I also like having PDFs for reference. I thought I would patch the PDF scan by scanning the pages of my hard copy corresponding to the pages missing in the PDF. However, the scan and the paperback version have the exact same errors.

Pages 1, 2, 31, and 32 are missing. Actually, pages 1 and 2 are included, but not where you’d expect. After page 396 comes page 2 and then page 1. Pages 31 and 32 are completely missing.

I went back to the Internet Archive and found a different scan of the same book, one that appears to be complete.

My guess is that the errors went back to the first printing in 1933, and were corrected in subsequent editions. Other than correcting the missing pages, I don’t know of any changes between editions.

***

[1] The author describes his motivation for writing the book as follows.

Two distinct reasons impelled me to undertake this work. First, in my lectures at Greenwich to junior members of the Royal Corps of Naval Constructors I have occasion to treat certain aspects of difference equations; secondly, the calculation of tables of elliptic functions and integrals, on which I have been recently engages, give rise to several interesting practical difficulties which had to be overcome.

Multiple angle identity for cotangent

The previous post discusses generalized replicative functions, functions that satisfy

a_n f(nx) + b_n = \sum_{i=0}^{n-1} f\left(x + \frac{i}{n} \right )

Donald Knuth says in Exercise 40 of Section 1.2.4 of TAOCP that the functions cot πx and csc² πx fall into this class of functions.

If this is true of cot πx then it follows by differentiation that it is true of csc² πx as well, so we will focus on the former.

I believe that for f(x) = cot(πx) we have an = n and bn = 0, but I don’t have a proof. If this is correct then

n \cot (n \pi x) = \sum_{i=0}^{n-1} \cot\left(\pi\left(x + \frac{i}{n} \right ) \right )

Here’s a Python script that draws values of n and x at random and verifies that the left and right sides above are equal to at least 10 decimal places. It’s not a proof, but it’s good evidence.

    import numpy as np

    def cot(x): return 1/np.tan(x)

    def f1(n, x): return n*cot(n*np.pi*x)

    def f2(n, x): return sum(cot(np.pi*(x + k/n)) for k in range(n))

    np.random.seed(20211011)
    for _ in range(1000):
        n = np.random.randint(1, 50)
        x = np.random.random()
        a, b = f1(n,x), f2(n, x)
        assert( abs(a - b) <= abs(a)*1e-10 )

Note that NumPy doesn’t have a cotangent function. See this post for quirks of how trig functions are supported in various environments.

An actual proof is left as an exercise for the reader. Knuth left it as an exercise too, but he didn’t say what the a‘s are, so I’m being nicer than he was. :)

If you come up with a proof, please share it in the comments.

Update

Here’s a proof combining the first two comments below. Start with the identity

\sin(nx) = 2^{n-1}\prod_{k=1}^{n-1}\sin\left(x+\frac{k\pi}{n}\right)

and substitute πx for x. Take the logarithm and then the derivative of both sides, divide by π, and you get the identity we’re trying to prove.

Multiplication theorem rabbit hole

When I started blogging I was reluctant to allow comments. It seems this would open the door to a flood of spam. Indeed it does, but nearly all of it can be filtered out automatically.

The comments that aren’t spam have generally been high quality. A comment on my post about the sawtooth function a few days ago has sent me down a fascinating rabbit hole.The sawtooth function belongs to a general class of functions called replicative functions, functions that satisfy the equation

f(nx) = \sum_{i=0}^{n-1} f\left(x + \frac{i}{n} \right )

for positive integers n.

I ran across this by reading Knuth’s TAOCP. Knuth gives an understated comment in section 1.2.4 exercise 40 that “it may be interesting to study also the more general class of functions” for which the left-hand side above includes constants that depend on n, i.e.

a_n f(nx) + b_n = \sum_{i=0}^{n-1} f\left(x + \frac{i}{n} \right )

where the constants an and bn depend on n but not on x.

Bernoulli polynomials

The first Bernoulli polynomial B1(x) = x – ½ is a replicative function, and the only continuous replicative functions are multiples of B1 (proved in exercise 40 referenced above).

The rest of the Bernoulli polynomials Bm(x) are generalized replicative functions, with an = nm-1 and bn = 0. This was discovered by Joseph Raabe in the nineteenth century.

Gamma and digamma functions

Gauss’s multiplication theorem for the gamma function says that for positive n,

This is equation 6.1.20 in A&S. If we take logs of both sides, we find that log Γ is a almost a generalized replicative function.

\log \Gamma(nz) = \frac{1-n}{2}\log(2\pi) + \left(nz - \frac{1}{2}\right)\log n + \sum_{k=0}^{n-1} \log \Gamma\left(z + \frac{k}{n} \right )

The factor of n log(n) z prevents the definition from holding, but if we take derivatives, this term goes away.

This means that the derivative of log Γ, a.k.a the ψ function, a.k.a. the digamma function, is a generalized replicative function with an = n and bn = –n log n.

The derivative of a generalized replicative function is another function in the same class, so the higher derivatives of log Γ, known as the polygamma functions, are also generalized replicative functions.

Hurwitz zeta function

The Hurwitz zeta function ζ(s, z) is defined by

\zeta(s, z) = \sum_{k=0}^\infty \frac{1}{(z+k)^s}

The Riemann zeta function is the special case of the Hurwitz zeta function with z = 1.

The Hurwitz zeta function is related to the polygamma functions by

\psi^{(m)}(z)= (-1)^{m+1} m! \zeta (m+1,z)

It follows from this that for fixed s, ζ(s, z) is a generalized replicative function of z. You can show that an = ns and bn = 0.

Trig functions

See the next post for a couple examples involving trig functions: cot πx and csc² πx .

Related posts

Gamma and the Pi function

The gamma function satisfies

Γ(n+1) = n!

for all non-negative integers n, and extends to an analytic function in the complex plane with the exception of singularities at the non-positive integers [1]. Incidentally, going back to the previous post, this is an example of a theorem that would have to be amended if 0! were not defined to be 1.

Wouldn’t it be simpler if the gamma function were defined so that it’s value at n, not at n+1, extended factorial? Well, some things would be simpler, but other things would not.

The Pi function defined by

Π(z) = Γ(z + 1)

extends factorial with no extra factor of 1 to keep up with, and some theorems are simpler to state in terms of the Pi function than the gamma function. For example, it’s simpler to say that the volume of a unit n-sphere is

(π/2)n/2 / Π(n/2)

than to say it’s

(π/2)n/2 / Γ(n/2 + 1).

The former has an easy-to-remember form, with lower case π on top and capital π on bottom.

The reflection identity is also a little simpler in the form

Π(z) Π(-z) = 1/sinc(z)

than in the form

Γ(z) Γ(1-z) = π / sin(πz)

The drawback to the former is that you have to remember to use the convention

sinc(z) = sin(πz) / πz

because some sources define sinc with the factor of π and some without. Neither convention makes a large majority of theorems simpler, so there’s no clear winner. [2]

Fortunately the Pi function has a different name and isn’t an alternative convention for defining the gamma function. That would be terribly confusing. [3]

Related posts

[1] The Pi function is the unique way to extend factorial to an analytic function, given some regularity assumptions. See Wielandt’s theorem and the Bohr-Mollerup theorems.

[2] Fourier transforms may be the worst example in mathematics of no clear convention winner. Multiple conventions thrive because each makes some some things easier to state. (See this page for a sort of Rosetta stone translating between the various Fourier transform conventions.) In fact, the lack of standardization for the sinc function is related to the lack of standardization for the Fourier transform: you’d like to define the sinc function so that it has the simplest Fourier transform under your convention.

[3] I’m fuzzy on the history of all this, but I think what we now call the Pi function was an alternative definition of the gamma function briefly, but fortunately that got ironed out quickly.

Sawtooth and replicative functions

Here’s something I ran across in Volume 2 of Donald Knuth’s The Art of Computer Programming.

Knuth defines the sawtooth function by

((x)) = x – (⌊x⌋ + ⌈x⌉)/2.

Here’s a plot.

This is an interesting way to write the sawtooth function, one that could make it easier to prove things about, such as the following theorem. For positive integers n, we have

((nx)) = ((x)) + \left(\left(x + \frac{1}{n}\right)\right) +\cdots+ \left(\left(x + \frac{n-1}{n}\right)\right)

It’s not at all obvious, at least to me, that this should be true.

Knuth does not give a proof but leaves the proof to Exercise 2 in Section 3.3.3. I suppose you could prove it by slugging it out with floor and ceiling identities, but I looked at the solutions to see whether Knuth has a more elegant proof.

The solution for Exercise 2 says to see exercises 38 and 39 in Volume 1, section 1.2.4. These exercises break the problem of proving the theorem above into smaller parts. They also generalize the problem, defining a function f to be replicative if it satisfies the equation

f(nx) = \sum_{i=0}^{n-1} f\left(x + \frac{i}{n} \right )

The theorem above says that the sawtooth function is replicative. Other examples of replicative functions include the floor function and the function that takes x to x – ½.

Related posts

Estimating normal tail extreme probabilities

In the previous post I said the probability of a normal distribution being 50 standard deviations from its mean was absurdly small, but I didn’t quantify it.

You can’t simply fire up something like R and ask for the probability. The actual value is smaller than can be represented by a floating point number, so you’ll just get zero.

Though there’s no practical reason to compute the probability of such events, we’ll do it anyway. The technique shown here is useful for less extreme cases that could be practical.

Let Z be a normal random variable with mean 0 and standard deviation 1. We need to estimatethe complementary CDF [1], given by

 \Phi^c(t) = P(Z > t) = \frac{1}{\sqrt{2\pi}} \int_t^\infty e^{-x^2/2}\, dx.

for t = 50. This post gives us the upper and lower bounds we need:

\frac{t}{t^2 +1} < \sqrt{2\pi} e^{t^2/2} \Phi^c(t) < \frac{1}{t}.

The same post also gives tighter bounds, which are useful for less extreme cases, but not necessary here.

The bounds above tell us the probability of Z taking on a value of 50 or more is approximately

\frac{e^{-1250}}{50\sqrt{2\pi}}

We’d like to write this value in scientific notation. We can’t directly evaluate it with most software because, as we said above, it’s too small. If you were to type exp(-2500) in Python, for example, it would return exactly 0. We can, however, evaluate the log base 10 of the expression above using Python.

    >>> from math import exp, pi, log10
    >>> -1250*log10(exp(1)) - 0.5*log10(2*pi) - log10(50)
    -544.9661623175799

So our probability is about 10-545. This is far smaller than the probability of selecting a particular elementary particle at random from the universe. (See here.)

Related posts

[1] Numerical software packages will provide functions for CDFs and CCDFS, e.g. for Φ and Φc = 1 – Φ. This may seem unnecessary, and it would be if we could compute with infinite precision. But it is possible, and often necessary, to compute Φc(x) for values of x so large that all precision would be lost when evaluating 1 – Φ(x). For less extreme values of x we wouldn’t lose all precision, but we’d lose some. See the discussion of erf and erfc and other functions that seem redundant here.