Happy Pi Day!

[Gerson W. Barbosa]

Let's celebrate it by computing the first two hundred digits of pi using a very unusual formula 
for that purpose, the Wallis Product. Normally it would not be suitable for this application, as
one million terms yield only five decimal digits of pi. However, a correction term in continued
fraction form will handle that nicely:


                                  /                                         \
       256 196 144 100 64 36 16 4 |                      -1                 |
pi ~ 2.---.---.---.---.--.--.--.-.| 1 - ----------------------------------- |
       255 195 143  99 63 35 15 3 |                         3               |
                                  |     32 + ------------------------------ |
                                  |                           1             |
                                  |          2 - -------------------------- |
                                  |                             5           |
                                  |              32 + --------------------- |
                                  |                               3         |
                                  |                   2 - ----------------- |
                                  |                                 7       |
                                  |                       32 + ------------ |
                                  |                                   5     |
                                  |                            2 - -------- |
                                  |                                     9   |
                                  |                               32 + ---  |
                                  |                                     2   |
                                  \                                         /



      349075614466048
pi ~ ----------------- ~ 3.1415926535(708367)
      111114219110895

Since we have 8 terms of the Wallis Product one of the constants in the continued fraction is 32, that is,

4 times 8. The other constant is always 2. The alternate numerators follow a simple pattern (-1, 1, 3, 5...
and 3, 5, 7, 9...). 

For 200 digits, we will need only 150 terms of the Wallis Products plus another 150 terms of the continued
fraction.


It will take about 9 and half minutes on the real HP 50g:

200

%%HP: T(3)A(R)F(.);
\<< PUSH RAD -105 CF -3 CF DUP 3 * 4 + 8 IDIV2 DROP DUP + DUP 4 OVER * SWAP 1 + DUP 4 - 0 1 \-> d n1 n2 c p
  \<< 2 / 1 SWAP
    FOR i i DUP 1 - OVER 256 * * 64 + DUP PICK3 SQ * UNROT 32 - OVER * 16 + * 3 - / p * 'p' STO n2 d n1 2 c 
      - / + / 'c' STO 'n2' 2 STO- 'n1' 2 STO-
    NEXT p 1 c - * 2 * EXPAND
  \>> FXND DUP SIZE R\->I ALOG OVER - PICK3 * SWAP IQUOT + \->STR DUP HEAD 0 I\->R \->STR TAIL + SWAP TAIL + 
  1 ROT 2 + SUB POP
\>>

EVAL ->

3.14159265358979323846264338327950288419716939937510
  58209749445923078164062862089986280348253421170679
  82148086513282306647093844609550582231725359408128
  48111745028410270193852110555964462294895493038196


Notice the continued fraction correction factor has been obtained empirically, so no proof is available.

Edited to fix a typo.

[pinkman]

Happy to read this here today!

[Gerson W. Barbosa]

314 decimal places might be a more adequate number of places, but it takes 3613 seconds on the HP 50g, or roughly one hour. 100 places takes only 70 seconds, but that's too few digits -- there are people who know those by heart -- so I chose 200, which takes reasonable 10 minutes. Well, reasonable time considering we are using Wallis Product here.

For more digits, you might want to try the algorithm on programming languages that have arbitrary precision extension. Here are 997 digits of pi using Decimal BASIC, which is limited to 1000 digits:

Code:


OPTION ARITHMETIC DECIMAL_HIGH
INPUT  PROMPT "number of digits: ":nd
LET tm = TIME 
LET n = 2*INT((3*nd + 4)/8)
PRINT "n = ";n
LET pr = 1
LET n1 = n + 1
LET n2 = n1 - 4
LET d = 4*n
LET c = 0
FOR i = 1 TO n/2
   LET t = 256*i*(i - 1) + 64
   LET pr = pr*i*i*t/(i*(i*(t - 32) + 16) - 3)
   LET c = n2/(d + n1/(2 - c))
   LET n2 = n2 - 2
   LET n1 = n1 - 2
