Stoneham's series

[ttw]

Stoneham has a few funny series. Here are two that converge to Sqrt(2) and Pi. These seem similar in pattern but one has a quadratic surd for a limit and the other a transcendental number.

Sqrt(2) = 2 * Prod (1-1/(4*d^2) taken over all d odd.

Pi= 4 * Prod( 1-1/d^2) taken over all d>1.

"Random Generators and Normal Numbers"
David H. Bailey and Richard E. Crandall
from
Experimental Mathematics, Vol. 11 (2002), No. 4

[KeithB]

I assume that prod() is the infinite product.
I don't see how the second converges to pi.
The first term is .75, and 4 * .75 is 3, and all the other terms get progressively closer to 1.

[Thomas Klemm]

(06-04-2024 02:20 PM)KeithB Wrote:  I don't see how the second converges to pi.

He forgot to mention "over all odd d > 1":

Code:

from math import prod

2 * prod(1 - 1/4/d**2 for d in range(1, 1_000_000, 2))

1.4142137391498133

Code:

4 * prod(1 - 1/d**2 for d in range(3, 1_000_000, 2))

3.1415942243865067

[Nigel (UK)]

(06-04-2024 02:20 PM)KeithB Wrote:  I assume that prod() is the infinite product.
I don't see how the second converges to pi.
The first term is .75, and 4 * .75 is 3, and all the other terms get progressively closer to 1.

You are correct; but looking at the original reference (page 2 of the pdf), the product should be over all odd d greater than one.

Nigel (UK)

[Thomas Klemm]

In both cases we can use the infinite product representation of the sinc function:

\(
\begin{align}
\frac {\sin \pi z}{\pi z}=\prod _{n=1}^{\infty }\left(1-{\frac {z^{2}}{n^{2}}}\right)
\end{align}
\)

For \(z=\frac{1}{2}\) we have:

\(
\begin{align}
\frac {\sin \pi z}{\pi z}
&=\frac {\sin \pi \frac{1}{2}}{\pi \frac{1}{2}} \\
\\
&=\frac{2}{\pi} \\
\\
&=\prod _{n=1}^{\infty }\left(1-{\frac {1}{(2n)^2}}\right) \\
\\
&=\prod _{n=1}^{\infty }\frac{(2n-1)(2n+1)}{(2n)(2n)} \\
\\
&=\frac{1 \cdot 3}{2 \cdot 2} \cdot \frac{3 \cdot 5}{4 \cdot 4} \cdot \frac{5 \cdot 7}{6 \cdot 6} \cdots \\
\end{align}
\)

And thus:

\(
\begin{align}
\frac{\pi}{2}&=\frac{2 \cdot 2}{1 \cdot 3} \cdot \frac{4 \cdot 4}{3 \cdot 5} \cdot \frac{6 \cdot 6}{5 \cdot 7} \cdots \\
\\
\pi
&=4 \cdot \frac{2 \cdot 4}{3 \cdot 3} \cdot \frac{4 \cdot 6}{5 \cdot 5} \cdot \frac{6 \cdot 8}{7 \cdot 7} \cdots \\
\\
&=4 \cdot \left(1 - \frac{1}{3^2}\right) \cdot \left(1 - \frac{1}{5^2}\right) \cdot \left(1 - \frac{1}{7^2}\right) \cdots \\
\end{align}
\)

Similarly for \(z=\frac{1}{4}\) we have:

\(
\begin{align}
\frac {\sin \pi z}{\pi z}
&=\frac {\sin \pi \frac{1}{4}}{\pi \frac{1}{4}} \\
\\
&=\frac{2 \sqrt{2}}{\pi} \\
\\
&=\prod _{n=1}^{\infty }\left(1-{\frac {1}{(4n)^2}}\right) \\
\\
&=\prod _{n=1}^{\infty }\frac{(4n-1)(4n+1)}{(4n)(4n)} \\
\\
&=\frac{3 \cdot 5}{4 \cdot 4} \cdot \frac{7 \cdot 9}{8 \cdot 8} \cdot \frac{11 \cdot 13}{12 \cdot 12} \cdots \\
\end{align}
\)

If we multiply both products we get:

\(
\begin{align}
\frac{\pi}{2} \cdot \frac{2 \sqrt{2}}{\pi}
&= \sqrt{2} \\
\\
&=\frac{2 \cdot 2}{1 \cdot 3} \cdot \frac{6 \cdot 6}{5 \cdot 7} \cdot \frac{10 \cdot 10}{9 \cdot 11} \cdots \\
\\
\sqrt{2} = \frac{2}{\sqrt{2}}
&= 2 \cdot \frac{1 \cdot 3}{2 \cdot 2} \cdot \frac{5 \cdot 7}{6 \cdot 6} \cdot \frac{9 \cdot 11}{10 \cdot 10} \cdots \\
\\
&=2 \cdot \left(1 - \frac{1}{4 \cdot 1^2}\right) \cdot \left(1 - \frac{1}{4 \cdot 3^2}\right) \cdot \left(1 - \frac{1}{4 \cdot 5^2}\right) \cdots \\
\end{align}
\)

---

This might remind you of the Wallis product for \(\pi\):
Proof using Euler's infinite product for the sine function

[Albert Chan]

(06-04-2024 07:57 PM)Thomas Klemm Wrote:  In both cases we can use the infinite product representation of the sinc function:

\(
\begin{align}
\frac {\sin \pi z}{\pi z}=\prod _{n=1}^{\infty }\left(1-{\frac {z^{2}}{n^{2}}}\right)
\end{align}
\)

Another way, identities derived from cos(x) infinite product.

Si(pi*z)    = (1-(z/1)²) * (1-(z/2)²) * (1-(z/3)²) * (1-(z/4)²) ...
Si(pi*z/2) = (1-(z/2)²) * (1-(z/4)²) * (1-(z/6)²) * (1-(z/8)²) ...

Si(pi*z/2) * cos(pi*z/2) = Si(pi*z)
cos(pi*z/2) = (1-(z/1)²) * (1-(z/3)²) * (1-(z/5)²) * (1-(z/7)²) ...

Set z=1/2, we have:

√2/2 = (1-¼/1²) * (1-¼/3²) * (1-¼/5²) * (1-¼/7²) ...

Set z=1, we have 0 = 0 * ?. We move RHS first term to the left, to find out what ? is.
(this is why pi/4 infinite product skipped over first odd number)

cos(pi*z/2) / (1-z²) = (1-(z/3)²) * (1-(z/5)²) * (1-(z/7)²) * (1-(z/9)²) ...

limit(LHS, z=1) = limit(-sin(pi*z/2)*(pi/2) / (-2z), z=1) = pi/4

pi/4 = (1-1/3²) * (1-1/5²) * (1-1/7²) * (1-1/9²) ...

[Thomas Klemm]

Grant Sanderson has a nice explanation:

[Johnh]

I got interested in the infinite product for Pi, presented above:

Pi= 4 * Prod( 1-1/d^2) taken over all d>1, for d = odd numbers.

So I made a small routine for an Hp15c:

Code:


 
   001 { 42 21 11 } f LBL A
   002 {        2 } 2
   003 { 44 40  1 } STO + 1
   004 {    45  1 } RCL 1
   005 {       15 } 1/x
   006 {    43 11 } g x^2
   007 {       16 } CHS
   008 {        1 } 1
   009 {       40 } +
   010 { 44 20  2 } STO * 2
   011 {    22 11 } GTO A

Put 1 into Memory 1 and 4 into Memory 2, and then run A

Memory 1 starts at 1 and increments by 2 at each cycle, and Memory 2 starts with 4 and gets shaved gradually down towards Pi by multiplication with factors just slightly less than 1.

It's not a very fast method to get close to Pi, so needs many loops. The routine as written above runs without pausing or stopping to display. But adding RCL 2 and PSE after line 10 will let the steps be observed.

So I tried this on my Hp15C-CE and on two Android emulators on my Samsung S22, each set to run at max speed. Each ran for 60 seconds.

The real Hp15C-CE, incremented d from 1 up to 12135 in steps of 2, and achieved a PI approximation of 3.14172210.

Touch RPN, running in Hp15 mode incremented to 62631 and estimated Pi as 3.141617733

Jovial JRPN15 incremented to 2,664,645 and got to a Pi estimate of 3.141607465.

Jovial wins the test, much faster, doing over 1.3 million loops in a minute!, though not much extra convergence is achieved with this formula, for all the extra cycles.

Also, Jovial JRPN15 provided my code listing above as a direct copy to the Clipboard, which is a very handy feature, plus a copy-paste of its two number results. Touch RPN can also copy result numbers, and saves memory contents in a text format which has key codes, though not the key names in the current version and its not so easily accessible. But I prefer it for it's look and feel.

Maybe more interesting than useful, but a good test of these three Hp15 versions.

[Albert Chan]

(06-07-2024 07:55 PM)Johnh Wrote:  Jovial JRPN15 incremented to 2,664,645 and got to a Pi estimate of 3.141607465.

Jovial wins the test, much faster, doing over 1.3 million loops in a minute!, though not much extra convergence is achieved with this formula for all the extra cycles.

The problem is (1 - 1/b^2), when b goes big ... it is rounded to 1.
Doing product directly over-estimated what the algorithm should offer.

(1 - a) * (1 - b) = 1 - a - b + a*b = 1 - (a + (1-a)*b)

Code:

function p(n) -- 1 - product((1-1/b^2), b = 3 .. n step 2)
    local a = 0
    for b=3,n,2 do a = a + (1-a)/(b*b) end
    return a
end

lua> p(2664645)
0.21460168922870704
lua> (1-_) * 4
3.141593243085172

Or, we do product backwards, for more accuracy.

Code:

function q(n) -- 1 - product((1-1/b^2), b = n-(n+1)%2 .. 3 step -2)
    local a = 0
    for b=n-(n+1)%2,3,-2 do a = a + (1-a)/(b*b) end
    return a
end

lua> q(2664645)
0.2146016892287097
lua> (1-_) * 4
3.141593243085161


For reference, this is without rounding error (40 digits shown)

3.141593243085161166771739970151562977338...

[Johnh]

Thanks Albert, that's clever and insightful algebra wrangling, as are many of your posts!

I can see that the basic formula can never keep accuracy past about four decimal places, without an unusual number of calculation digits compared to many calculators. I might try it on Free42 with its 34 digits, just for fun, and to learn more about Free42.

But it works well enough to explore though, and I hadn't seen a method for Pi before like that, being based fully on products.

[Johnh]

....further on this, I put it into Free42, similar code that I posted above, with a 'View' line added to monitor the convergence towards Pi. I let it run 30 minutes, by which time it had completed 84 million loops, which is more than twice as many per minute as the Jovial hp15 sim!

The first 4 or 5 decimal places are almost instant, and when I stopped it it was at 3.14159266291 and counting down maybe one digit every few seconds in the 11th decimal place, slowing down exponentially, with an error of one digit remaining in the 8th decimal place. Maybe itd get that one given long enough, or maybe the accumulation of rounding error would throw it off.

I was very impressed with the Free42 App, my first time coding in it. Nice and clear and simple but powerful with handy alpha-numerics, and really great performance.

All more interesting than what I'm supposed to be doing this Saturday afternoon (marking engineering exams...yuk!)

[Albert Chan]

(06-08-2024 05:02 AM)Johnh Wrote:  I put it into Free42, similar code that I posted above, with a 'View' line added to monitor the convergence towards Pi. I let it run 30 minutes, by which time it had completed 84 million loops ... stopped it it was at 3.14159266291 and counting down maybe one digit every few seconds in the 11th decimal place, slowing down exponentially

Convergence is slow, similar to Basel problem, error ≈ 1/n
see thread Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B]

