Π day

[robve]

Happy \( \pi\ {\cal D}\alpha\gamma ! \)

Sharp's series

\( \frac{\pi}{6} = \frac{1}{\sqrt3}\left(1-\frac{1}{3^1\cdot3}+\frac{1}{3^2\cdot5}-\frac{1}{3^3\cdot7}+\cdots\right) \)

on a Sharp:

10 S=0,P=1
20 FOR N=0 TO 16
30 S=S+P/(3^N*(2*N+1)),P=-P
40 NEXT N
50 PRINT "Π=";6*S/√3
RUN
Π=3.141592654

and on a very nice non-Sharp:

10 DESTROY ALL @ S=0 @ P=1
20 FOR N=0 TO 20
30 S=S+P/(3^N*(2*N+1)) @ P=-P
40 NEXT N
50 DISP "PI=";6*S/SQR(3)
RUN
PI=3.14159265359

What's your favorite pi series?

- Rob

[Dave Britten]

Personally I kind of like the absurdity of using the Monte Carlo method.

https://www.wikihow.com/Calculate-Pi-by-...n-Hot-Dogs

[Gerson W. Barbosa]

(03-14-2022 03:35 AM)robve Wrote:  Happy \( \pi\ {\cal D}\alpha\gamma ! \)

Sharp's series

Sharp's series

\( \frac{\pi}{6} = \frac{1}{\sqrt3}\left(1-\frac{1}{3^1\cdot3}+\frac{1}{3^2\cdot5}-\frac{1}{3^3\cdot7}+\cdots\right) \)

on a Sharp:

10 S=0,P=1
20 FOR N=0 TO 16
30 S=S+P/(3^N*(2*N+1)),P=-P
40 NEXT N
50 PRINT "?=";6*S/√3
RUN
?=3.141592654

That’s a nice one. Interestingly, it was first derived in India by Madhava some three centuries earlier:

https://en.wikipedia.org/wiki/Madhava_of..._of_π_(pi)


Happy \(\pi\) Day you all!

[Gerson W. Barbosa]

(03-14-2022 11:44 AM)Dave Britten Wrote:  Personally I kind of like the absurdity of using the Monte Carlo method.

https://www.wikihow.com/Calculate-Pi-by-...n-Hot-Dogs

That’s literally the wurst method to compute \(\pi\), but it’s an interesting one nonetheless.

[EdS2]

> Sharp's series

It is a nice one, and a nice connection between Sharp and Sharp.

> Happy π Day to you all!

And you!

Bit of a spoiler, but here's a related quiz:
Pick a random number between 0 and 1, then pick another and divide the first by the second. Round off the result to an integer. What's the probability that it's even?

> What's your favorite pi series?

There's an unusual one by Valentin... 3 years ago, which seems appropriate.

It feels worthwhile also to showcase one of Gerson's [but see downthread!] pandigital approximations ("good to 17 digits"):

That is: (ln{[2×5!+(8-1)!]^(sqrt(9))+4!+(3!)!})/(sqrt(67))

[robve]

Π day really needs some RPN/FORTH too. So here it goes with FORTH500:

: Sharp ( iter -- )
  0e
  0 DO
    I 1 AND 2* 1- NEGATE S>F
    3e I S>F F** I 2* 1+ S>F F* F/ F+
  LOOP
  6e F* 3e FSQRT F/ ." PI=" F.
;
24 Sharp
PI=3.1415926535897

Let's use double precision float to produce pi in glorious 20 decimals:

: Sharp ( iter -- )
  0d
  0 DO
    I 1 AND 2* 1- NEGATE S>F
    3d I S>F F** I 2* 1+ S>F F* F/ F+
  LOOP
  6d F* 3d FSQRT F/ ." PI=" F.
;
38 Sharp
PI=3.1415926535897932385

Edit note for clarification: FORTH500 has a separate float stack, like many 80s Forth's do, which means that 0 DO takes the loop limit from the stack, i.e. the value passed to the Sharp word.

