Hi,
Once again, I had another discussion with my cousin about PI (it's either PI or prime numbers...).
He told me about the beauty of the Viète's formula :
[
attachment=8585]
History here:
link
This video shows the formula entered in CAS, and yes! I find it beautiful:
https://youtu.be/BEYDBl94UpM
Regards,
Thibault
Cool! Really beautiful to see in the CAS the pi approximation in action! :-)
very nice, thank you!
tiny tip: in CAS, for algebraic or fractional input, the a.b/c key also runs approx(); one less key-stroke than shift-enter.
(06-18-2020 12:51 PM)pinkman Wrote: [ -> ]Compared to Wallis formula (https://www.hpmuseum.org/forum/thread-14...ght=wallis), it's absolutely fast!
Yes, but it pales in comparison to the
Wallis-Wasicki formula :-)
+---+---------------------+---------------------+
| 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 |
+---+---------------------+---------------------+
Only 7 iterations (or 14, depending on how you implement the algorithm) for 18 correct decimal digits
(21697209162666264236130304/6906436179074198667175275 = 3.141592653589793238[633]...).
The Pascal code is available
here. It should translate easily into the Prime programming language, but probably no joy with only 12 significant digits...
(06-18-2020 12:51 PM)pinkman Wrote: [ -> ]And the convergence speed is acceptable : for p(18) we have 10 digits.
p(18) means Σ(x,x,1,18) = 171 iterations.
For faster convergence, you might want to try this algorithm:
Code:
Program Viete;
Uses Crt;
var k: Byte;
n: LongInt;
c: Char;
d, p, t, v: Extended;
begin
ClrScr;
WriteLn('+---+----------------------+----------------------+');
WriteLn('| k | 2*v | ~ pi |');
WriteLn('+---+----------------------+----------------------+');
v := 1;
n := 1;
d := 0;
for k := 1 to 15 do
begin
n := n + n;
d := Sqrt(d + 2);
v := v*d;
t := 6*v*v - 2/3;
p := 2*n/v*(t + 1)/t;
WriteLn('|',k:2,' |',2*n/v:21:18,' |',p:21:18,' |')
end;
WriteLn('+---+----------------------+----------------------+');
c := ReadKey
end.
+---+----------------------+----------------------+
| k | 2*v | ~ pi |
+---+----------------------+----------------------+
| 1 | 2.828427124746190098 | 3.077994223988500989 |
| 2 | 3.061467458920718174 | 3.137427049842311165 |
| 3 | 3.121445152258052286 | 3.141329731135014592 |
| 4 | 3.136548490545939264 | 3.141576182507746574 |
| 5 | 3.140331156954752913 | 3.141591623553754287 |
| 6 | 3.141277250932772868 | 3.141592589203296386 |
| 7 | 3.141513801144301077 | 3.141592649565492850 |
| 8 | 3.141572940367091385 | 3.141592653338272210 |
| 9 | 3.141587725277159701 | 3.141592653574073140 |
|10 | 3.141591421511199975 | 3.141592653588810732 |
|11 | 3.141592345570117743 | 3.141592653589731833 |
|12 | 3.141592576584872667 | 3.141592653589789401 |
|13 | 3.141592634338562990 | 3.141592653589793000 |
|14 | 3.141592648776985670 | 3.141592653589793224 |
|15 | 3.141592652386591346 | 3.141592653589793238 |
+---+----------------------+----------------------+
Quote:Yes, but it pales in comparison to the Wallis-Wasicki formula :-)
Here is a quick PPL port of your Wallis-Wasicki implementation:
Code:
EXPORT WallisWasicki()
BEGIN
LOCAL i, k, m, n;
LOCAL c, d, t, w;
PRINT("+---+-------------+------------+");
PRINT("| N | 2*W | 2*W*C |");
PRINT("+---+-------------+------------+");
FOR k := 1 TO 7 DO
n := 2*k;
w := 1;
c := 0;
d := 4*n;
m := n + 1;
FOR i := 1 TO k DO
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;
PRINT("|" + n + " |" + 2*w + " |" + 2*w*c + " |");
END;
PRINT("+---+-------------+------------+");
RETURN 2*w*c;
END;
The terminal output is not well formatted, but still easy to read.
I’m also porting your 314 pi digits algorithm, I had to stop (work...) but it will be ready in a few hours (if I find the time).
(06-23-2020 06:39 PM)Gerson W. Barbosa Wrote: [ -> ] (06-18-2020 12:51 PM)pinkman Wrote: [ -> ]And the convergence speed is acceptable : for p(18) we have 10 digits.
p(18) means Σ(x,x,1,18) = 171 iterations.
For faster convergence, you might want to try this algorithm:
[...]
Yes! But... I really love to see the CAS in action, even if it reaches its limit in terms of recursivity.
(06-23-2020 06:39 PM)Gerson W. Barbosa Wrote: [ -> ]...
For faster convergence, you might want to try this algorithm:
...
The following appears to be better (more tests required) . Only one line has been changed. Anyway, here is it again:
Code:
Program Viete;
Uses Crt;
var k: Byte;
n: LongInt;
c: Char;
d, p, t, v: Extended;
begin
ClrScr;
WriteLn('+---+----------------------+----------------------+');
WriteLn('| k | 2*v | ~ pi |');
WriteLn('+---+----------------------+----------------------+');
v := 1;
n := 1;
d := 0;
for k := 1 to 15 do
begin
n := n + n;
d := Sqrt(d + 2);
v := v*d;
t := 6*v*v - 27/10 - 1/(n*n);
p := 2*n/v*(t + 1)/t;
WriteLn('|',k:2,' |',2*n/v:21:18,' |',p:21:18,' |')
end;
WriteLn('+---+----------------------+----------------------+');
c := ReadKey
end.
+---+----------------------+----------------------+
| k | 2*v | ~ pi |
+---+----------------------+----------------------+
| 1 | 2.828427124746190098 | 3.140960508696045357 |
| 2 | 3.061467458920718174 | 3.141593674140638978 |
| 3 | 3.121445152258052286 | 3.141592707164788547 |
| 4 | 3.136548490545939264 | 3.141592654569488290 |
| 5 | 3.140331156954752913 | 3.141592653605653769 |
| 6 | 3.141277250932772868 | 3.141592653590043215 |
| 7 | 3.141513801144301077 | 3.141592653589797153 |
| 8 | 3.141572940367091385 | 3.141592653589793300 |
| 9 | 3.141587725277159701 | 3.141592653589793240 |
|10 | 3.141591421511199975 | 3.141592653589793239 |
+---+----------------------+----------------------+
P.S.:
This yields 1.8 digits per iteration, three times as much when compared to the plain Viète's formula.
P.P.S.:
This new formula is easier to program and will yield slightly more than two digits per term:
\(\pi \approx 2\left ( \frac{4}{3} \times \frac{16}{15}\times \frac{36}{35}\times\frac{64}{63} \times \cdots \times \frac{ 4n ^{2}}{ 4n ^{2}-1}\right )\cdot \left ( 1+\frac{2}{8n+3+\frac{3}{8n+4+\frac{15}{8n+4+ \frac{35}{8n+4 + \frac{63}{\dots\frac{\ddots }{8n+4+\frac{4n^{2}-1}{8n+4}}} }}} } \right )\)
TurboBCD program:
Code:
Program Wallis_Wasicki;
var d, n, i: Byte;
t: Integer;
c, w: Real;
begin
ClrScr;
WriteLn('+---+---------------------+---------------------+');
WriteLn('| N | 2*W | 2*W*C/(C-2) |');
WriteLn('+---+---------------------+---------------------+');
for n := 1 to 8 do
begin
c := 0;
d := 8*n + 4;
w := 1;
for i := n downto 1 do
begin
t := 4*i*i;
w := w*t/(t - 1);
c := (t - 1)/(c + d)
end;
c := c + d + 1;
WriteLn('|',n:2,' |',2*w:20:17,' |',2*c*w/(c - 2):20:17,' |')
end;
WriteLn('+---+---------------------+---------------------+')
end.
+---+---------------------+---------------------+
| N | 2*W | 2*W*C/(C-2) |
+---+---------------------+---------------------+
| 1 | 2.66666666666666666 | 3.14074074074074073 |
| 2 | 2.84444444444444446 | 3.14159848961611078 |
| 3 | 2.92571428571428570 | 3.14159260997123044 |
| 4 | 2.97215419501133786 | 3.14159265392705764 |
| 5 | 3.00217595455690690 | 3.14159265358714120 |
| 6 | 3.02317019200136082 | 3.14159265358981426 |
| 7 | 3.03867362888341912 | 3.14159265358979309 |
| 8 | 3.05058999605551094 | 3.14159265358979325 |
+---+---------------------+---------------------+
In this table
N is the number of terms,
W is the Wallis Product evaluated to
N terms and
C/(C - 2) is the correction factor.
(06-23-2020 09:58 PM)pinkman Wrote: [ -> ]Quote:Yes, but it pales in comparison to the Wallis-Wasicki formula :-)
Here is a quick PPL port of your Wallis-Wasicki implementation:
...
I’m also porting your 314 pi digits algorithm, I had to stop (work...) but it will be ready in a few hours (if I find the time).
Thank you very much for the PPL port!
Perhaps it's time I should get myself a Prime. But I think I will wait until a good arbitrary precision package is available, either built-in or third-party.
Gerson.
The port of the 200-314-997 decimals will need the CAS, which has capabilities for manipulating big integers but not arbitrary precision decimal numbers.
I’ll check if I can avoid using decimals by calculating on a base of estimated_pi * 10^200 (or less, depending on CAS limit).
Need a few hours more!
I wish the extremely fast HP Prime had the LongFloat Library integrated.
- -
VP
(06-23-2020 09:58 PM)pinkman Wrote: [ -> ]Here is a quick PPL port of your Wallis-Wasicki implementation:
[...]
Gerson sent me a PM to thank me for having posted this code to
hpcalc.org, but in fact I did not post anything, I guess Eric did.
Funny, and good idea, but the credits come to Gerson and his continuous fraction quick convergence for Wallis product.