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

[Albert Chan]

(06-23-2024 03:00 AM)Gerson W. Barbosa Wrote:  \(\pi \approx \pi{'} \left[1+\frac{1}{6}\left(\frac{\pi{'}}{2^{n+1} }\right)^2\right]\left[1+\frac{3}{40}\left(\frac{\pi{'}}{2^{n+1} }\right)^4\right]\left[1+\frac{9}{280}\left(\frac{\pi{'}}{2^{n+1} }\right)^6\right] \)

Correction factor had error of O(z^8), where z=pi'/2^(n+1)
If we use Horner's scheme, we might as well use correct z^8 coefficient.
It made code shorter (less terms), also more accurate, error = O(z^10)

XCas> expand( (1+1/6*z^2) * (1+3/40*z^4) * (1+9/280*z^6) )

1 + 1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 3/560*z^8 + 27/11200*z^10 + 9/22400*z^12

Replaced this with asin(z)/z, upto z^8 only

XCas> expand( 1+1/6*z^2 * (1+3²/(4*5)*z^2 * (1+5²/(6*7)*z^2 * (1+7²/(8*9)*z^2))) )

1 +1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 35/1152*z^8

To make code stack-only, we can't store n, and had to recover 2^(n+1) factor.
Code was using 2^(n+1) = IP((355/113)/z), but this will fail if n is too big.

Instead, I use log to recover n from z (Yes, this is stupid, all we need is to store n)

z = pi' / 2^(n+1)
log2(z) = log2(pi') - (n+1)
-(n+1) ≈ log2(z) - 1.6515
-(n+1) = IP(log2(z)) - 2

n=1 --> IP(log2(z)) - 2 = IP(-0.50) - 2 = -2
n=9 --> IP(log2(z)) - 2 = IP(-8.35) - 2 = -10
n=99 --> IP(log2(z)) - 2 = IP(-98.35) - 2 = -100
n=999 --> IP(log2(z)) - 2 = IP(-998.35) - 2 = -1000
...

Code:

00 { 79-Byte Prgm }
01▸LBL "Vh"
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          @ z = sin(pi/2^(n+1)) = pi'/2^(n+1)
14 X↑2
15 ENTER
16 ENTER
17 ENTER
18 49
19 ×
20 72
21 ÷
22 ×
23 +
24 25
25 ×
26 42
27 ÷
28 ×
29 +
30 .45
31 ×
32 ×
33 +
34 6
35 ÷            @ eps ≈ asin(z)/z - 1
36 X<>Y
37 SQRT
38 ENTER
39 LN           @ ln(z)   z   eps
40 2
41 LN
42 ÷            @ log2(z) z   eps   eps
43 IP
44 2
45 -
46 LASTX
47 X<>Y
48 Y↑X
49 ÷            @ pi'   eps   eps   eps
50 ×
51 LASTX
52 +
53 END

  1 XEQ "Vh" --> 3.138316846974900201597524772067293
10 XEQ "Vh" --> 3.141592653589793238462643383274432
11 XEQ "Vh" --> 3.141592653589793238462643383279499
12 XEQ "Vh" --> 3.141592653589793238462643383279504