- Rob

[robve]

(03-14-2022 12:25 PM)Gerson W. Barbosa Wrote:  

(03-14-2022 11:44 AM)Dave Britten Wrote:  Personally I kind of like the absurdity of using the Monte Carlo method.

https://www.wikihow.com/Calculate-Pi-by-...n-Hot-Dogs

That’s literally the wurst method to compute \(\pi\), but it’s an interesting one nonetheless.

Worst means hot dog (or generally sausage) in Dutch. So yeah, definitely the worst way to compute \( \pi \) in Anglo-Dutch. I tried it, but after one throw the dogs grabbed the hot dogs. Not even one digit of \( \pi \). Comte de Buffon would not approve.

Perhaps the worst way to compute \( \pi \) is to exploit an interesting relationship between (co-)primes and \( \pi \). If we take \( n \) random integer pairs and determine the fraction \( R_n \) that are co-prime, then we have \( \pi \approx \sqrt{\frac{6}{R_n}} \) where \( \displaystyle \lim_{n\rightarrow\infty} \sqrt{\frac{6}{R_n}} = \pi \). It will take one million random pairs to get about 1 or 2 decimal places of \( \pi \) and 100,000,000 random pairs to get 5 decimal places.

- Rob

[Gerson W. Barbosa]

(03-14-2022 01:55 PM)EdS2 Wrote:  .

It feels worthwhile also to showcase one of Gerson's pandigital approximations ("good to 17 digits"):

That is: (ln{[2×5!+(8-1)!]^(sqrt(9))+4!+(3!)!})/(sqrt(67))

Thanks you, but it’s not mine. Actually Ramanujan’s, I think. I’ve only written it a (not so elegant) pandigital form. I once asked Eric to add a remark, but apparently he didn’t read my message.

[Dave Britten]

(03-14-2022 12:25 PM)Gerson W. Barbosa Wrote:  

(03-14-2022 11:44 AM)Dave Britten Wrote:  Personally I kind of like the absurdity of using the Monte Carlo method.

https://www.wikihow.com/Calculate-Pi-by-...n-Hot-Dogs

That’s literally the wurst method to compute \(\pi\), but it’s an interesting one nonetheless.

Have you ever sausage ingenuity, though?

[Gerson W. Barbosa]

(03-14-2022 02:10 PM)robve Wrote:  Π day really needs some RPN/FORTH too.

Not the fastest algorithm in town, quite the contrary! And yet faster than this BASIC equivalent by two orders of magnitude.

[Thomas Klemm]

(03-14-2022 03:35 AM)robve Wrote:  \( \frac{\pi}{6} = \frac{1}{\sqrt3}\left(1-\frac{1}{3^1\cdot3}+\frac{1}{3^2\cdot5}-\frac{1}{3^3\cdot7}+\cdots\right) \)

That's an alternating series so we can use Valentin's program:
(11C) Summation of infinite, alternating series

We write the program B to calculate the \(k\)-th term:

\(
\begin{align*}
a_k &= \frac{1}{3^k \, (2k + 1)}
\end{align*}
\)

Code:

084 - 42,21,12  LBL B
085 -        3  3
086 -       34  x<>y
087 -       14  y^x
088 -    43 36  LSTx
089 -        2  2
090 -       20  x
091 -        1  1
092 -       40  +
093 -       20  x
094 -       15  1/x
095 -    43 32  RTN

We use \(\text{PSum} = 10\), \(\text{NDif} = 7\) for maximum accuracy:

10 ENTER 7 A

After some time, we get the following result in the display:

0.9068996822

The correct result \(\frac{\pi}{\sqrt{12}}\) on the HP-11C is:

0.9068996821

---

Now you might be wondering: Isn't that like cracking a nut with a sledgehammer?

Let me use this program for the HP-42S instead:

Code:

00 { 17-Byte Prgm }
01▸LBL 00
02 1/X
03 LASTX
04 X<> ST Z
05 3
06 ÷
07 -
08 X<>Y
09 2
10 -
11 X>0?
12 GTO 00
13 R↓
14 END

We sum the terms backwards to improve accuracy.
It turns out that a starting value of \(43\) is good enough:

CLST
43
XEQ 00

We get:

0.906899682117

If we multiply that by \(\sqrt{12}\) we get:

3.14159265359

---

I must admit that I cheated a bit since I used this simulator for the HP-11C and Free42 instead of the real calculators.
Thus the results might be slightly off.

[robve]

(03-14-2022 01:55 PM)EdS2 Wrote:  There's an unusual one by Valentin... 3 years ago, which seems appropriate.

Here is the unusual pi computation in Haskell, which self-applies \( f\,(h,p)=(p\lceil\frac{h}{p}\rceil,p-1) \) iteratively \( n-1 \) times starting with \( (1,n) \) to return \( \frac{n^2}{h} \) as an approximation of \( \pi \):

unusual_pi n =
  let (h,p) = head (drop (n-2) (iterate (\(h,p) -> (p*ceiling (fromIntegral h/fromIntegral p), p-1)) (1,n)))
  in (fromIntegral (n*n))/(fromIntegral h)
unusual_pi 100000
3.141526183050601

The reason why this works is explained here: "Beginning with any positive integer \( n \), round up to the nearest multiple of \( n–1 \), then up to the nearest multiple of \( n–2 \), and so on, up to the nearest multiple of 1. Let \( f(n) \) denote the result. For example, \( f(10) = 34 \). Interestingly, the ratio \( n^2/f(n) \) approaches \( \pi \) (i.e., 3.14159...) as \( n \) increases."

(03-14-2022 01:55 PM)EdS2 Wrote:  It feels worthwhile also to showcase one of Gerson's [but see downthread!] pandigital approximations ("good to 17 digits"):
That is: (ln{[2×5!+(8-1)!]^(sqrt(9))+4!+(3!)!})/(sqrt(67))

Yes, very interesting and impressive work by Gerson.

- Rob

[Gerson W. Barbosa]

(03-14-2022 09:35 PM)robve Wrote:  

(03-14-2022 01:55 PM)EdS2 Wrote:  It feels worthwhile also to showcase one of Gerson's [but see downthread!] pandigital approximations ("good to 17 digits"):
That is: (ln{[2×5!+(8-1)!]^(sqrt(9))+4!+(3!)!})/(sqrt(67))

Yes, very interesting and impressive work by Gerson.

- Rob

Thank you, but that one is not totally mine. I’ve only put it in pandigital form.

The following is fully mine. Not so many digits, though:

\(\rm{e}\times\sqrt[12]{\rm{e}^{-3\times4}+5.67890}\)

[ttw]

One can often apply various sequence transformations to speed up these series. A couple of new methods are given here.

https://arxiv.org/pdf/1702.07199.pdf

https://people.mpim-bonn.mpg.de/zagier/f...lltext.pdf


In most record-setting attempts, the desired result is a fraction with really big integers, even bigger than 355/113. An (to me, anyway) amusing feature is that most of the time is taken by the final long division. Third order, first order, fourth order, etc., methods all take about the same amount of time as the time for arithmetic grows pretty fast.

[robve]

(03-14-2022 11:06 PM)ttw Wrote:  One can often apply various sequence transformations to speed up these series. A couple of new methods are given here.

https://arxiv.org/pdf/1702.07199.pdf

https://people.mpim-bonn.mpg.de/zagier/f...lltext.pdf


In most record-setting attempts, the desired result is a fraction with really big integers, even bigger than 355/113. An (to me, anyway) amusing feature is that most of the time is taken by the final long division. Third order, first order, fourth order, etc., methods all take about the same amount of time as the time for arithmetic grows pretty fast.