(06-07-2024 08:54 PM)Albert Chan Wrote:  Or, we do product backwards, for more accuracy.

Code:

function q(n) -- 1 - product((1-1/b^2), b = n-(n+1)%2 .. 3 step -2)
    local a = 0
    for b=n-(n+1)%2,3,-2 do a = a + (1-a)/(b*b) end
    return a
end

With b going big, (1-a) does not change much.
If we assume (1-a) = 1, product turned into summation.

a = 1/3^2 + 1/5^2 + 1/7^2 + ... + 1/n^2 = (pi^2/8 - 1) if n→∞

Because this only sum odd square reciprocal, error ≈ .5/n

We are using (1-a)*4 for pi estimate, so we expected error ≈ -2/n
Of course (1-a) ≈ pi/4 ≠1 --> my guess is error around (-pi/2)/n

lua> n = 84e6*2+1
lua> (1-q(n)) * 4
3.1415926629397712
lua> _ + (-pi/2)/n
3.141592653589793

[Gerson W. Barbosa]

(06-08-2024 05:02 AM)Johnh Wrote:  ....further on this, I put it into Free42, similar code that I posted above, with a 'View' line added to monitor the convergence towards Pi. I let it run 30 minutes, by which time it had completed 84 million loops, which is more than twice as many per minute as the Jovial hp15 sim!

