[VA] SRC #013 - Pi Day 2023 Special

[J-F Garnier]

So this time we got two solutions from Valentin, plus as often interesting comments and references.
The first solution is nice and quite unexpected, but my preference goes to the second solution that has no limitation due to memory.

My last solution posted just a few hours before Valentin's final solutions probably didn't catch much attention.
Yet I was quite happy with the tricks I used to speed up my code by almost a factor of 2.

One trick was to skip the unnecessary evaluation of the Möbius function for all multiples of 4, i.e. one number over 4 or a 25% gain.
But this was not enough to beat Valentin's timings, and even more with his last 5% optimization of his fastest version.

So here is another attempt, pushing the optimisation further towards speed with just a slightly larger code, but no huge memory requirements.

5 ! SRC13 verE
10 INPUT N
15 T=TIME @ M=(SQR(N)+3) DIV 4*4 @ S=0
20 FOR I=M TO 5 STEP -1
25 I=I-1 @ R=I @ C=-1
30 P=PRIM(R) @ IF NOT P THEN S=S+C*(N DIV (I*I)) ELSE IF MOD(R,P*P) THEN C=-C @ R=R/P @ GOTO 30
35 I=I-1 @ R=I/2 @ C=1
40 P=PRIM(R) @ IF NOT P THEN S=S+C*(N DIV (I*I)) ELSE IF MOD(R,P*P) THEN C=-C @ R=R/P @ GOTO 40
45 I=I-1 @ R=I @ C=-1
50 P=PRIM(R) @ IF NOT P THEN S=S+C*(N DIV (I*I)) ELSE IF MOD(R,P*P) THEN C=-C @ R=R/P @ GOTO 50
55 NEXT I
60 S=S-(N DIV 9)-(N DIV 4)+N
65 T=TIME-T @ DISP N;S;SQR(6*N/S);T


Execution times are now closer ... no, better (slightly) than Valentin's results :-)

 N      Count  Timings (HP-71B)
  12345   7503 7.4s
 100000  60794 24s
 567890 345237 64s
1000000 607926 87s

More results on Emu71/Win, in fast and Authentic modes:

N         Count     PI approx.  Timings (Emu71 fast, auth.)
1E7       6079291  3.14158749068 0.3s 5min11s
1E8      60792694  3.14159307180 0.9s 18min8s
1E9     607927124  3.14159259637 3.1s 63min1s
1E10   6079270942  3.14159267337  11s N/A
1E11  60792710280  3.14159265115  39s N/A
1E12 607927102274  3.14159265250 144s N/A

Do we have to stop at 1E12?
No! The computing loop uses numbers from 2 to SQR(N) which is no problem.
Just the count must be carried carefully, starting from smallest terms, to the highest that can exceed 1E12.
Also some operations must be written carefully, for instance replacing IP(N/X) with N DIV X.
This explains the changes in the program above.

Of course, the final count will have only 12 significant digits, but if is correct within a few ULPs as expected, then we can get an approximation of pi with the same accuracy.

Results on Emu71/Win full speed only:

N         Count         PI approx.  Timings (Emu71 fast)
1E13 6.07927101830E12 3.14159265365 9min27s
1E14 6.07927101860E13 3.14159265357 39min
1E15 6.07927101854E14 3.14159265359 166min

J-F