Pi digits (again) - an unbounded spigot program

[EdS2]

There's a previous mention of Jeremy Gibbons' unbounded spigot not too far away (1, 2) but I don't think I've seen a running version here.

Anyone care to convert the following Basic code, or something close to this JavaScript code, to an HP-71B, or perhaps to a Free42 based calculator, or a Prime, or perhaps a Sharp Pocket Computer?

Strictly, for this spigot to work properly, we need a bignum package, but it turns out to do something interesting even if you just use ordinary floating point numbers. We might get 40 digits of pi on Free42, for example - or we might not!

Because this code will only print some number of correct digits in any given environment, it might be worth reformulating with a FOR NEXT loop that counts out the correct digits and then stops. (Ideally, also print the first incorrect digit too, so we know we've reached the right point.)

Code:

   10 q = 1: r = 0: t = 1: k = 1: n = 3: l = 3
   20 REPEAT
   30   IF (q * 4 + r - t < n * t) ELSE GOTO 90
   40   PRINT ;n;
   50   nr = (r - n * t) * 10
   60   n = INT((q * 3 + r) * 10 / t) - n * 10
   70   q = q * 10: r = nr
   80 UNTIL FALSE
   90 nr = (q * 2 + r) * l
  100 nn = INT((q * k * 7 + 2 + r * l)/(t * l))
  110 q = q * k: t = t * l
  120 l = l + 2: k = k + 1
  130 n = nn: r = nr
  140 UNTIL FALSE

You can see this code running in a BBC Micro context as a one-liner here. (Edit: this link works as expected in Chrome, not in Safari, sorry about that. Not sure about other browsers. If you don't see a pi program running, just copy the text from the box above and paste it in as a replacement to the left pane of owlet.)

[rprosperi]

(05-12-2024 05:03 PM)EdS2 Wrote:  You can see this code running in a BBC Micro context as a one-liner here.

Well, close but not the same program; this runs a multi-colored HELLO WORLD, but I the way you can send the actual code to the site via the link is clever!

[EdS2]

Bother - if the owlet link doesn't work for you (it works for me) please paste the program into the left pane.

Edit: or perhaps try https://bbcmic.ro/?t=9SPqt

[Thomas Klemm]

(05-12-2024 07:43 PM)EdS2 Wrote:  if the owlet link doesn't work for you

It is a well known issue that MyBB breaks long links: Image hyperlink insertion - LTR marks (?) automatically inserted then http 404

[rprosperi]

It works right with the new link, but it halts after 31459265 with "Too big at line 10", where it exceeded what the old machine could handle. Funny that with all the commands in a single line, that error message is not really so helpful, lol...

[Thomas Klemm]

(05-12-2024 05:03 PM)EdS2 Wrote:  Anyone care to convert the following Basic code, or something close to this JavaScript code, to an HP-71B, or perhaps to a Free42 based calculator, or a Prime, or perhaps a Sharp Pocket Computer?

Cf. (35S) spigot algorithm for the digits of PI

---

Cf. Many Digits of Pi (HP32sii, HP-16c, HP-12c, HP-12cp, HP-19c, HP-30b, HP-15c)

Quote:based on:
-------------------------------------------------------------------------
pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)]
= 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......

I could be wrong, but I assume that Katie Wasserman used the same spigot algorithm.

[EdS2]

Thanks Thomas! I think that's the earlier Rabinowitz (and Wagon) bounded spigot (although I'm not sure) - this on the other hand is Gibbons' unbounded spigot. Different, in some sense, except that we don't have bignums here. So it will in fact come out less capable than the bounded spigot...

[toml_12953]

(05-12-2024 05:03 PM)EdS2 Wrote:  Anyone care to convert the following Basic code, or something close to this JavaScript code, to an HP-71B, or perhaps to a Free42 based calculator, or a Prime, or perhaps a Sharp Pocket Computer?

Code:

   10 q = 1: r = 0: t = 1: k = 1: n = 3: l = 3
   20 REPEAT
   30   IF (q * 4 + r - t < n * t) ELSE GOTO 90
   40   PRINT ;n;
   50   nr = (r - n * t) * 10
   60   n = INT((q * 3 + r) * 10 / t) - n * 10
   70   q = q * 10: r = nr
   80 UNTIL FALSE
   90 nr = (q * 2 + r) * l
  100 nn = INT((q * k * 7 + 2 + r * l)/(t * l))
  110 q = q * k: t = t * l
  120 l = l + 2: k = k + 1
  130 n = nn: r = nr
  140 UNTIL FALSE