next i
LET c = 2 - 2*c
LET p = pr*c
LET r = TIME - tm
LET r$ = STR$(p - INT(p))
PRINT
PRINT STR$(INT(p));
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT REPEAT$(" ",1 + LEN(STR$(INT(p))));
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "##.##": r;
PRINT " seconds"
END

number of digits: 997
n =  748 

3.1415926535 8979323846 2643383279 5028841971 6939937510 
  5820974944 5923078164 0628620899 8628034825 3421170679 
  8214808651 3282306647 0938446095 5058223172 5359408128 
  4811174502 8410270193 8521105559 6446229489 5493038196 
  4428810975 6659334461 2847564823 3786783165 2712019091 
  4564856692 3460348610 4543266482 1339360726 0249141273 
  7245870066 0631558817 4881520920 9628292540 9171536436 
  7892590360 0113305305 4882046652 1384146951 9415116094 
  3305727036 5759591953 0921861173 8193261179 3105118548 
  0744623799 6274956735 1885752724 8912279381 8301194912 
  9833673362 4406566430 8602139494 6395224737 1907021798 
  6094370277 0539217176 2931767523 8467481846 7669405132 
  0005681271 4526356082 7785771342 7577896091 7363717872 
  1468440901 2249534301 4654958537 1050792279 6892589235 
  4201995611 2129021960 8640344181 5981362977 4771309960 
  5187072113 4999999837 2978049951 0597317328 1609631859 
  5024459455 3469083026 4252230825 3344685035 2619311881 
  7101000313 7838752886 5875332083 8142061717 7669147303 
  5982534904 2875546873 1159562863 8823537875 9375195778 
  1857780532 1712268066 1300192787 6611195909 2164201

Runtime:   .47 seconds

Edited to fix a typo

PS: Some algorithm and program optimizations

-----------------------------------

HP 50g:

314

%%HP: T(3)A(R)F(.);
\<< PUSH RAD -105 CF -3 
CF 1 DUP PICK3 3 * 8 
+ 8 IDIV2 DROP DUPDUP 
+ 1 + OVER 8 * 0 \-> m 
d c
  \<<
    FOR i i DUP 16 * 
8 - * DUP SQ SWAP DUP 
2 - * 3 - / * m 4 - d 
m 2 c - / + / 'c' STO 
'm' 2 STO-
    NEXT 2 * 1 c - *
  \>> EXPAND FXND DUP 
SIZE R\->I ALOG OVER - 
PICK3 * SWAP IQUOT + 
\->STR DUP HEAD 0 I\->R 
\->STR TAIL + SWAP TAIL 
+ 1 ROT 2 + SUB POP
\>>

TEVAL ->

3.14159265358979323846264338327950288419716939937510
  58209749445923078164062862089986280348253421170679
  82148086513282306647093844609550582231725359408128
  48111745028410270193852110555964462294895493038196
  44288109756659334461284756482337867831652712019091
  45648566923460348610454326648213393607260249141273
  72458700660631

s:3265.5887 (real HP 50g)

337 bytes, CK = # F36Bh


-----------------------------------

TurboBCD:

Code:


Program Wallis;

var i, k, m, n: integer;
    c, d, t, w: real;

begin
  ClrScr;
  WriteLn('+---+---------------------+---------------------+');
  WriteLn('| N |         2*W         |        2*W*C        |');
  WriteLn('+---+---------------------+---------------------+');
  for k := 1 to 7 do
    begin
      n := 2*k;
      w := 1;
      c := 0;
      d := 4*n;
      m := n + 1;
      for i := 1 to k do
        begin
          t := i*(16*i - 8);
          w := w*t*t/(t*(t - 2) - 3);
          c := (m - 4)/(d + m/(2 - c));
          m := m - 2
        end;
      c := 1 - c;
      WriteLn('|',n:2,' |',2*w:20:17,' |',2*w*c:20:17,' |')
    end;
  WriteLn('+---+---------------------+---------------------+')
end.

