Pi digits (again) - an unbounded spigot program

[Guenter Schink]

I must admit that I didn't think about any optimization. Lazy me, simply copied Thomas' code and added some lines for the timing. Just to quickly show something that works in Python on Prime.

Günter

[Thomas Klemm]

(05-19-2024 12:34 PM)Guenter Schink Wrote:  I must admit that I didn't think about any optimization.

If speed were our primary concern, we should probably use a different algorithm:

Quote:Like Rabinowitz and Wagon’s algorithm, our proposals are not competitive with state-of-the-art arithmetic-geometric mean algorithms for computing \(\pi\).

What I like about the spigot algorithm is that one digit after the other appears instead of showing the entire result at the end.

[Thomas Klemm]

(05-18-2024 08:03 AM)EdS2 Wrote:  I never had a clue how to rewrite Haskell into something I could understand.

These are translations to Python that are closer to the original Haskell programs.
They use typing and itertools which is missing in MicroPython.
So I had to remove them in the program for the Prime.

I'm using version 3.12.2.
There were some changes in the area of typing, so it might not work with older versions of Python.

Rabinowitz and Wagon

Code:

from fractions import Fraction
from itertools import count, islice
from sys import argv
from typing import Callable, Generator

State = tuple[int, int, int, int]


def main(n: int, k: int):

    def stream(
        next: Callable[[State], int],
        safe: Callable[[State, int], bool],
        prod: Callable[[State, int], State],
        cons: Callable[[State, State], State],
        z: State,
        lfts: Generator[State, None, None],
    ) -> Generator[int, None, None]:
        while True:
            y = next(z)
            if safe(z, y):
                yield y
                z = prod(z, y)
            else:
                [x] = list(islice(lfts, 1))
                z = cons(z, x)

    def next(z: State) -> int:
        return int(extr(z, 3))

    def safe(z: State, n: int) -> bool:
        return n == int(extr(z, 4))

    def extr(lft: State, x: int) -> Fraction:
        """
        | q r | | x | = | q*x+r |
        | s t | | 1 |   | s*x+t |
        """
        q, r, s, t = lft
        return Fraction(q * x + r, s * x + t)

    def prod(z: State, n: int) -> State:
        """
        | 10 -10n | * z -> z
        |  0   1  |
        """
        return comp((10, -10 * n, 0, 1), z)

    def comp(a: State, b: State) -> State:
        """
        | q r | | u v | = | q*u+r*w q*v+r*x |
        | s t | | w x |   | s*u+t*w s*v+t*x |
        """
        q, r, s, t = a
        u, v, w, x = b
        return q * u + r * w, q * v + r * x, s * u + t * w, s * v + t * x

    cons = comp
    init = (1, 0, 0, 1)
    lfts = ((k, 4 * k + 2, 0, 2 * k + 1) for k in count(1))
    pi = stream(next, safe, prod, cons, init, lfts)
    for _ in range(n):
        print("".join(map(str, islice(pi, k))))


if __name__ == "__main__":
    n, k, *_ = argv[1:]
    main(int(n), int(k))

Lambert

Code:

from fractions import Fraction
from itertools import count, islice
from sys import argv
from typing import Callable, Generator

LFT = tuple[int, int, int, int]
State = tuple[LFT, int]


def main(n: int, k: int):

    def stream(
        next: Callable[[State], int],
        safe: Callable[[State, int], bool],
        prod: Callable[[State, int], State],
        cons: Callable[[State, LFT], State],
        z: State,
        lfts: Generator[LFT, None, None],
    ) -> Generator[int, None, None]:
        while True:
            y = next(z)
            if safe(z, y):
                yield y
                z = prod(z, y)
            else:
                [x] = list(islice(lfts, 1))
                z = cons(z, x)

    def next(z: State) -> int:
        [[q, r, s, t], i] = z
        x = 2 * i - 1
        return int(Fraction(q * x + r, s * x + t))

    def safe(z: State, n: int) -> bool:
        [[q, r, s, t], i] = z
        x = 5 * i - 2
        return n == int(Fraction(q * x + 2 * r, s * x + 2 * t))

    def prod(s: State, n: int) -> State:
        """
        | 10 -10n | * z -> z
        |  0   1  |
        """
        z, i = s
        return comp((10, -10 * n, 0, 1), z), i

    def cons(s: State, x: LFT) -> State:
        z, i = s
        return comp(z, x), i + 1

    def comp(a: LFT, b: LFT) -> LFT:
        """
        | q r | | u v | = | q*u+r*w q*v+r*x |
        | s t | | w x |   | s*u+t*w s*v+t*x |
        """
        q, r, s, t = a
        u, v, w, x = b
        return q * u + r * w, q * v + r * x, s * u + t * w, s * v + t * x

    init = ((0, 4, 1, 0), 1)
    lfts = ((2 * i - 1, i**2, 1, 0) for i in count(1))
    pi = stream(next, safe, prod, cons, init, lfts)
    for _ in range(n):
        print("".join(map(str, islice(pi, k))))


