[VA] SRC #012a - Then and Now: Probability

[J-F Garnier]

(10-16-2022 09:20 AM)Albert Chan Wrote:  Here is optimized C.Ret version, without "neighbor hunting" 2 inner loops.

>RUN
[VA]SRC012A R,S= 30,60
...
67.71     30 60      9.51234350207E-6

From this thread so far, this is the most elegant code, and the fastest!

And perfectly exact too, down to the last place! Well maybe also by a little bit of chance.

Also thanks to the slightly optimized memory required for Q(R+1,CEIL(R/2)+1), it fits in a 24k-RAM HP-75 !

Here are the few simple changes to run your code on the HP-75:
7 OPTION BASE 0 @ REAL W(30,30),P(30,30),Q(31,16)
10 INPUT "[VA]SRC012A R,S= "; R,S @ T0=TIME
20 REDIM W(R,R),P(R,R)
30 REDIM Q(R+1,CEIL(R/2)+1)
35 MAT P=ZER @ MAT Q=ZER @ MAT W=ZER ! vars must be initialized
The rest is unchanged.

HP-75 result is: 9.51234350244E-6

Slightly OT: why is the HP-75 result different, and less accurate (I like to investigate these kind of questions) ?
After all, it's all about additions, multiplications and divisions, and the HP-75 and HP-71 share the same algorithms and 12-digit (15 internally) accuracy.

Well, almost the same algorithms.
The Saturn machine algorithms (from the HP-71) introduce a very small improvement, known as the "round-to-even" or "banker's rounding" rule.
When an internal 15-digit result is rounded to the user 12-digit form, and it ends exactly with ...500, then it is rounded to the closest even value, instead of being rounded upward.
It can be demonstrated for instance with 1/(2^18) = 3.81469726562(500) e-6
HP-71: 3.81469726562e-6 (round to even)
HP-75: 3.81469726563e-6 (round up)

This difference occurs in the first steps of the calculation, and introduces a small bias that is enough to bring a different and less accurate answer at the end.

J-F