I can't even run this in JavaScript. I get an error on the q in line 10:

Expected ";" but found "q"

[toml_12953]

(05-12-2024 05:03 PM)EdS2 Wrote:  Anyone care to convert the following Basic code, or something close to this JavaScript code, to an HP-71B, or perhaps to a Free42 based calculator, or a Prime, or perhaps a Sharp Pocket Computer?

Code:

   10 q = 1: r = 0: t = 1: k = 1: n = 3: l = 3
   20 REPEAT
   30   IF (q * 4 + r - t < n * t) ELSE GOTO 90
   40   PRINT ;n;
   50   nr = (r - n * t) * 10
   60   n = INT((q * 3 + r) * 10 / t) - n * 10
   70   q = q * 10: r = nr
   80 UNTIL FALSE
   90 nr = (q * 2 + r) * l
  100 nn = INT((q * k * 7 + 2 + r * l)/(t * l))
  110 q = q * k: t = t * l
  120 l = l + 2: k = k + 1
  130 n = nn: r = nr
  140 UNTIL FALSE

Here's the code in ANSI BASIC:

Code:

100 DECLARE EXTERNAL SUB digits
110 CALL digits
120 END
130 !
140 EXTERNAL SUB digits
150 LET q = 1
160 LET r = 0
170 LET t = 1
180 LET k = 1
190 LET n = 3
200 LET l = 3
210 DO
220    IF q * 4 + r - t < n * t THEN
230       PRINT n;
240       LET nr = (r - n * t) * 10
250       LET n = INT((q * 3 + r) * 10 / t) - n * 10
260       LET q = q * 10
270       LET r = nr
280    ELSE
290       LET nr = (q * 2 + r) * l
300       LET nn = INT((q * k * 7 + 2 + r * l) / (t * l))
310       LET q = q * k
320       LET t = t * l
330       LET l = l + 2
340       LET k = k + 1
350       LET n = nn
360       LET r = nr
370    END IF
380 LOOP
390 END SUB

[EdS2]

I fired up an HP-71B in DOSbox and typed the program in...
... result is 3.14159265358979...
with further digits being inaccurate.

   

[toml_12953]

Here's a translation of a FORTRAN program. It will give you as many digits as you want with no bignum needed. It runs in a small memory machine with no large array needed. Just increase the 50000 in the FOR statement to increase the number of digits. Be warned, though, it takes quite a while to calculate 200,000 digits!

Code:

DIM vect(3350)
SET MARGIN 55
LET more = 0
LET total=0
MAT vect = 2 * CON  ! Fill array vect with the number 2 in all positions.
FOR n = 1 TO 50000
   LET karray = 0
   FOR l = 3350 TO 1 STEP -1
      LET num = 100000*vect(l) + karray*l
      LET karray = INT(num/(2*l - 1))
      LET vect(l) = num - karray*(2*l - 1)
   NEXT l
   LET k = INT(karray/100000)
   WHEN EXCEPTION IN  ! If number is too big to print in 5 positions, use 6.
      PRINT USING$("%%%%%",more+k);    ! print more+k with leading zeros
   USE
      PRINT USING$("%%%%%%",more+k);    ! print more+k with leading zeros
   END WHEN
   LET total=total+5
   LET more = karray - k*100000
NEXT n
PRINT
PRINT total-5;"decimal places"
END

[EdS2]

(05-13-2024 07:21 PM)toml_12953 Wrote:  Here's the code in ANSI BASIC:
...

Thanks - just converted back to BBC Basic and it checks out!

(05-14-2024 01:45 PM)toml_12953 Wrote:  Here's a translation of a FORTRAN program. It will give you as many digits as you want with no bignum needed. It runs in a small memory machine with no large array needed.
...

Hmm - that's a Rabinowitz spigot, I think, and it uses an array for the base conversion. Granted, I expect it's much more memory efficient than the Gibbons spigot, but I haven't tried to measure that. (Bearing in mind, the point of my initial post here was to show an unbounded spigot, rather than the usual bounded spigot.)

[John Keith]

(05-12-2024 08:57 PM)EdS2 Wrote:  ... this on the other hand is Gibbons' unbounded spigot. Different, in some sense, except that we don't have bignums here.

We do if we are using the HP 49/50 or Prime. The spigot programs use only integer math and thus do not require arbitrary-precision floating-point software, making them a good fit for the newer HP's. Here is a straightforward RPL implementation. The program stores the digits in a string for convenience but could be easily adapted to display each digit if desired. Given an integer s on the stack, the program will return approximately 0.3*s digits. This is because only about 1 in 4 iterations actually return a new digit.