+---+---------------------+---------------------+
| N |         2*W         |        2*W*C        |
+---+---------------------+---------------------+
| 2 | 2.84444444444444444 | 3.14385964912280701 |
| 4 | 2.97215419501133786 | 3.14158816337302932 |
| 6 | 3.02317019200136082 | 3.14159266276745771 |
| 8 | 3.05058999605551092 | 3.14159265357083669 |
|10 | 3.06770380664349896 | 3.14159265358983256 |
|12 | 3.07940134316788626 | 3.14159265358979314 |
|14 | 3.08790206983111306 | 3.14159265358979321 |
+---+---------------------+---------------------+

-----------------------------------

Decimal BASIC:

Code:


OPTION ARITHMETIC DECIMAL_HIGH
INPUT  PROMPT "number of digits: ":nd
LET tm = TIME 
LET n = 2*INT((3*nd + 8)/8)  
POPTION ARITHMETIC DECIMAL_HIGH
INPUT  PROMPT "number of digits: ":nd
LET tm = TIME 
LET n = 2*INT((3*nd + 8)/8)  
PRINT "n = ";n
LET w = 1
LET c = 0
LET d = 4*n
LET m = n + 1
FOR i = 1 TO n/2
   LET t = i*(16*i - 8)
   LET w = w*t*t/(t*(t - 2) - 3)
   LET c = (m - 4)/(d + m/(2 - c))
   LET m = m - 2
next i
LET c = 1 - c
LET p = 2*w*c
LET r = TIME - tm
LET r$ = STR$(2*w*c - INT(2*w*c))
PRINT
PRINT STR$(INT(p));
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT REPEAT$(" ",1 + LEN(STR$(INT(p))));
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "##.##": r;
PRINT " seconds"
END

[Dave Britten]

Last month: "Darn, I can't wear my Pi Day t-shirt to work because it'll be on Saturday."

This month: "Darn, I can't wear my Pi Day t-shirt to work because we're in the middle of the apocalypse."

[Gerson W. Barbosa]

(03-14-2020 06:19 AM)pinkman Wrote:  Happy to read this here today!

Yes, let’s not forget this is ultimately the result of a conversation you and your cousin started one month ago. I don’t think any of my brothers and cousins would ever engage in such inspiring conversations (most of them lawyers) :-)

https://www.hpmuseum.org/forum/thread-14476.html

Thanks again!

Gerson.

[rprosperi]

(03-14-2020 02:53 PM)Gerson W. Barbosa Wrote:  

(03-14-2020 06:19 AM)pinkman Wrote:  Happy to read this here today!

Yes, let’s not forget this is ultimately the result of a conversation you and your cousin started one month ago. I don’t think any of my brothers and cousins would ever engage in such inspiring conversations (most of them lawyers) :-)

https://www.hpmuseum.org/forum/thread-14476.html

Thanks again!

Gerson.

I'm so sorry Gerson... sad news indeed...

On the other hand, if I was to discuss the number of digits in PI with my siblings, they would no doubt assume I was talking about eating pie with my fingers.

[Jim Horn]

(03-14-2020 05:04 PM)rprosperi Wrote:  I'm so sorry Gerson... sad news indeed...

On the other hand, if I was to discuss the number of digits in PI with my siblings, they would no doubt assume I was talking about eating pie with my fingers.

As a Gyro vendor's sign once noted, "Fingers are better eaten separately."

I'll spare the readership my 22 Pi-kus...

[BruceH]

