The 71B/Turtle and the 49G+/Hare: The Harder 2nd Half [LONG] Message #1 Posted by Valentin Albillo on 5 Mar 2004, 9:15 p.m.
Hi again,
For the purpose of testing the speed and accuracy of Saturnbased HP calculators such as the HP71B, HP28S, HP42S, and HP48/49 series (most specially the allnew, superfast 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 halfsuite, 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 HP71B are provided for all of them, as well as some relevant comments where necessary:
 Test 6  Complexvalued Matrix operations:
Set up the test by creating a complexvalued 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.
HP71B setup code:
OPTION BASE 1
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
HP71B code:
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

 y(x).dx = 1/3

/0
where y(x) is an ultraradical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E3 and 1E6 for the integral.
HP71B code:
X=FNROOT(0,1,INTEGRAL(O,FVAR,1E3,FNROOT(0,1,FVAR^5+FVARIVAR))1/3)
this gives:
X = 0.854136725005 in 433 seconds (precision = 1E3)
= 0.854138746461 in 771 seconds (precision = 1E6)

Test 8  Integrating a recursivelydefined function:
Find the value of the definite integral:
/PI

I =  F(x).dx

/0
using precisions 1E3 and 1E6, where F(x) is recursively defined as follows:
if ABS(x) < 1E6 then F(x) = x
else F(x) = 2*F(x/2)*Sqrt(1F(x/2)^2)
HP71B code:
10 DEF FNF(X) @ IF ABS(X) < 1E6 THEN FNF = X @ END
20 DIM Y @ Y = FNF(X/2) @ FNF = 2*Y*SQR(1Y*Y) @ END DEF
30 I = INTEGRAL(0, PI, 1E3, FNF(IVAR))
this gives
I = 2.00000022021 in 85 seconds (precision = 1E3)
= 2.00000000021 in 171 seconds (precision = 1E6)

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
+ 252*x^5 + 210*x^4 + 120*x^3 + 45*x^2
+ 10*x + 1 = 0
HP71B code:
10 OPTION BASE 1 @ DIM C(11) @ COMPLEX R(10)
20 DATA 1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1
30 READ C @ MAT R=PROOT(C) @ MAT DISP R
This gives (in 26 seconds) the following computed roots:
R(1) = (0.999999998810, 0)
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.33448992061E3)
R(7) = (1.005339673390, 9.33448992061E3)
R(8) = (0.994615440783, 9.24298394994E3)
R(9) = (0.994615440783, 9.24298394994E3)
R(10)= (1.010783214010, 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).
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)
Product of roots = ( 0.999999999994, 7.00472767309E15)
extremely close to the correct values, (10, 0) and (1,0) respectively.
 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
   
D = 4/5 *     Sqrt(x^2 + y^2 + (zw)^2) .dw.dx.dy.dz
   
/0 /0 /0 /0
/1 /1 /1 /1
   
+ 1/5 *     Sqrt(1 + (yu)^2 + (zw)^2) .du.dw.dy.dz
   
/0 /0 /0 /0
HP71B code:
10 DEF FNF(X,Y,Z,W) = SQR(X*X + Y*Y + (ZW)*(ZW))
20 DEF FNJ(Y,Z,U,W) = SQR(1 + (YU)*(YU) + (ZW)*(ZW))
30 DEF FNG(X,Y,Z) = INTEGRAL(0, 1, 1E3, FNF(X,Y,Z, IVAR))
40 DEF FNK(Y,Z,U) = INTEGRAL(0, 1, 1E3, FNJ(Y,Z,U, IVAR))
50 DEF FNH(X,Y) = INTEGRAL(0, 1, 1E3, FNG(X,Y, IVAR))
60 DEF FNL(Y,Z) = INTEGRAL(0, 1, 1E3, FNK(Y,Z, IVAR))
70 DEF FNI(X) = INTEGRAL(0, 1, 1E3, FNH(X, IVAR))
80 DEF FNM(Y) = INTEGRAL(0, 1, 1E3, FNL(Y, IVAR))
90 D = 4/5*INTEGRAL(0,1,1E3,FNI(IVAR)) + 1/5*INTEGRAL(0,1,1E3,FNM(IVAR))
using a precision of 1E3 for all integrals, this gives:
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 MonteCarlobased simulation in HP71B code:
10 INPUT "NUM. POINTS="; N @ RANDOMIZE 1 @ S=0
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((X2X1)*(X2X1)+(Y2Y1)*(Y2Y1)+1)
50 D2=SQR((X3X1)*(X3X1)+Y1*Y1+Z3*Z3)
60 S=S+D1/5+4/5*D2
70 NEXT I
80 DISP "PROBABILITY: ";S/N
running this gives:
with N = 10 > 0.940
N = 100 > 0.924
which certainly approaches the theoretical value.
Thanks for any and all results you can contribute, and
Best regards from Valentin Albillo
