This is a description of how to use Spigot algorithms to compute mathematical constants. It is based on Rabinowitz and Wagon's paper "A Spigot Algorithm for the Digits of π", but also adds several innovations and generalizations that might perhaps not be in the literature.

This is probably the most efficient way to solve the holes for the mathematical constants described in this article for languages that don't have built-in bignum support, and in particular seems most efficient for Assembly.

If you do have built-in bignum support, then it's probably more efficient to compute $10^k C$ directly in bignums, perhaps using the series expressions detailed here (unless the language has built-in arbitrary-precision special functions), using add/multiply loops if you choose the "sum of product" versions.

Note that depending on how exactly you implement and optimize the loops in your language, you might need to replace $j$ with $j \pm 1$ or $j \pm 2$ in $n_j$ and $d_j$.

The idea of spigot algorithms is to express numbers as this sum, known as a "mixed-base representation", which is an extension of the customary base-N notation (the special case of $n_j = 1$ and $d_j = N$):

$$ A = a_0 + \frac{n_1}{d_1} (a_1 + \frac{n_2}{d_2} (a_2 + \frac{n_3}{d_3} \dots $$

where $|a_j| \lt d_j$, $d_j \gt 0$.

Or more formally:

$$ A = \lim\limits_{P\to\infty} \sum\limits_{k=0}^{P} a_k \prod\limits_{j=1}^k \frac{n_j}{d_j} $$

Note that if $A = \sum\limits_{k=0}^{\infty} A_k$, then the above holds with $a_j = A_0$, and $n_j$ and $d_j$ such that $\frac{n_j}{d_j} = \frac{A_j}{A_{j-1}}$ (by product telescoping).

The key advantage is that it is easy to print this number in base 10 (or any base) with this algorithm, which is just a minor adaptation of the algorithm to print arbitrary-precision fractional bignums, as long as this condition always holds for any possible intermediate value produced by the algorithm:

$$ 0 \le A - a_0 = \sum\limits_{k=1}^{P} a_k \prod\limits_{j=1}^k \frac{n_j}{d_j} \lt 1 $$

In case that $n_j \ge 0$ and $a_j \ge 0$, since we have that $|a_j| < d_j$, the condition holds if this inequality is true, which can be verified by symbolic computation:

$$ \sum\limits_{k=1}^{P} (d_j - 1) \prod\limits_{j=1}^k \frac{n_j}{d_j} \lt 1 $$

# initialize a0 and a_j for the constant to print, make sure that a_j < d_j

print(to_string(a0) + ".")

for _ in 0..OUTPUT_DIGITS {
  # set A = 10 (A - a0), and reduce to maintain |a_j| < d_j
  C := 0 # actually any value works if P is big enough
  for j in P..1 {
    t = a_j * 10 + C
    q = t // d_j
    r = t % d_j
    a_j := r
    C := q * n_j
    # note it's also possible to multiply by n_j at the start of the loop
    # and in a0 = C (but usually n_1 = 1, so the latter is a no-op)
  }
  a0 = C
  # here we still have |a_j| < d_j due to the modulus operation,
  # and a_j >= 0 is preserved if a_j >= 0 and n_j >= 0 initially

  print(a0)
}

This algorithm works because $a_0$ is indeed the integer part of A and thus the correct digit because the pre-condition states that $A - a_0 \lt 1$; the algorithm thus prints $a_0$ and then recursively prints $10(A - a_0)$, which contains the rest of the digits.

$$ e = \sum\limits_k^{\infty} \frac{1}{k!} = \lim\limits_{P\to\infty} \sum\limits_{k=0}^{P} 1 \prod\limits_{j=1}^k \frac{1}{j} $$

This can be expressed by $a_j = 1$, $n_j = 1$, $d_j = j$.

We have that $A - a_0 = \sum\limits_{k=1}^{P} a_k \prod\limits_{j=1}^k \frac{n_j}{d_j} \lt \sum\limits_{k=1}^{P} (d_k - 1) \prod\limits_{j=1}^k \frac{n_j}{d_j} = \sum \limits_{k=1}^{P} (k - 1) \prod\limits_{j=1}^k \frac{1}{j} = \sum \limits_{k=1}^{P} \frac{k - 1}{k!} \lt \sum \limits_{k=1}^{\infty} \frac{k - 1}{k!} = 1$ so the algorithm works.

For most other constants such as $\pi$, the algorithm cannot be applied directly because condition $\sum\limits_{k=1}^{P} (d_k - 1) \prod\limits_{j=1}^k \frac{n_j}{d_j} \lt 1$ does not hold with the simplest way to express the constant.

The Rabinowitz paper presents "predigits" as a mechanism to handle this, and while it works it's not efficient for code golfing purposes (and probably not for any other purpose either) since it requires a lot of extra code.

Instead, we modify the constant by adding a "buffer" term:

$A = a_0 + \frac{1}{B} (a_1 + \frac{n'_1}{d'_1} (a_2 + \frac{n'_2}{d'_2} (a_3 + \frac{n'_3}{d'_3} \dots $

The idea is that instead of printing digits right away, we hold them in the buffer term $a_1$, and only output them once the value is larger than $B$, which for large enough $B$ will result in the digits not "changing" and thus not requiring a pre-digit mechanism.

More formally, we express the constant $\frac{C}{10^r}$ as $a_0 = 0$, $a_j = \frac{B}{10^r} a^\prime_{j-1}$, $n_1 = 1$, $d_1 = B$, $n_j = n^\prime_{j-1}$, $d_j = d'_{j-1}$, where $a'_j$, $n'_j$, $d'_j$ is the basic expression for the constant $C$.

We then have that as required:

Now the condition that needs to hold becomes:

This is equivalent to $S' \lt B - a_1$. As long as $S'$ converges to value that fits in an integer data type, since $S'$ is constant and the $a_1$ values will be pseudo-random, this is highly likely to hold for large enough $B$, and thus it's possible to apply the algorithm to the modified sequence (there may be a way to mathematically prove that there are infinite values of $B$ that satisfy the inequality at least for specific constants, but for code golfing finding it for specific cases is good enough).

Depending on the language used, it may be most efficient to set $B = 10^r$ for some $r$, so that the $a_j$ are small, or set a specific value of $r$, so that the algorithm output is most efficiently adjusted to the proper output.

$$ \pi = 2 \sum_{k=0}^\infty\frac{k!}{(2k+1)!!} = \lim\limits_{P\to\infty} \sum\limits_{k=0}^{P} 2 \prod\limits_{j=1}^k \frac{j}{2j + 1} $$

Hence $a_j = 2$, $n_j = j$, $d_j = 2j + 1$ expresses $\pi$. For $\tau$, set $a_j = 4$ instead. You'll need to apply buffering as described above.

$$ \frac{\ln 2}{10} = \frac{1}{10} \sum_{n=1}^\infty \frac{1}{2^{n}n} = \frac{1}{10} \lim\limits_{P\to\infty} \sum\limits_{k=0}^{P} 5 \prod\limits_{j=1}^k \frac{j}{2j + 2} $$

Hence $a_j = 5$, $n_j = j$, $d_j = 2j + 2$ expresses $\frac{\ln 2}{10}$. You'll need to apply buffering as described above.

From symbolic computation:

$$ \sqrt{C} = \lim\limits_{P\to\infty} \sum\limits_{k=0}^{P} \prod\limits_{j=1}^k \frac{(C - 1)(2j - 1)}{2Cj} $$

Hence $a_j = 1$, $n_j = 2j - 1$, $d_j = 4j$ expresses $\sqrt{2}$. You'll need to apply buffering as described above.

Setting $C=\frac{5}{4}$, $a_j = 1$, $n_j = 2j - 1$, $d_j = 10j$ expresses $\frac{\sqrt{5}}{2}$. You'll need to apply buffering as described above.

$\phi = \frac{\sqrt{5}}{2} + \frac{1}{2}$ which is easily obtained by changing a decimal digit from $\frac{\sqrt{5}}{2}$.

Buffering can also supports cases where the $n_j$ or $a_j$ can be negative.

In this case, we need to calculate $S'$ using $a_j = \pm (d_j - 1)$ choosing the sign so that all terms are positive, and then $\frac{1}{B} (a_1 - S') \ge 0$ also needs to hold (except for the first iterations, where we can verify that it holds manually), which is equivalent to $a_1 \ge S'$. With the same argument used before, it's usually possible to find $B$ such that all $a_1$ values satisfy this relation.

Hence $a_j = 1$, $n_j = -j^3(3j+2)$, $d_j = (2j+1)^3(3j-1)$ expresses $G$. You'll need to apply buffering as described above.

Hence $a_j = 125$, $n_j = -j^3$, $d_j = 2(2j+1)(j+1)^2$ expresses $\zeta(3)$. You'll need to apply buffering as described above.

The Euler-Mascheroni constant $\gamma$ can also be computed by a Spigot algorithm, but it's more complicated.

This identity holds for exponential integrals:

$$ E_1(z) = -\gamma - \ln z - \sum_{k=1}^{\infty} \frac{(-z)^k}{k \cdot k!} $$

Rearranging and setting $z = 2^n$,

$$ \gamma = - n \ln 2 - \sum\limits_{k=1}^{\infty} \frac{(-2^n)^k}{k \cdot k!} - E_1(2^n) $$

We have $\lim_{x\to\infty} E_1(x) = 0$, so

$$ \gamma = \lim_{n\to\infty} -n \ln 2 -\sum\limits_{k=1}^{\infty} \frac{(-2^n)^k}{k \cdot k!} $$

The digits of $n \ln 2$ can be computed like $\ln 2$, but multiplying $a_j$ by $n$, and the other expression can only be computed as described below. It is then possible to combine the two spigot algorithms by adding or subtracting the "buffered" term from the first one to the second one and using the second one to output digits, and it's usually possible to find a $B$ such that this works properly.

For the second part, we have

However, we can't directly use $n_j = -2^n j$ and $d_j = (j+1)^2$ (unless we have bignum support in the language), because while the sum $S' = \sum\limits_{k=1}^{\infty} \frac{k^2-1}{2^n} \frac{(2^n)^k}{k \cdot k!}$ converges, it converges to a huge value (approximately $e^{2^n}$) that doesn't allow to find a value $B$ for buffering that fits in a 64-bit integer.

So, instead, we use the simpler mixed-base representation $n_j = j$ and $d_j = (j+1)^2$ (which gives $S' = 1.4004$), and compute the sum using arithmetic in the mixed-base representation.

For this purpose, we manipulate the expression like this:

This can be computed by starting with $a_j = 0$ and then executing a loop with $v$ from $P$ to $1$, where in each iteration we add 1 to $a_v$, and then multiply the value $A$ by $-2^n$ (using the same multiplication algorithm that is used to multiply by 10 in the printing algorithm).

All this together allows to compute the Euler-Mascheroni constant with spigot algorithms. For 1000 digits, $n = 12$ should be sufficient.