Code:


\<< 1 0 1 1 3 3 \-> s q r t k n l
  \<< ""                                   @ Empty string on stack
    WHILE k s <
    REPEAT
      IF q 4 * r + t - n t * <
      THEN n +                             @ Append n to string
        r n t * - 10 *                     @ nr
        q 3 * r + 10 * t IQUOT
        n 10 * - 'n' STO
        'q' 10 STO* 'r' STO
      ELSE q 2 * r + l *                   @ nr
        q k * 7 * r l * + 2 + t l * IQUOT  @ nn
        'q' k STO* 't' l STO*
        'l' 2 STO+ 'k' 1 STO+
        'n' STO 'r' STO
      END
    END
  \>>
\>>

If we assume that the conjecture in section 6 of Gibbons' paper is correct, we can use the last program in the paper called 'piG3' which is much more efficient. Not only does it return a new digit for every iteration, it also has features that make the program run faster. It only requires one integer divide per digit, and those are slow in RPL. Also, since there is no test clause, the control structure is simpler and we can replace the WHILE loop with a FOR loop which is faster. Here, then, is an RPL implementation of 'piG3'. It runs about 7 times faster than the previous program and is about 100 bytes smaller. Given an integer on the stack, the program will return that many digits.

Code:


\<< 1 180 60 1 1 \-> q r t u y
  \<< ""                                       @ Empty string on stack
    2 ROT 1 +                                  @ k from 2 to n+1
    FOR k k 27 * 12 - q * r 5 * + t 5 * IQUOT  @ y
      DUP 'y' STO +                            @ Append y to string
      k 3 * 1 + DUP 1 + * 3 * 'u' STO          @ u
      q k * 10 * k 2 * 1 - *                   @ Pre-compute q
      k 5 * 2 - q * r + y t * - 10 * u *       @ r
      'r' STO 'q' STO 't' u STO*               @ Update r, q, t
    NEXT
  \>>
\>>

[EdS2]

Thanks John, that's splendid! I never had a clue how to rewrite Haskell into something I could understand.

