[VA] SRC #013 - Pi Day 2023 Special

[Valentin Albillo]

  
Hi, J-F Garnier and EdS2,

J-F Garnier Wrote:So this time we got two solutions from Valentin, plus as often interesting comments and references.

Thanks for your continued appreciation, J-F. I created first solution 1, as doing a sieve was the fastest way to accomplish the goal, but as HP didn't see fit to implement BYTE and/or BOOLEAN arrays in 71 BASIC, like many other high-level languages do, the memory consumption of INTEGER arrays (3 bytes per element) required lots of RAM for large N, and thus I created solution 2, which in the long run gets increasingly slower but has insignificant memory requirements.

J-F Garnier Wrote: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.

Oh yes, it did. I noticed it at once but refrained from replying something outright as I intended to include a version of my solution 1 able to run in puny BASIC dialects without using any of the HP-71 BASIC enhancements provided by the Math and JPC ROMs, as EdS2 requested here for using it on BBC BASIC or other "ordinary BASIC" (his words).

J-F Garnier Wrote: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.

I like the way you never let go until you succeed, and your "tricks" (I'd call them "techniques") are indeed clever and I'm sure interested people will benefit from studying them.

J-F Garnier Wrote:So here is another attempt, pushing the optimisation further towards speed with just a slightly larger code, but no huge memory requirements. [...]

Ditto. However, your use of PRIM means, as I've said, that in the long run that approach can't beat sieve-based methods, because finding a factor (or ensuring that there's none) for large enough numbers is quite costly.

J-F Garnier Wrote:Execution times are now closer ... no, better (slightly) than Valentin's results :-) [...] More results on Emu71/Win, in fast and Authentic modes

Congratulations, but I think you mentioned before that your Emu71/Win is running at ~ 1,500x, I quote:

      "I estimate a peak ratio of ~1500x on my latest Ryzen5 laptop"

while my timings are for a ~ 976x instance, so it would be proper to take that speed difference into account for accurate comparisons of the Emu71/Win timings.

J-F Garnier Wrote:Do we have to stop at 1E12? No! [...] 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.

Impressive ! ... ... and that last approximation for N = 1E15, namely 3.14159265359, is a sight to behold ! Congratulations, again !

---

EdS2 Wrote:I'm thinking, on a platform which lacks the very handy PRIM predicate, one might proceed first with a sieve and then lookup into that for the PRIM calls. [...] It might be that saving the primes in this way is very similar to Valentin's program which records the Möbius values, although in that case I think we still need an implementation of FPRIM.

Here you are, the promised code for your first "related challenge", to run my solution 1 for the HP-71B but using an "unaugmented [sic] 71B, or a similarly ordinary Basic" (your words):

  EDS2    10 lines, 389 bytes

   10  DESTROY ALL @ INPUT T @ SETTIME 0 @ N=IP(SQR(T)) @ GOSUB 70
   20  INTEGER M(N) @ FOR K=1 TO N @ M(K)=1 @ NEXT K @ L=1 @ P=2
  *30  S=P*P @ FOR K=S TO N STEP S @ M(K)=0 @ NEXT K @ FOR K=P TO N STEP P
   40  M(K)=-M(K) @ NEXT K @ L=L+1 @ IF L<=W THEN P=F(L) @ GOTO 30
   50  S=T @ FOR K=2 TO N @ IF M(K) THEN S=S+M(K)*(T DIV (K*K))
   60  NEXT K @ R=SQR(6*T/S) @ DISP T;S;R;ABS(PI-R);TIME @ END

  *70  H=CEIL(N/2) @ INTEGER F(H) @ F(1)=1 @ FOR S=2 TO (SQR(2*H)+1)/2 @ IF F(S) THEN 90
   80  FOR K=2*S*(S-1)+1 TO H STEP 2*S-1 @ F(K)=1 @ NEXT K
  *90  NEXT S @ W=1 @ F(1)=2 @ FOR K=1 TO H @ IF NOT F(K) THEN W=W+1 @ F(W)=2*K-1
  100  NEXT K @ INTEGER F(W) @ RETURN

Description:

● Lines 10 to 60 are essentially the same as in my solution 1 for the HP-71B, with changes (described below in Notes) to adapt it to run in plain-vanilla BASIC dialects. In line 10, before beginning the sieve to compute the Möbius function values en masse, the subroutine beginning at line 70 is first called to obtain at once all prime numbers required

● Lines 70 to 100, nonexistent in solution 1, are now necessary to implement the functionality that the FPRIM keyword allowed there. This subroutine (called just once, before performing the Möbius sieve) computes and stores in integer array F all prime numbers required by performing yet another sieve. Once completed, the prime numbers identified are placed consecutively in array F, which is then redimensioned to reclaim memory, thus leaving more available to later dimension integer array M for the second sieve.

● As was the case with the Möbius sieve, the new prime sieve doesn't use any factoring, multiplications or divisions, for maximum speed. As will be seen in the timings below, this adapted program takes only 24% longer than the original one, which used the assembler keyword FPRIM, so in this battle between my BASIC code vs. JPC Assembler code, I declare myself the winner !

Notes:

This is based on my solution 1 for the HP-71B but with the following changes:

● I've replaced MAT M=CON, WHILE .. END WHILE, FPRIM and RES by functionally equivalent code, and I've deleted USING "image". No keywords from the Math or JPC ROMs remain. The loop at lines 30-40 can be reimplemented using a WHILE .. WEND construct, if available.

● I'm not sure if DIV (integer division) exists in BBC BASIC or other "ordinary" dialects. If it doesn't, A DIV B can be replaced by INT(A/B).

● Same thing with H=CEIL(N/2) at line 70. If CEIL isn't available, it will need to be replaced by trivial equivalent code.

● For convenience while running and timing the code, I've left in the above code the HP-71B BASIC keywords DESTROY ALL, SETTIME 0 and TIME, which should be either deleted or replaced by equivalent code in other BASIC dialects.

● The integer array F is first dimensioned in line 70 as INTEGER F(H), with a variable number of elements H. In BASIC dialects which do not admit dimensioning arrays with a variable number of elements given by an expression, F will probably need to be dimensioned to the fixed maximum number of elements to use. Same thing happens at line 20 INTEGER M(N).

● Also, the same array F is later redimensioned in line 100 as INTEGER F(W), reducing its size (because W is always < H) to reclaim memory. If the plain-vanilla dialect can't redimension arrays, or if it can but it doesn't keep intact the unaffected elements (initializes or destroys the whole array's contents,) then this redimensioning will not be possible and memory won't be reclaimed, which will limit the maximum argument N being input to the program.

The memory savings by redimensioning can be huge. For instance, for N = 100,000,000, redimensioning F reduces it from 5,000 elements to just 1,229 elements.

● As it stands now, the program recomputes the array of prime numbers each time it's run. If there's a RAM disk available or other sufficiently fast way to store and retrieve the array, it could be computed to its maximum size and stored just once, to be retrieved in whole or in part when the program is executed afterwards. This would reduce the running time considerably, depending on the retrieval speed vs. the computation speed.

Results:

Upon running the program with the following arguments on a physical/virtual HP-71B, we get these results and timings:

    N       Count     π approx       |Error|       go71b  Emu71/Win  Physical
                                                   @128x    @976x     HP-71B 
------------------------------------------------------------------------------
   12345      7503  3.14198204634  0.00038939275   0.12"    0.02"        15"
  100000     60794  3.14155932716  0.00003332643   0.34"    0.04"        44"
  567890    345237  3.14158683822  0.00000581537   0.83"    0.11"     1' 46"
 1000000    607926  3.14159550063  0.00000284704   1.11"    0.15"     2' 22"

10000000   6079291  3.14158749068  0.00000516291   3.61"    0.47"     7' 42"
25000000  15198180  3.14159239999  0.00000025360   5.77"    0.76"    12' 19"
33000000  20061593  3.14159276017  0.00000010658   6.65"    0.87"    14' 11"

   1E8    60792694  3.14159307180  0.00000041821  11.71"    1.54"    24' 59"
   1E9   607927124  3.14159259637  0.00000005722  38.00"    4.98"    1h  21'

That's all. Hope it's useful, or at least enjoyable. Tell me how it goes when running it on your beloved BBC BASIC.

V.