On A Spigot Algorithm for the Digits of π (Rabinowitz and Wagon)

Two equivalent expressions of a standard decimal number:

A mixed-radix number is somewhat like having a different base for each digit. The primary objective of this paper is stating mixed-radix representation’s of $e$ and $\pi$, and converting them to more useful bases.

Example of mixed-radix number:

Which expands to

When $a_i=1$ for all $i$, this is the well-known power series for $e$.

Euler’s number spigot

A. H. J. Sale doi 10.1093/comjnl/11.2.229.

Well-known $e$ power series:

Given in mixed-radix form:

Also can be stated as a recursion:

The paper then multiplies the fractional part of the expression by $\frac{10}{10}$, and distributes the numerator into the recursion:

The goal of this manipulation is to extract base-10 digits from this mixed-radix form, one digit at a time. We first take this recursion to some specific depth. $5$ in this example.

From right to left, we consecutively reduce the inner values mod the outer divisor, moving out the quotient in the process. An example is more clear:

This method of distributing the divisor, diving with remainder, and factoring is repeated:

This shows clearly that the 2nd digit is $7$. To get further digits, we again multiply the mixed-radix part by $10/10$.

And repeat the process:

And again:

And again:

Which now produces an incorrect result.

Why? Well, it is reasonable to assume that, because we limited the recursion to a finite number of levels, the remaining levels begin to become significant as the powers of $10$ rise.

Looking at the last reasonable result, we see at the end:

Let’s call $(\ldots)$ to be $r$ , we see

Which expands to

Multiplying in the $10^3$ creates a value less than $1$, and so it does not effect the calculations.

However, when the power of ten is increased to $10^4$, we see:

Which is greater than one and will likely cause a deleterious effect when it is missing from the calculations. So, the question is: how many terms does the recursion need to be expanded to to obtain $d$ digits? Instead of thinking about the algorithm as algebraic manipulations on a recursive equation, I find it easier to think of as an array of numbers. The first step of the above example becomes:

Where:

Inside the array, there are 5 slots (columns). Each slot represents a level of recursion. More slots => more digits extractable.

Let $d$ represent the number of digits being extracted, and $s$ represent the number of slots. A calculation using $s$ slots will yield $d$ digits when the following is true:

Infinite sums are cumbersome so instead a simpler form is desired. So long as the value of the sum overestimated, a simpler form will still work.

The following form is valid according to Desmos and WolframAlpha, so I’ll trust those sites on that.

Substituting this simpler expression in:

Solve for digits:

Plotting this, we see

The linear approximation of

Fits well under the entire domain. This approximation becomes poorer as more slots are added, but never overestimates. When $s=70$, $d \le \approx100$

A polynomial or piecewise approximation could provide improvement. The original paper provides a different approximation entirely. I have not reviewed that approximation.

When in one of the previous examples the result of $2.7180\ldots$ was obtained, only 5 slots were used, and so only two digits after the decimal place should have been expected with confidence to be correct.

Back again to thinking of the problem algorithmically:

Where:

To find the second digit, update the table such that:

• $v_i=10r_i$ for all $i$
• $q_i$ is cleared for all $i$

Then apply the same rules, stopping after $i=1$:

And so digit #2 is $1$: $2.71\ldots$.

An implementation in Rust is provided.

fn main() {
    for i in 1..20 {
        println!("{:?}", e(i));
    }
}

fn e(d: u32) -> Vec<u8> {
    let mut digits = vec![0; d as usize];
    let slots = d as usize + 3;
    let mut tmp = vec![10; slots];
    for stop in 0..d as usize {
        let mut carry = 0;
        for i in (0..slots).rev() {
            let n = i as u32 + 2;
            let s = carry + tmp[i];
            let q = s / n;
            let r = s % n;
            tmp[i] = r * 10;
            carry = q;
        }
        digits[stop] = carry as u8;
    }
    digits
}

Back to π

After a slight detour, the study of π is resumed.

Todo….