Here’s an equivalent expression plus a correction factor which allows for 33 correct digits with only fifteen loops:

\[\pi \approx \left(4\prod _{k=1}^{n}{\frac {(2k+1)^2-1}{(2k+1)^2}}\right) \left({1-\frac{2}{8n+9+\frac{1\times3}{8n+8+\frac{3\times5}{8n+8+… +\frac{\left({2n}\right)^2-1}{8n+8}}}}}\right)\]

Code:


00 { 29-Byte Prgm }
01▸LBL "Sπ"
02 4
03▸LBL 00
04 RCL ST Y
05 STO+ ST X
06 SIGN
07 RCL+ ST L
08 X↑2
09 ABS
10 DSE ST X
11 RCL÷ ST L
12 ×
13 DSE ST Y
14 GTO 00
15 END

1332322 XEQ Sπ ->

3.141593243085161166771739970151(294) (~ 3 seconds on Free42)

Code:


00 { 68-Byte Prgm }
01▸LBL "SWπ"
02 4
03 STO 01
04 STO+ ST X
05 RCL× ST Y
06 RCL+ ST L
07 RCL ST X
08▸LBL 01
09 RCL ST Z
10 STO+ ST X
11 X↑2
12 DSE ST X
13 X<>Y
14 ÷
15 RCL+ ST Y
16 RCL ST Z
17 STO+ ST X
18 SIGN
19 RCL+ ST L
20 X↑2
21 ABS
22 DSE ST X
23 RCL÷ ST L
24 STO× 01
25 R↓
26 DSE ST Z
27 GTO 01
28 1
29 +
30 1/X
31 STO+ ST X
32 +/-
33 1
34 +
35 RCL× 01
36 END

