Viète à grande Vitesse or Accelerated Viète (kind of)

[Gerson W. Barbosa]

.
In an earlier post I presented both an HP-42S/Free42 program and a correction factor to the Viète's formula for \(\pi\).

Starting with Viète's formula

\(\frac2\pi = \frac{\sqrt 2}2 \cdot \frac{\sqrt{2+\sqrt 2}}2 \cdot \frac{\sqrt{2+\sqrt{2+\sqrt 2}}}2 \cdots\)

we isolate \(\pi\), limit the expression to n terms and add a correction factor c₁ to it:

\(\pi \approx \frac {2^{n+1}}{ \sqrt 2 \cdot \sqrt{2+\sqrt 2} \cdot \sqrt{2+\sqrt{2+\sqrt 2}} }\cdot c_{1}\)

The correction factor c₁ is

\(c_{1} = 1 + \frac {\pi^{2}}{96\times 4^{n-1}}\)

Instead of \(\pi\) we use its approximate value obtained after the evaluation of the first n terms in the formula, which is enough to improve the original result.

For example, when n = 3

\(\pi \approx \frac {2^{4}}{ \sqrt 2 \cdot \sqrt{2+\sqrt 2} \cdot \sqrt{2+\sqrt{2+\sqrt 2}} }\)

\(\pi \approx 3.12144515227\)

At this point we can use the approximate value of \(\pi\) thus obtained to improve it a bit more:

\(\pi \approx 3.12144515227 \cdot c_{1}\)

\(\pi \approx 3.12144515227 \cdot \left(1 + \frac {3.12144515227^{2}}{96\times 4^{2}}\right)\)

\(\pi \approx 3.14124564095 \)

Or, on the HP-42S:

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

3 XEQ "VWπ" ->

Y: 3.12144515227
X: 3.14124564095

As we can see, despite the initial poor \(\pi\) approximation this correction factor has provided two extra correct digits compared to the ones obtained with the plain Viète's formula. That's nice, but it can be better. Why not find a second correction factor? Indeed, again the errors form a geometric progression, except that now the ratio is 16 rather than 4. The second correction factor is easy to find, we only need to compute the first element of the GP. This leads to

\(c_{2} = 1 + \frac {3\cdot\pi^{4}}{10240\times 16^{n-1}}\)

Now we have

\(\pi \approx \left(\frac {2^{n+1}}{ \sqrt 2 \cdot \sqrt{2+\sqrt 2} \cdot \sqrt{2+\sqrt{2+\sqrt 2}} }\cdot c_{1}\right)c_{2}\)

Let's check this out on Free42 and on the HP-42S:

Code:


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

20 XEQ "VV" ->

Z: 3.141592653588618236614208590772408
Y: 3.141592653589793238462642196715985
X: 3.141592653589793238462643383279502


On the HP-42S:

1 XEQ "VV" ->

Z: 2.82842712475
Y: 3.06412938514
X: 3.14326223323

2 XEQ "VV" ->

Z: 3.06146745894
Y: 3.13619104715
X: 3.14174642842

3 XEQ "VV" ->

Z: 3.12144515227
Y: 3.14124564095
Z: 3.14159565930

4 XEQ "VV" ->

Z: 3.13654849054
Y: 3.14157081551
X: 3.14159270299

5 XEQ "VV" ->

Z: 3.14033115695
Y: 3.14159128636
X: 3.14159265437

9 XEQ "VV" ->

Z: 3.14158772526
Y: 3.14159265356
X: 3.14159265359

I haven't investigated this any further, but it looks like we have something here. Most likely a third correction factor can be determined, and a fourth, and a fifth... If they follow a pattern than this might make for a viable method for calculation \(\pi\) to many digits. Take, for instance, the integer constants that appear in the first two correction factors, 96 and 10240. They can be written as 96 = 8² + 2⁵ and 96² + 2¹⁰, respectively. Task for another day...