I recall there was also a HP Forum post on accelerating alternating sums with Albert Chan suggesting to use Aitken extrapolation correction. In addition, I observed that summing the terms in reverse order may sometimes improve accuracy.

Here is the code I wrote a year ago (for my collection):

' SUMALT - summation of infinite alternating series
' Euler transformation:
' https://en.wikipedia.org/wiki/Series_acceleration
' https://mathworld.wolfram.com/EulersSeri...ation.html
' Adapted from:
' https://albillo.hpcalc.org/programs/HP%2...Series.pdf
' Combined with Aitken extrapolation correction applied to series S2:
' https://en.wikipedia.org/wiki/Aitken%27s...ed_process
' Combined with reverse term summation
' Performs N+D+2 function evaluations and D*(D+1)/2+2*(D+1) table operations
' The term y(i) to sum is defined in line 200 using I as index and returning value Y

' VARIABLES
' N number of terms
' D order of differences, e.g. 7 should suffice
' S S1 and final sum_i=0^inf (-1)^i*y(i)
' A() auto array with difference terms t[0] to t[d] stored in A(27) to A(27+D)
' I sum iterator
' J iterator
' X S2
' Y value of y(i)
' Z scratch

' driver program
10 N=10,D=7
20 INPUT "N=";N,"D=";D
30 GOSUB 100: PRINT S: END

' init differences t[0] to t[d]
100 FOR I=N+1 TO N+D+1: GOSUB 200: A(26+I-N)=Y: NEXT I
' apply Euler transformation to compute differences t[0] to t[d] and store in A(27) to A(27+D)
110 FOR I=1 TO D: FOR J=D TO I STEP -1: A(27+J)=A(27+J)-A(26+J): NEXT J: NEXT I
120 Z=2
130 FOR I=0 TO D: A(27+I)=A(27+I)/Z,Z=-Z-Z: NEXT I
' apply Aitken extrapolation as a correction to initialize S2
140 S=0,X=-A(27+D)*A(27+D)/(A(27+D)-A(26+D)),Z=1
' compute S = S1 = sum_0^n (-1)^i*y(i)
150 FOR I=0 TO N: GOSUB 200: S=S+Z*Y,Z=-Z: NEXT I
' compute X = S2 = sum_0^d t[i] in reverse order to retain accuracy
160 FOR I=D TO 0 STEP -1: X=X+A(27+I): NEXT I
' S = S1+S2 if N is odd, S = S1-S2 otherwise, where Z=1 if N is odd or -1
170 S=S+Z*X
180 RETURN
' y(i)
200 Y=1/(3^I*(2*I+1)): RETURN

Line 200 is the alternating \( \displaystyle \frac{1}{3^i\,(2i+1)} \) term.

This method only takes N=5 to produce 8 decimal places of \( \pi \) with order of differences D=3:

RUN
N=5
D=3
         9.068996707E-01
9.068996707E-01*6/√3
3.141592614

Not bad!

- Rob

Edit note/warning: there is no reason to push this method hard with N=10 and D=7, which will require N+D+2=18 function evaluations and D*(D+1)/2+2(D+1)=44 tabulation and summation steps! That's far more CPU power than the N=17 steps for the Sharp series to converge to 10 decimal places.

Edit 2 The double precision version for the Sharp PC-1475 to compute up to 20 decimal places:

' driver program
10 DEFDBL: DEFDBL A,S,X-Z: N=10,D=7
20 INPUT "N=";N,"D=";D
30 ERASE A: DIM A(D)
40 GOSUB 100: PRINT S: END