15 XEQ SWπ ->

3.14159265358979323846264338327950(4)

5 XEQ SWπ ->

3.1415926535(8) (~ 2 seconds on the HP-42S)

[Namir]

Very cool using only stack registers!!!

Namir

[Gerson W. Barbosa]

(06-09-2024 10:55 PM)Namir Wrote:  Very cool using only stack registers!!!

Namir

Except for the second program, but I’ll leave that as an exercise to the reader :-)

[Gerson W. Barbosa]

(06-09-2024 11:04 PM)Gerson W. Barbosa Wrote:  Except for the second program, but I’ll leave that as an exercise to the reader :-)

No takers?

Code:


00 { 72-Byte Prgm }
01▸LBL "SWπ"
02 8
03 0
04 8
05 RCL× ST T
06 RCL+ ST Z
07▸LBL 00
08 R↑
09 NOT
10 +/-
11 STO+ ST X
12 DSE ST X
13 STO× ST X
14 STO÷ ST T
15 DSE ST X
16 STO× ST T
17 X<> ST L
18 RCL+ ST X
19 STO× ST X
20 DSE ST X
21 X<>Y
22 STO+ ST Z
23 X<> ST Z
24 STO÷ ST Y
25 X<> ST L
26 R↓
27 X<>Y
28 DSE ST T
29 GTO 00
30 +
31 1
32 +
33 1/X
34 +/-
35 0.5
36 +
37 ×
38 END

15 XEQ SWπ ->

3.14159265358979323846264338327950(1)

This lacks some optimization, though.

[Albert Chan]

(06-09-2024 10:36 PM)Gerson W. Barbosa Wrote:  \[\pi \approx \left(4\prod _{k=1}^{n}{\frac {(2k+1)^2-1}{(2k+1)^2}}\right) \left({1-\frac{2}{8n+9+\frac{1\times3}{8n+8+\frac{3\times5}{8n+8+… +\frac{\left({2n}\right)^2-1}{8n+8}}}}}\right)\]