if __name__ == "__main__":
    n, k, *_ = argv[1:]
    main(int(n), int(k))

Gosper

Code:

from itertools import count, islice
from sys import argv
from typing import Callable, Generator

LFT = tuple[int, int, int, int]
State = tuple[LFT, int]


def main(n: int, k: int):

    def stream(
        next: Callable[[State], int],
        prod: Callable[[State, int], State],
        cons: Callable[[State, LFT], State],
        u: State,
        lfts: Generator[LFT, None, None],
    ) -> Generator[int, None, None]:
        while True:
            [x] = list(islice(lfts, 1))
            v = cons(u, x)
            y = next(v)
            yield y
            u = prod(v, y)

    def next(z: State) -> int:
        [[q, r, s, t], i] = z
        x = 27 * i + 15
        return (q * x + 5 * r) // (s * x + 5 * t)

    def prod(s: State, n: int) -> State:
        """
        | 10 -10n | * z -> z
        |  0   1  |
        """
        z, i = s
        return comp((10, -10 * n, 0, 1), z), i

    def cons(s: State, x: LFT) -> State:
        z, i = s
        return comp(z, x), i + 1

    def comp(a: LFT, b: LFT) -> LFT:
        """
        | q r | | u v | = | q*u+r*w q*v+r*x |
        | s t | | w x |   | s*u+t*w s*v+t*x |
        """
        q, r, s, t = a
        u, v, w, x = b
        return q * u + r * w, q * v + r * x, s * u + t * w, s * v + t * x

    init = ((1, 0, 0, 1), 1)
    lfts = (
        (i * (2 * i - 1), j * (5 * i - 2), 0, j)
        for i, j in ((i, 3 * (3 * i + 1) * (3 * i + 2)) for i in count(1))
    )
    pi = stream(next, prod, cons, init, lfts)
    for _ in range(n):
        print("".join(map(str, islice(pi, k))))


if __name__ == "__main__":
    n, k, *_ = argv[1:]
    main(int(n), int(k))

You may have noticed though, that I replaced the recursive call to stream by an infinite loop.
Well, Python lacks tail call optimization, so we could end up with a stack overflow.

---

Archive:  unbounded-spigot.zip
  Length      Date    Time    Name
---------  ---------- -----   ----
     1790  05-20-2024 11:17   pi-rabinowitz-wagon.py
     1840  05-20-2024 11:17   pi-lambert.py
     1630  05-20-2024 11:17   pi-gosper.py
---------                     -------
     5260                     3 files

[Thomas Klemm]

(05-19-2024 12:34 PM)Guenter Schink Wrote:  Lazy me, simply copied Thomas' code and added some lines for the timing.

We can use that the product of two upper triangular matrices is upper triangular:

\(
\begin{bmatrix}
q & r \\
0 & s
\end{bmatrix}
\begin{bmatrix}
u & v \\
0 & w
\end{bmatrix}
=
\begin{bmatrix}
qu & qv + rw \\
0 & sw
\end{bmatrix}
\)

This saves us some operations and a variable of the state.

And then inline all function calls and use a simple for-loop:

Code:

