Hi again,
For the purpose of testing the speed and accuracy of Saturn-based HP calculators such as the HP-71B, HP28S, HP42S, and HP-48/49 series (most specially the all-new, super-fast HP49G+), here's the second (and harder) part of my original advanced test suite, which adds the remaining tests 6 to 10 to the already posted half-suite, tests 1 to 5.
I would be very obligued to those of you who own one of those models (or a suitable emulator) and are willing to run these tests on their machines and take due note of the times and accuracy achieved. Ideally, a minimum of two independent tests per model would be needed to confirm the results. Of course, you can try other HP models instead, but the tests are hard enough to require significant computing muscle which usually restricts the choice to the specified models.
That said, these are the tests comprising the second, more difficult half of the suite. Code and results for the HP-71B are provided for all of them, as well as some relevant comments where necessary:
- Test 6 - Complex-valued Matrix operations:
Set up the test by creating a complex-valued matrix Z with the specified dimensions, and fill it up with different random complex values between (0,0) and
(1,1). Then create a second complex matrix W with the same dimensions and perform the specified matrix operations.HP-71B setup code:
OPTION BASE 1
HP-71B code:
COMPLEX Z(15,15), W(15,15)
FOR I=1 TO 15 @ FOR J=1 TO 15 @ Z(I,J)=(RND,RND) @ NEXT J @ NEXT I
MAT W = INV(Z) 185.3 sec
MAT W = Z*Z 107.2 sec.
MAT W = Z+Z 1.8 sec
- Test 7 - Solving a definite integral of an implicit function:
Find X in [0,1] such that
/X
where y(x) is an ultra-radical function (a member of the family of elliptic functions) implicitly defined by:
|
| y(x).dx = 1/3
|
/0
5
using precisions of 1E-3 and 1E-6 for the integral.
y + y = x
HP-71B code:
X=FNROOT(0,1,INTEGRAL(O,FVAR,1E-3,FNROOT(0,1,FVAR^5+FVAR-IVAR))-1/3)
this gives:
X = 0.854136725005 in 433 seconds (precision = 1E-3)
= 0.854138746461 in 771 seconds (precision = 1E-6)
Test 8 - Integrating a recursively-defined function:
Find the value of the definite integral:
/PI
using precisions 1E-3 and 1E-6, where F(x) is recursively defined as follows:
|
I = | F(x).dx
|
/0
if ABS(x) < 1E-6 then F(x) = x
HP-71B code:
else F(x) = 2*F(x/2)*Sqrt(1-F(x/2)^2)
this gives10 DEF FNF(X) @ IF ABS(X) < 1E-6 THEN FNF = X @ END
20 DIM Y @ Y = FNF(X/2) @ FNF = 2*Y*SQR(1-Y*Y) @ END DEF
30 I = INTEGRAL(0, PI, 1E-3, FNF(IVAR))
I = 2.00000022021 in 85 seconds (precision = 1E-3)
= 2.00000000021 in 171 seconds (precision = 1E-6)
Test 9 - Polynomial Solver for roots of High Multiplicity
Find all roots real and complex of:
x^10 + 10*x^9 + 45*x^8 + 120*x^7 + 210*x^6
HP-71B code:
+ 252*x^5 + 210*x^4 + 120*x^3 + 45*x^2
+ 10*x + 1 = 0
10 OPTION BASE 1 @ DIM C(11) @ COMPLEX R(10)
This gives (in 26 seconds) the following computed roots:
20 DATA 1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1
30 READ C @ MAT R=PROOT(C) @ MAT DISP R
R(1) = (-0.999999998810, 0)
All roots should be equal to (-1, 0). As it always happens with roots of high multiplicity (10 in this case), the accuracy deteriorates, but nevertheless all roots are computed accurate to nearly 3 digits, and in some cases even more accurately (for example, R1 is correct to 8 digits).
R(2) = (-0.999998488294, 0)
R(3) = (-0.999973753013, 0)
R(4) = (-1.000003962220, 0)
R(5) = (-0.989330355306, 0)
R(6) = (-1.005339673390, -9.33448992061E-3)
R(7) = (-1.005339673390, -9.33448992061E-3)
R(8) = (-0.994615440783, -9.24298394994E-3)
R(9) = (-0.994615440783, -9.24298394994E-3)
R(10)= (-1.010783214010, 0)
This is so because the polynomial is extremely flat near X = -1 and evaluates to 0 for most any X within (-1.01 .. -0.99). You can globally check the accuracy of the roots by computing their sum and product. In this case, we get
Sum of the roots = (-9.99999999999, 0)
extremely close to the correct values, (-10, 0) and (1,0) respectively.
Product of roots = ( 0.999999999994, 7.00472767309E-15)
- Test 10 - A probabilistic theoretical application
The expected distance D between two random points on different faces of the unit cube is given by the following expression, involving the sum of two quadruple definite integrals:
/1 /1 /1 /1
HP-71B code:
| | | |
D = 4/5 * | | | | Sqrt(x^2 + y^2 + (z-w)^2) .dw.dx.dy.dz
| | | |
/0 /0 /0 /0/1 /1 /1 /1
| | | |
+ 1/5 * | | | | Sqrt(1 + (y-u)^2 + (z-w)^2) .du.dw.dy.dz
| | | |
/0 /0 /0 /0
10 DEF FNF(X,Y,Z,W) = SQR(X*X + Y*Y + (Z-W)*(Z-W))
using a precision of 1E-3 for all integrals, this gives:
20 DEF FNJ(Y,Z,U,W) = SQR(1 + (Y-U)*(Y-U) + (Z-W)*(Z-W))
30 DEF FNG(X,Y,Z) = INTEGRAL(0, 1, 1E-3, FNF(X,Y,Z, IVAR))
40 DEF FNK(Y,Z,U) = INTEGRAL(0, 1, 1E-3, FNJ(Y,Z,U, IVAR))
50 DEF FNH(X,Y) = INTEGRAL(0, 1, 1E-3, FNG(X,Y, IVAR))
60 DEF FNL(Y,Z) = INTEGRAL(0, 1, 1E-3, FNK(Y,Z, IVAR))
70 DEF FNI(X) = INTEGRAL(0, 1, 1E-3, FNH(X, IVAR))
80 DEF FNM(Y) = INTEGRAL(0, 1, 1E-3, FNL(Y, IVAR))
90 D = 4/5*INTEGRAL(0,1,1E-3,FNI(IVAR)) + 1/5*INTEGRAL(0,1,1E-3,FNM(IVAR))
D = 0.92639 (correct to all 5 digits shown)
in 5h 49 min (23,325.95 seconds). [For the purpose of testing emulators, the value returned is, to 12 digits, 0.926385069656]
In case your machine can't compute quadruple integrals conveniently, here is a MonteCarlo-based simulation in HP-71B code:
10 INPUT "NUM. POINTS="; N @ RANDOMIZE 1 @ S=0
running this gives:
20 FOR I=1 TO N
30 X1=RND @ Y1=RND @ Z1=0 @ X2=RND @ Y2=RND @ Z2=1 @ X3=RND @ Y3=0 @ Z3=RND
40 D1=SQR((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)+1)
50 D2=SQR((X3-X1)*(X3-X1)+Y1*Y1+Z3*Z3)
60 S=S+D1/5+4/5*D2
70 NEXT I
80 DISP "PROBABILITY: ";S/N
with N = 10 -> 0.940
which certainly approaches the theoretical value.
N = 100 -> 0.924
Thanks for any and all results you can contribute, and
Best regards from Valentin Albillo