[VA] SRC#006- Pi Day 2020 Special: A New Fast Way to Compute Pi

[J-F Garnier]

Hi Valentin,

It's a pleasure to read you again in the forum!

Your iterative algorithm is very clever, and is (as always with your contributions) the opportunity to think and learn new things.

So basically, your iterative algorithm solves X=X+SIN(X) around X=3, that is SIN(X)=0 and the solution is of course X=PI.
Now, let me discuss your claim that it is an efficient way to calculate PI:

You anticipated the objections
> Now you may be thinking something like: "Hey, you're using a trigonometric function, namely sin(x), to help compute the value of Pi and that's circular [pun intended]. Come to that, you might simply compute Pi = 4*arctan(1), say."
and you justified your claim with:
>we only need a fast sin(x) custom implementation, which only needs to be able to compute sin(x) for x in the very short range [3 .. Pi], say. Pi itself isn't needed in the implementation as no argument reduction is ever performed, which a general implementation of sin(x) would require.

In the algorithms used by HP (and probably by others, but I know better the HP ones), there is an explicit argument reduction that used a stored value of PI. In the Saturn implementation, the PI value used for argument reduction is stored on 31 digits, as you known, and this is that permits to calculate SIN(3.14159265358) accurately as 9.79323846264E-12 (the next 12 digits of PI).

The fundamental question is: is it possible to efficiently calculate the SIN(X) function around X=PI, actually closer and closer to PI at each iteration, without using argument reduction (that would require PI itself).
You mentioned a reference to an article that states that sin(x) can be performed in O(M(n)log(n)) operations, but doesn't it assume explicitly or implicitly that the argument is reduced to a small range such as [0..PI/4], or if not that the PI value is known with the same n-digits precision?

So here is the challenge I would like to propose:
write an efficient program that calculates, for instance, the 2nd iteration starting from X=3+SIN(3)=3.141120008... and returns the result accurate to let's say 11 digits as: 3.1415926535(7)
without using the PI value itself (so don't use trig operations).
Here is my attempt on the 71B, based on the sin expansion series with a few refinements: the terms are stored in an array, then summed with 15-digit accuracy using the Math ROM DOT operation, starting from the smallest terms.
Even with this refinements, I had to sum terms of the sin(x) expansion up to x^21 (11 terms).

10 !
20 OPTION BASE 1
30 N=100
40 DIM A(N),B(N)
60 X=3 @ X=X+SIN(X)
70 S=0
80 X1=1
90 X2=X*X
100 MAT A=ZER
110 K=N
120 FOR I=1 TO 21 STEP 2
130 A(K)=X1/FACT(I)
140 K=K-1
150 X1=-X1*X2
160 NEXT I
170 MAT B=(X)
180 S=DOT(A,B)
190 DISP X;S;SIN(X)
200 DISP X+S

> RUN
3.14112000806 4.7264551803E-4 4.72645512196E-4
3.14159265358

J-F