def main(n, k):
    q, r, s = 1, 0, 1
    line = ""
    for i in range(1, n):
        j = 3 * (3 * i + 1) * (3 * i + 2)
        u, v, w = i * (2 * i - 1), j * (5 * i - 2), j
        # | q r | * | u v | =  | qu  qv + rw |
        # | 0 s |   | 0 w |    |  0       sw |
        q, r, s = q * u, q * v + r * w, s * w
        x = 27 * i + 42
        y = (q * x + 5 * r) // (5 * s)
        line += str(y)
        if i % k == 1:
            print(line)
            line = ""
        # | 10  -10y | * | q r | =  | 10q  10(r - ys) |
        # |  0     1 |   | 0 s |    |   0           s |
        q, r, s = 10 * q, 10 * (r - y * s), s


if __name__ == "__main__":
    n, k = 1_002, 25
    main(n, k)

The digits are printed line by line.
But you may easily change that.

[EdS2]

Thanks Thomas, everyone, some interesting contributions there for sure.

I'm especially interested in Albert's spigot from the ABC documentation - is it described anywhere else, is it attached to any name? It seems simple and effective.

As with Thomas, I like spigots because they produce results incrementally. I think I'm right in saying that the older bounded spigots will produce the first digits slowly and then progressively speed up until their final output. They allocate some fixed (knowable) amount of memory before they begin. The unbounded spigots are the opposite: they produce the first digit immediately and then progressively slow down, allocating (or requiring) ever more memory as they proceed.

(The non-spigot algorithms come in a couple of classes, I think: one class produces the whole result right at the end, the other class creates successive approximations to pi but without necessarily guarantees about how many digits are final in each iteration.)

Other than the time taken to produce some number of digits, perhaps another figure of merit is how many digits (or bytes) of storage is needed. From a very quick look, it seems like the unbounded spigots need around 10 digits per digit produced.

[Thomas Klemm]

(05-20-2024 02:39 PM)EdS2 Wrote:  I'm especially interested in Albert's spigot from the ABC documentation - is it described anywhere else, is it attached to any name? It seems simple and effective.

From the linked document:

Quote:The program works by repeatedly refining a so-called continued fraction that represents pi. The first approximation is 4/1, the second is 4/(1+1/3), which is 12/4. In the next, the 3 is replaced by 3+4/5, and in the next after that, the 5 is replaced by 5+9/7 and so on, so that each part of the fraction is k²/(2k+1+...) for k = 0, 1, 2, ... .

So it appears to use Lambert's formula:

\(
\begin{align}
\pi = \frac{4}{
1 + \frac{1^2}{
3 + \frac{2^2}{
5 + \frac{3^2}{
7 + \cdots}}}}
\end{align}
\)

[EdS2]

Ah, thanks!

[DavidM]

(05-18-2024 01:55 PM)John Keith Wrote:  

(05-18-2024 08:03 AM)EdS2 Wrote:  How fast is a Prime, or an HP49 or HP50, in computing say 1000 digits of pi this way?

Not very fast on the HP 50. I haven't tested up to 1000 digits on a physical calculator but 302 digits (the number returned by Gibbons' original program with input of 1000) takes about 198 seconds. I timed 1000 digits on EMU48 but I'm not at home now and i don't remember the exact time.

Here's some timings for completing 1000 digits of PI:
\begin{array}{|l|r|}
\hline
\textbf{Platform} & \textbf{Time (seconds)} \\
\hline
\text{John's } \textbf{piG3 } \text{on real 50g} & 1960 \\
\text{SysRPL version based on John's } \textbf{piG3} & 1520 \\
\text{Emu48 on Desktop running } \textbf{SysRPL piG3} & 7 \\
\text{LongFloat's } \textbf{FPI } \text{on real 50g} & 43 \\
\text{LongFloat's } \textbf{FPI } \text{on Emu48} & 0.2 \\
\hline
\end{array}

(05-18-2024 01:55 PM)John Keith Wrote:  By way of comparison, the LongFloat function FPI is between 25 and 50 times as fast depending on the number of digits. Gibbons' paper does mention that spigot algorithms are not competitive with state-of-the-art conventional algorithms.

FWIW: the documentation for the LongFloat library attributes the computation of FPI to Werner Huysegom's implementation of Chudnovsky.

[EdS2]

Thanks for the timings!

[EdS2]

There's a handy web page with examples of various spigots in python:
Spigot algorithms for the digits of pi in Python

Quote:Understanding some of the most common spigot algorithms for the digits of pi and their different implementations in Python found across the web.

(Hat-tip to [hoglet] for that)