Nice!

I have trouble doing correction for √2 ... it just seems too complicated.

\(\sqrt{2} \approx \displaystyle
\left(2\prod _{k=1}^{n} \left( 1 - \frac {1/4}{(2k-1)^2} \right) \right)
\left({1-\frac{2}
{1+32n+\frac{(4^2-1)^2/(2^2-1)}
{32n+\frac{(8^2-1)^2/(4^2-1)}
{32n+… +\frac{{((4n)^2-1)^2/((2n)^2-1)}}{32n}}}}}
\right)\)

[Gerson W. Barbosa]

(06-11-2024 05:12 PM)Albert Chan Wrote:  I have trouble doing correction for √2 ... it just seems too complicated.

\(\sqrt{2} \approx \displaystyle
\left(2\prod _{k=1}^{n} \left( 1 - \frac {1/4}{(2k-1)^2} \right) \right)
\left({1-\frac{2}
{1+32n+\frac{(4^2-1)^2/(2^2-1)}
{32n+\frac{(8^2-1)^2/(4^2-1)}
{32n+… +\frac{{((4n)^2-1)^2/((2n)^2-1)}}{32n}}}}}
\right)\)

Very nice! I presume finding the pattern of the numerators may have been particular troublesome, as they are not trivial. Your correction factor gives about two or perhaps 25/12 correct decimal digits per iteration:

. n |2*Product(((4k - 2)^2 - 1)/(4k - 2)^2,{k,1,n})*(1 + -2/(32n + 1 + ContinuedFractionK[((4k)^2-1)^2/((2k)^2-1),32n,{k,1,n}]))
.---+---------------------------------------------------------------------------------------------------------------------------
. 1 | 1.4151193633952254641910 | 1067/754
. 2 | 1.4142087816013794108361 | 465863/329416
. 3 | 1.4142135946217028898547 | 43179041/30532192
. 4 | 1.4142135621369493256490 | 99253807/70183040
. 5 | 1.4142135623748906144019 | 664950544227/470191038976
. 6 | 1.4142135623730811143616 | 153235360737661/108353762695168
. 7 | 1.4142135623730951582795 | 3602252403912943/2547177102352384
. 8 | 1.4142135623730950479347 | 620602650775941779/438832342786015232
. 9 | 1.4142135623730950488086 | 15082567717216397339/10664985910549020672
.10 | 1.4142135623730950488016 | 15606594740681145876391/11035528972365946421248


Of course, this is not adequate for computing the decimal digits of π or any other mathematical constant for that matter. Anyway, just as another mathematical curiosity, I thought of combining the Wallis Product and the Stoneham's product for √2/2:

π/2=4/3.16/15.36/35.64/63.100/99.144/143. …

√2/2=3/4.35/36.99/100. …

π/2.√2/2=4/3.3/4.16/15.36/35.35/36.64/63.100/99.99/100.144/143. …

π√2/4=16/15.64/33.144/143. …

π=2.√2.Π[k=1,inf,(4k)^2/((4k)^2-1)]

and

π~2.√2.Π[k=1,n,(4k)^2/((4k)^2-1)]*c,

where 'c' is a correction factor.

Thanks for saving me the trouble to find a correction factor for this one! it turns out that yours for the √2 product almost exactly matches it:

. n |2√2*Product((4k)^2/((4k)^2 - 1),{k,1,n})*(1 + 2/(32n + 15 + ContinuedFractionK[((4k)^2-1)^2/((2k)^2-1),32n + 16,{k,1,n}]))
.---+--------------------------------------------------------------------------------------------------------------------------
. 1 | 3.1412407295336494743353 | 25888/11655*√2
. 2 | 3.1415952167705760328966 | 1863108608/838692855*√2
. 3 | 3.1415926338285612147454 | 2607005335552/1173564727335*√2
. 4 | 3.1415926537455103768956 | 86592236552192/38980201708305*√2
. 5 | 3.1415926535885530258848 | 42678860747884199936/19212237343166192325*√2
. 6 | 3.1415926535898031757455 | 1474308092195358988304384/663671815222966253879175*√2
. 7 | 3.1415926535897931585427 | 25754557013712857799386464256/11593623947422059715601253825*√2
. 8 | 3.1415926535897932391070 | 16248219196766078634894964031488/7314268425672804441304956579225*√2
. 9 | 3.1415926535897932384574 | 3482380024930917685304486445187072/1567621777752447501676104398875275*√2
.10 | 3.1415926535897932384627 | 897305682114824147754072834090261807104/403929473094228867300191362519982221425*√2