On a related note, those nice people at Wolfram sent me a link to a special offer - 31.41% discount on a Mathematica Home/Individual licence just because it's Pi day. Unfortunately the link doesn't work. :-(

I've sent them an email - I'll let you know if it eventually happens.

[Eddie W. Shore]

Happy Pi Day Plus One! (or really Happy Pi Day Rounded To Two Decimal Places Plus One Hundredth)

[DM48]

Gerson, below is part of a python spigot algo my son uses to calculate PI. The range 10,000,000 would need to be reduced a lot for the HP48G. Any thoughts on the best way to program this up in RPL?


q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
for j in range(10,000,000):
if 4 * q + r - t < m * t:
yield m
q, r, t, k, m, x = 10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x
else:
q, r, t, k, m, x = q*k, (2*q+r)*x, t*x, k+1, (q*(7*k+2)+r*x)//(t*x), x+2

[Albert Chan]

Hi, DM48

A better PI spigot is with continued fraction:

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 (*)

Below Python code produce 7656 pi digits, in 6 sec on my laptop.
To produce same output, your code (adjusted to range(33080)) take 31 sec.
(It take about 2 sec to print digits, speed ratio ≈ 29/4 ≈ 7)

Code:

p, q = 1, 3
a, b, a1, b1 = 4, 1, 12, 4
for i in range(10000):
    p, q = p+q, q+2        # next convergent
    a, b, a1, b1 = a1, b1, p*a+q*a1, p*b+q*b1
    while 1:
        d = a//b
        r = a1 - d*b1      # break if a//b != a1//b1
        if r < 0 or r >= b1: break
        print d,           # common digits
        a, a1 = 10*(a-d*b), 10*r

(*) Python is descendant of ABC, both developed by Guido van Rossum

[Gerson W. Barbosa]

(12-12-2021 06:44 PM)DM48 Wrote:  Gerson, below is part of a python spigot algo my son uses to calculate PI. The range 10,000,000 would need to be reduced a lot for the HP48G. Any thoughts on the best way to program this up in RPL?
...

Hi,

Pi Day 2022 is still a bit away, but this is not a problem. The problem is I do not understand Python. Anyway, I may give it a try later if I figure out how it works.
In Pi Day 2009 Egan Ford presented a comprehensive collection of programs to calculate Pi, including "HP-71B Spigot Forth" (1000 digits in 6720.22 seconds). Perhaps it should be easier to convert that one to RPL.

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=148399

Regarding the formula in the OP, I've found a nicer one. An RPL program is also offered in the vein "it's slow but it's mine" :-) (357 decimal digits in 3820.207 seconds)

[Gil]

A nice program I copied somewhere — possibly in this forum — in my HP50G :

\<< -105 CF 0 \-> N R
\<< 280000 N ALOG * DUP 'R' STO 2 \-> y
\<<
DO y * y 1 + IQUOT 50 IQUOT DUP 'R' STO+ 2 'y' STO+
UNTIL DUP 0 ==
END DROP
\>> 30336 N ALOG * DUP 'R' STO+ 2 \-> y
\<<
DO 9 * y * 6250 IQUOT y 1 + IQUOT DUP 'R' STO+ 2 'y' STO+
UNTIL DUP 0 ==
END DROP R
\>>
\>>
\>>

Fast and works fine.

[DM48]

(12-13-2021 12:41 AM)Gerson W. Barbosa Wrote:  

(12-12-2021 06:44 PM)DM48 Wrote:  Gerson, below is part of a python spigot algo my son uses to calculate PI. The range 10,000,000 would need to be reduced a lot for the HP48G. Any thoughts on the best way to program this up in RPL?
...

Hi,

Pi Day 2022 is still a bit away, but this is not a problem. The problem is I do not understand Python. Anyway, I may give it a try later if I figure out how it works.
In Pi Day 2009 Egan Ford presented a comprehensive collection of programs to calculate Pi, including "HP-71B Spigot Forth" (1000 digits in 6720.22 seconds). Perhaps it should be easier to convert that one to RPL.

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=148399

Regarding the formula in the OP, I've found a nicer one. An RPL program is also offered in the vein "it's slow but it's mine" :-) (357 decimal digits in 3820.207 seconds)

Albert's Python script may be better to use since he found a significant speed increase. Thanks Albert!

I am assuming the HP42s / DM42/ Free 42 wouldn't be appropriate for this type of calculation due to the nature of storing a lot of digits of Pi. Is that a good assumption Albert?

Pi day is a ways off still, but we can't be too prepared. :-)

Gil, thank you for the RPL program. I am going to input this into my 48G and see what happens. I see at least one instruction I am not familiar with. Hopefully, it works. I love to run these on my phone in iHP48 inside a 48GX to see just how fast they can go.

[Gerson W. Barbosa]

(12-13-2021 01:21 AM)Gil Wrote:  A nice program I copied somewhere — possibly in this forum — in my HP50G :
…

Yes, from this very old thread:

https://www.hpmuseum.org/cgi-sys/cgiwrap...read=53988 (Message #17).