' init differences t[0] to t[d]
100 FOR I=N+1 TO N+D+1: GOSUB 200: A(I-N-1)=Y: NEXT I
' apply Euler transformation to compute differences t[0] to t[d] and store in A(0) to A(D)
110 FOR I=1 TO D: FOR J=D TO I STEP -1: A(J)=A(J)-A(J-1): NEXT J: NEXT I
120 Z=2
130 FOR I=0 TO D: A(I)=A(I)/Z,Z=-Z-Z: NEXT I
' apply Aitken extrapolation as a correction to initialize S2
140 S=0,X=-A(D)*A(D)/(A(D)-A(D-1)),Z=1
' compute S = S1 = sum_0^n (-1)^i*y(i)
150 FOR I=0 TO N: GOSUB 200: S=S+Z*Y,Z=-Z: NEXT I
' compute X = S2 = sum_0^d t[i] in reverse order to retain accuracy
160 FOR I=D TO 0 STEP -1: X=X+A(I): NEXT I
' S = S1+S2 if N is odd, S = S1-S2 otherwise, where Z=1 if N is odd or -1
170 S=S+Z*X
180 RETURN

RUN
N=28
D=1
0.9068996821171089253
0.9068996821171089253#*6/√3
3.14159265358979793285

[Eddie W. Shore]

My favorite is the one where I memorized the digits: let's see...

3.1415926535923846...

HAPPY PI DAY!

I have been alive for 45 years today, 16,436 days.

[Ren]

How 3
I 1
wish 4
I 1
could 5
calculate 9
pi 2

[EdS2]

(03-14-2022 06:06 PM)Gerson W. Barbosa Wrote:  Thanks you, but it’s not mine.

Oops, noted.

Quote:Actually Ramanujan’s, I think.

Yikes, what a substitution!

I don't think I have a favourite series, but I am intrigued by the many and various "Machin-like formulas" each of which must somehow be discovered.

Katie's pi program, recently mentioned here, uses one of them, by Euler I think:

Code:

4 REM Uses the formula:  pi = 20*arctan(1/7) + 8*arctan(3/79)

and I'm wondering about converting it to use another and see what difference it makes. Gauss' looks like a good one.

[rprosperi]

(03-15-2022 01:09 AM)Eddie W. Shore Wrote:  My favorite is the one where I memorized the digits: let's see...

3.1415926535923846...

HAPPY PI DAY!

I have been alive for 45 years today, 16,436 days.

Happy Birthday Eddie. Another Pi baby, how appropriate!

[robve]

(03-14-2022 09:17 PM)Thomas Klemm Wrote:  That's an alternating series so we can use Valentin's program:
(11C) Summation of infinite, alternating series

We use \(\text{PSum} = 10\), \(\text{NDif} = 7\) for maximum accuracy:

10 ENTER 7 A

After some time, we get the following result in the display:

0.9068996822

The correct result \(\frac{\pi}{\sqrt{12}}\) on the HP-11C is:

0.9068996821

The number of function evaluations with n=10 and d=7 is n+d+2=19, but only 18 function evaluations are needed with Sharp's series to get 10 decimal places. So I am not sure if this is the best choice of parameterization for this particular alternating series.

To demonstrate what I mean, I did a bit of experimentation using a BASIC version of the algorithm with the PC-1475 double precision float (up to 20 digits, see my prior post). I found the following results to compute 19 decimal places:

Direct computation of Sharp's series to get 19 decimal places of pi:
35 function evaluations

Sum of alternating series with Aitken extrapolation and summing in reverse order to get 19 decimal places of pi:
n=26 d=1: 29 function evaluations + 5 table operations
n=25 d=3: 30 function evaluations + 14 table operations
n=24 d=4: 30 function evaluations + 20 table operations
n=23 d=5: 30 function evaluations + 27 table operations
n=22 d=6: 30 function evaluations + 35 table operations
n=21 d=7: 30 function evaluations + 44 table operations

The table operations are cheaper to compute compared to the function evaluations, but the overhead of increasing d grows quadratically. These computations are not negligible, so the smallest d should be chosen even if that requires increasing n by the same amount or slightly more.

The n=26 d=1 choice is optimal for 19 decimal places and slightly better than direct computation of the series. Testing with fewer places of pi show that d=1 is still the best choice.

This result is only applicable to the Sharp's alternating series. Other series will give different results and require different n and d parameterizations.

- Rob