[Gerson W. Barbosa]

(06-12-2024 07:01 PM)Gerson W. Barbosa Wrote:  . n |2√2*Product((4k)^2/((4k)^2 - 1),{k,1,n})*(1 + 2/(32n + 15 + ContinuedFractionK[((4k)^2-1)^2/((2k)^2-1),32n + 16,{k,1,n}]))
.---+--------------------------------------------------------------------------------------------------------------------------
. 1 | 3.1412407295336494743353 | 25888/11655*√2
. 2 | 3.1415952167705760328966 | 1863108608/838692855*√2

This result is also a consequence of a connection of a series of Wallis-like products to the Viète's formula for \(\pi\), which I had never noticed before:

\(
\begin{align}
\space\space\space\space \prod _{n=1}^{\infty}{\frac {(2^kn)^2}{(2^kn)^2-1}}
\\
\\ &=\frac{\pi}{2} \space\space\space\space,\space k = 1 \\
\\ &=\frac{\pi}{2\sqrt{2}} \space\space\space\space,\space k = 2\\
\\ &=\frac{\pi}{4\sqrt{2-\sqrt{2}}} \space\space\space\space,\space k = 3\\
\\ &=\frac{\pi}{8\sqrt{2-\sqrt{2+\sqrt{2}}}} \space\space\space\space,\space k = 4\\
\\ &=\frac{\pi}{16\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2}}}}} \space\space\space\space,\space k = 5\\
\\ &\cdots
\end{align}
\)

This HP-42S/Free42 program computes Viète's results for k=3 and above. The program uses the original version of Viètes's formula in order to avoid cancellation errors. The steps 22 through 45 are optional. That's just an experimental method for obtaining a few extra digits of \(\pi\), kept in the stack register Y at program exit.

Code:


00 { 82-Byte Prgm }
01▸LBL "VWπ"
02 2
03 -
04 LASTX
05 STO 01
06 STO+ 01
07 RCL ST X
08 SQRT
09 RCL ST X
10▸LBL 00
11 RCL+ ST Z
12 SQRT
13 STO× ST Y
14 X<> 01
15 STO+ ST X
16 X<> 01
17 DSE ST T
18 GTO 00
19 R↓
20 RCL÷ 01
21 1/X
22 RCL 01
23 X↑2
24 16
25 ÷
26 RCL ST X
27 64
28 ÷
29 STO- ST Y
30 8
31 ÷
32 STO+ ST Y
33 2
34 ÷
35 STO- ST Y
36 4
37 ÷
38 STO+ ST Y
39 4
40 ÷
41 -
42 RCL× ST Y
43 1/X
44 RCL+ ST Y
45 X<>Y
46 END

44 XEQ "VWπ" ->

Y: 3.141592653589793238462643383279502
X: 3.141592653589793238462643366581728

—————

P. S.: This uses a better correction term for Viète's formula and works also for k=2:

Code:


00 { 79-Byte Prgm }
01▸LBL "VWπ"
02 DSE ST X
03 STO 01
04 0.5
05 1
06 0
07▸LBL 00
08 RCL× ST Z
09 RCL+ ST Z
10 SQRT
11 STO× ST Y
12 DSE ST T
13 GTO 00
14 2
15 RCL÷ ST Z
16 16
17 10
18 1/X
19 2
20 +
21 -8
22 Y↑X
23 3
24 -
25 10
26 ÷
27 X↑2
28 X↑2
29 -
30 15
31 X<>Y
32 ÷
33 2
34 RCL 01
35 Y↑X
36 SQRT
37 ×
38 X↑2
39 X↑2
40 1/X
41 RCL+ ST Y
42 X<>Y
43 END

44 XEQ "VWπ" ->

Y: 3.141592653589793238462643383279504
X: 3.141592653589793238462643366581725

On the HP-42S:

10 XEQ "VWπ" ->

Y: 3.14159265358
X: 3.14158772527

—————

P. P. S.: Stack-only, k≥2:

Code:



00 { 70-Byte Prgm }
01▸LBL "VWπ"
02 DSE ST X
03 2
04 0
05 N!
06▸LBL 00
07 2
08 STO× ST Z
09 RCL+ ST L
10 SQRT
11 ×
12 DSE ST Z
13 GTO 00
14 RCL÷ ST Y
15 1/X
16 X<>Y
17 112.5
18 ×
19 X↑2
20 2.1
21 -8
22 Y↑X
23 3
24 -
25 10
26 ÷
27 X↑2
28 X↑2
29 16
30 -
31 X↑2
32 X↑2
33 X<>Y
34 ÷
35 RCL+ ST Y
36 X<>Y
37 END

43 XEQ "VWπ" ->

Y: 3.141592653589793238462643383279500
X: 3.141592653589793238462643316488386


On the HP-42S:

2 XEQ "VWπ" ->

Y: 3.15140913388
X: 2.82842712476

3 XEQ "VWπ" ->

Y: 3.14221296121
X: 3.06146745893

10 XEQ "VWπ" ->

Y: 3.14159265358
X: 3.14158772527

[Gerson W. Barbosa]

(06-14-2024 05:42 PM)Gerson W. Barbosa Wrote:  P. P. S.: Stack-only, k≥2:

Code:



00 { 70-Byte Prgm }
01▸LBL "VWπ"
02 DSE ST X
03 2
04 0
05 N!
06▸LBL 00
07 2
08 STO× ST Z
09 RCL+ ST L
10 SQRT
11 ×
12 DSE ST Z
13 GTO 00
14 RCL÷ ST Y
15 1/X
16 X<>Y
17 112.5
18 ×
19 X↑2
20 2.1
21 -8
22 Y↑X
23 3
24 -
25 10
26 ÷
27 X↑2
28 X↑2
29 16
30 -
31 X↑2
32 X↑2
33 X<>Y
34 ÷
35 RCL+ ST Y
36 X<>Y
37 END

43 XEQ "VWπ" ->

Y: 3.141592653589793238462643383279500
X: 3.141592653589793238462643316488386


On the HP-42S:

2 XEQ "VWπ" ->

Y: 3.15140913388
X: 2.82842712476

3 XEQ "VWπ" ->

Y: 3.14221296121
X: 3.06146745893

10 XEQ "VWπ" ->

Y: 3.14159265358
X: 3.14158772527

I've finally found the correction factor for the Vieta's formula for \(\pi\) I was looking for, which turns the previous programs obsolete. The correction factor is


\(c = 1 + \frac {1}{\frac{96}{\pi^2 }4^\left(n-1\right)}\)

This simplifies to

\(c = 1 + \frac {\pi^2\times 2^\left(-2n-3\right)}{3}\)

where n is the number of iterations and \(\pi\) is the approximated value obtained at loop exit

The program uses a more compact version of the latter for size-optimization.

Code:


00 { 39-Byte Prgm }
01▸LBL "VWπ"
02 2
03 0
04 SIGN
05▸LBL 00
06 2
07 STO× ST Z
08 RCL+ ST L
09 SQRT
10 ×
11 DSE ST Z
12 GTO 00
13 1/X
14 STO× ST Y
15 X↑2
16 6
17 +
18 RCL÷ ST L
19 RCL× ST Y
20 END

32 XEQ "VWπ" ->

Y: 3.141592653589793238462607815493113
X: 3.141592653589793238462643383279503

Notice that from 22 correct decimal digits the correction factor gets another eleven, a 50 % increase in this case.


On the HP-42S:

1 XEQ "VWπ" ->

Y: 2.82842712475
X: 3.06412938514

2 XEQ "VWπ" ->

Y: 3.06146745894
X: 3.13619104715

3 XEQ "VWπ" ->

Y: 3.12144515227
X: 3.14124564095

4 XEQ "VWπ" ->

Y: 3.13654849054
X: 3.14157081551

5 XEQ "VWπ" ->

Y: 3.14033115695
X: 3.14159128636

9 XEQ "VWπ" ->

Y: 3.14158772526
X: 3.14159265356