I quite like a pi program which needs an unproved conjecture (but is known to be correct for at least a thousand digits (Edit: and I've since verified it to 100k digits) (Edit: and now to a million digits))

In BBC Basic:

Code:

BBC BASIC for MacOS Console v0.46
(C) Copyright R. T. Russell, 2024
>LIST
   10 REM unbounded pi spigot, Gibbons piG3
   20 REM thanks to John Keith on hpmuseum
   30 q = 1: r = 180: t = 60: u = 1: y = 1
   40 FOR k=2 TO 30
   50   y = INT((q * (27*k-12) + (5*r))/(5*t))
   60   PRINT ;y;
   70   u = 3 * (3*k+1) * (3*k+2)
   80   r = 10*u*(q*(5*k-2) + r - y*t)
   90   q = 10*q*k * (2*k-1)
  100   t = t*u
  110 NEXT
>RUN
31415926535897932384626433812>

In Owlet (an 8 bit Basic with 5 byte floats) we get

Code:

314159265358
Too big at line 80

which is rather nicely more digits than can fit in a float.

How fast is a Prime, or an HP49 or HP50, in computing say 1000 digits of pi this way?

[Albert Chan]

(12-12-2021 11:04 PM)Albert Chan Wrote:  PI = 4/(1+1/(3+2^2/(5+3^2/(7+4^2/(9+5^2/(11 ...

This keep integers *much* smaller ⇒ code run faster.
see ABC Programmer's Handbook, Pi example

Translated for HP71B

10 P=1 @ Q=3 @ A=4 @ B=1 @ A1=12 @ B1=4
20 P=P+Q @ Q=Q+2 @ A0=A @ B0=B @ A=A1 @ B=B1 @ A1=P*A0+Q*A1 @ B1=P*B0+Q*B1
30 D=A DIV B @ R=A1-D*B1 @ IF R<0 OR R>=B1 THEN 20
40 DISP D; @ A=10*(A-D*B) @ A1=10*R @ GOTO 30

>run
3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 5 6 2 ...

12 decimals float give 20 good pi digits.

[Thomas Klemm]

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

Here's the Haskell program for Gosper's series rewritten in Python:

Code:

def main(n, k):
    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 range(1, n))
    )
    pi = stream(init, lfts)
    i = 0
    for d in pi:
        print(d, end="")
        if i % k == 0:
            print()
        i += 1

def stream(u, lfts):
    for x in lfts:
        v = cons(u, x)
        y = next(v)
        yield y
        u = prod(v, y)

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

def prod(s, n):
    z, i = s
    return comp((10, -10 * n, 0, 1), z), i

def cons(s, x):
    z, i = s
    return comp(z, x), i + 1

def comp(a, b):
    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

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

It runs with MicroPython. Therefore, it should run on the Prime as well.

(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?

I don't have a Prime to test it, but it takes 0.053s on a MacBook M2 Pro:

Code:

3
1415926535897932384626433
8327950288419716939937510
5820974944592307816406286
2089986280348253421170679
8214808651328230664709384
4609550582231725359408128
4811174502841027019385211
0555964462294895493038196
4428810975665933446128475
6482337867831652712019091
4564856692346034861045432
6648213393607260249141273
7245870066063155881748815
2092096282925409171536436
7892590360011330530548820
4665213841469519415116094
3305727036575959195309218
6117381932611793105118548
0744623799627495673518857
5272489122793818301194912
9833673362440656643086021
3949463952247371907021798
6094370277053921717629317
6752384674818467669405132
0005681271452635608277857
7134275778960917363717872
1468440901224953430146549
5853710507922796892589235
4201995611212902196086403
4418159813629774771309960
5187072113499999983729780
4995105973173281609631859
5024459455346908302642522
3082533446850352619311881
7101000313783875288658753
3208381420617177669147303
5982534904287554687311595
6286388235378759375195778
1857780532171226806613001
9278766111959092164201989
python micro-pi.py  0.04s user 0.01s system 88% cpu 0.053 total

[John Keith]

(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.

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.

[Guenter Schink]

How fast is PRIME?

I copied Thomas' Python code to the Prime. Some adjustments for timing:

Code:

from hpprime import *

 
def main(n, k):
    eval('PRINT()')
    a=eval('time')
    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 range(1, n))
    )
    pi = stream(init, lfts)
    i = 0
    for d in pi:
        print(d, end="")
        if i % k == 0:
            print()
        i += 1
    a=round(eval('time')-a,4)
    print("\n",a)

def stream(u, lfts):
    for x in lfts:
        v = cons(u, x)
        y = next(v)
        yield y
        u = prod(v, y)

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

def prod(s, n):
    z, i = s
    return comp((10, -10 * n, 0, 1), z), i

def cons(s, x):
    z, i = s
    return comp(z, x), i + 1

def comp(a, b):
    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

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

output is the same as Thomas'
time is between 5 and 6.5 seconds with several runs

Günter

[Thomas Klemm]

(05-18-2024 03:19 PM)Guenter Schink Wrote:  How fast is PRIME?
(…)
time is between 5 and 6.5 seconds with several runs

Not sure if it helps, but I remembered Writing fast and efficient MicroPython and made the functions local to main:

Code:

def main(n, k):

    def stream(u, lfts):
        for x in lfts:
            v = cons(u, x)
            y = next(v)
            yield y
            u = prod(v, y)

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

    def prod(s, n):
        z, i = s
        return comp((10, -10 * n, 0, 1), z), i

    def cons(s, x):
        z, i = s
        return comp(z, x), i + 1

    def comp(a, b):
        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 range(1, n))
    )
    pi = stream(init, lfts)
    i = 0
    for d in pi:
        print(d, end="")
        if i % k == 0:
            print()
        i += 1


if __name__ == "__main__":
    n, k = 1002, 100
    main(n, k)

[komame]

(05-18-2024 03:19 PM)Guenter Schink Wrote:  output is the same as Thomas'
time is between 5 and 6.5 seconds with several runs

If performance is important, this code can be optimized by buffering the result in a variable and then outputting it to the terminal all at once, instead of outputting partial results multiple times.

Code:

from hpprime import *

def main(n, k):
    eval('PRINT')
    s = ''
    a=eval('time')
    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 range(1, n))
    )
    pi = stream(init, lfts)
    i = 0
    for d in pi:
        s += str(d)
        if i % k == 0:
            s += '\n'
        i += 1
    a=round(eval('time')-a,4)
    print(s)
    print("\n",a)

def stream(u, lfts):
    for x in lfts:
        v = cons(u, x)
        y = next(v)
        yield y
        u = prod(v, y)

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

def prod(s, n):
    z, i = s
    return comp((10, -10 * n, 0, 1), z), i

def cons(s, x):
    z, i = s
    return comp(z, x), i + 1

def comp(a, b):
    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

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

This reduces the execution time of this task on the HP Prime to below 1 second (significant acceleration).