The Museum of HP Calculators

HP Forum Archive 17

 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #1 Posted by Egan Ford on 24 Apr 2007, 11:32 p.m. The 50g, 48GX, and the 71B all have the same Saturn processor (50g emulated). All 3 have the same machine epsilon (the smallest n where n+1=1, n>0) and the same random number generator. First generate a large (order 20+) random square matrix and vector with values ranging from -1/2 to 1/2. Next solve for x (Ax=b): 50g/48GX: ```b A / (x on stack) ``` 71B: ```MAT x=SYS(A,b) ``` Compare x on both systems. They are exactly the same, as are A and b (assuming the random seeds are the same). Next check the accuracy of the matrix operations (Ax-b): 50g/48GX: ```A x * b - (Z on stack) ``` 71B: ```MAT Z=A*X @ MAT Z=Z-B ``` Compare Z on both systems. Some of the values differ by 1 ULP. Why? Output of my test: ``` 50g/48GX 71B --------------- --------------- .000000000001 .000000000001 -.000000000007 -.000000000007 .000000000002 .000000000002 -.000000000009 -.000000000009 .000000000007 .000000000007 -.000000000005 -.000000000005 -.000000000003 -2.9E-12 != -.00000000001 -.00000000001 .000000000009 .000000000009 .000000000006 .000000000006 .000000000015 .000000000015 .000000000005 .000000000005 .00000000001 .00000000001 -.00000000001 -.00000000001 -1.3E-12 -1.4E-12 != 2.1E-12 2.1E-12 -.000000000017 -.000000000016 != -.000000000012 -.000000000012 .000000000001 .000000000001 .000000000004 .000000000004 ``` 50g/48GX code (put N on stack. A, b, x will be saved): ```%%HP: T(3)A(R)F(.); \<< MEM DROP MEM TICKS 0. 0. 0. 0. 1. \-> N M1 T1 M2 T2 T3 S E \<< 1. RDZ 1. N FOR I .5 RAND - NEXT N \->ARRY 1. N FOR I 1. N FOR J .5 RAND - NEXT N \->ARRY NEXT N ROW\-> MEM 'M2' STO 'A' STO 'b' STO TICKS 'T2' STO b A / 'x' STO TICKS 'T3' STO x OBJ\-> OBJ\-> DROP \->LIST \GSLIST 'S' STO 0. FIX "ORDER: " N + "EST MEM USED: " N N 1. + * 8. * + "MEM USED: " M1 M2 - + "FREE MEM: " M2 + 2. FIX "ARRAY TIME: " T2 T1 - B\->R 8192. / + "SOLVE TIME: " T3 T2 - B\->R 8192. / + "FLOPS: " 2. 3. / N 3. ^ * 3. 2. / N 2. ^ * + T3 T2 - B\->R 8192. / / \->NUM + "TOTAL TIME: " T3 T1 - B\->R 8192. / + STD "SUM(x): " S + "eps: " WHILE E 1. + \->NUM 1. \=/ REPEAT E 2. / \->NUM 'E' STO END E + 6. FIX "RESIDUAL 1: " A x * b - RNRM E A CNRM N * * / + "RESIDUAL 2: " A x * b - RNRM E A CNRM x CNRM * * / + "RESIDUAL 3: " A x * b - RNRM E A RNRM x RNRM * * / + \>> 13. \->LIST 1. \<< 10. CHR + \>> DOSUBS \GSLIST 'RM.TXT' STO STD \>> ``` 71B code (Enter N when prompted. A, B, X, Z in RAM): ```10 STD @ DESTROY ALL @ OPTION BASE 1 @ RANDOMIZE 1 20 H=13 @ DIM L\$(H)[80] 30 M=INT(SQR(MEM/2/8))-1 40 INPUT "N=",STR\$(M);N 50 M1=MEM @ REAL A(N,N),B(N,1) @ M2=MEM @ REAL X(N,1) @ M3=MEM 60 DISP "INITIALIZING ARRAYS..." 70 T1=TIME 80 FOR I=1 TO N @ B(I,1)=.50-RND @ NEXT I 90 FOR I=1 TO N @ FOR J=1 TO N @ A(I,J)=.50-RND @ NEXT J @ NEXT I 100 DISP "SOLVE..." 110 T2=TIME @ MAT X=SYS(A,B) @ T3=TIME @ BEEP 120 S=0 @ FOR I=1 TO N @ S=S+X(I,1) @ NEXT I 130 E=1 140 IF E+1#1 THEN E=E/2 @ GOTO 140 150 REAL Z(N,1) 160 MAT Z=A*X @ MAT Z=Z-B 170 L\$(1)="ORDER:"&STR\$(N) 180 L\$(2)="EST MEM USED:"&STR\$(N*(N+1)*8) 190 L\$(3)="MEM USED:"&STR\$(M1-M3) 200 L\$(4)="FREE MEM:"&STR\$(M2) 210 L\$(5)="ARRAY TIME:"&STR\$(T2-T1) 220 L\$(6)="SOLVE TIME:"&STR\$(T3-T2) 230 FIX 2 240 L\$(7)="FLOPS:"&STR\$(N^3/(T3-T2)) 250 STD 260 L\$(8)="TOTAL TIME:"&STR\$(T3-T1) 270 L\$(9)="SUM(x):"&STR\$(S) 280 L\$(10)="eps:"&STR\$(E) 290 FIX 6 ! STD ! FIX 6 300 V\$=CHR\$(124)&CHR\$(124) 310 U\$=CHR\$(95) 320 L\$(11)=V\$&"Ax-b"&V\$&U\$&"oo/(eps*"&V\$&"A"&V\$&U\$&"1*N)= " 330 L\$(11)=L\$(11)&STR\$(RNORM(Z)/(E*CNORM(A)*N)) 340 L\$(12)=V\$&"Ax-b"&V\$&U\$&"oo/(eps*"&V\$&"A"&V\$&U\$&"1*"&V\$&"x"&V\$&U\$&"1)= " 350 L\$(12)=L\$(12)&STR\$(RNORM(Z)/(E*CNORM(A)*CNORM(X))) 360 L\$(13)=V\$&"Ax-b"&V\$&U\$&"oo/(eps*"&V\$&"A"&V\$&U\$&"oo*"&V\$&"x"&V\$&U\$&"oo)=" 370 L\$(13)=L\$(13)&STR\$(RNORM(Z)/(E*RNORM(A)*RNORM(X))) 380 STD 385 FOR I=1 TO H @ DISP L\$(I) @ NEXT I @ GOTO 400 390 CALL SCROLLW(L\$(),H) 400 ! DESTROY A,B,X,Z,L 410 END ``` Edited: 25 Apr 2007, 12:16 a.m.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #2 Posted by Egan Ford on 25 Apr 2007, 12:54 a.m.,in response to message #1 by Egan Ford More info. If you use the swapdisk MATHROM it will be less like the 50G/48GX. Perhaps my MATHROM has bug fixes (but not enough fixed :-). ``` 50g/48GX 71B (my MATHROM) 71B (swapdisk MATHROM) --------------- ---------------- -------------- .000000000001 .000000000001 -.000000000001 -.000000000007 -.000000000007 -.000000000007 .000000000002 .000000000002 .000000000003 -.000000000009 -.000000000009 -.000000000005 .000000000007 .000000000007 .000000000002 -.000000000005 -.000000000005 -.000000000002 -.000000000003 -2.9E-12 7.E-13 -.00000000001 -.00000000001 -.000000000007 .000000000009 .000000000009 .000000000004 .000000000006 .000000000006 .000000000006 .000000000015 .000000000015 .000000000003 .000000000005 .000000000005 .000000000006 .00000000001 .00000000001 .000000000004 -.00000000001 -.00000000001 -.000000000009 -1.3E-12 -1.4E-12 4.1E-12 2.1E-12 2.1E-12 -4.3E-12 -.000000000017 -.000000000016 -.000000000018 -.000000000012 -.000000000012 -.000000000014 .000000000001 .000000000001 -.000000000002 .000000000004 .000000000004 .000000000008 ```

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #3 Posted by Rodger Rosenbaum on 25 Apr 2007, 2:03 a.m.,in response to message #2 by Egan Ford Is your MATHROM the official hardware version sold by HP, or is it burned into an EEPROM?

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #4 Posted by Egan Ford on 25 Apr 2007, 2:07 a.m.,in response to message #3 by Rodger Rosenbaum Burned, came on EPROM.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #5 Posted by Rodger Rosenbaum on 25 Apr 2007, 2:19 a.m.,in response to message #1 by Egan Ford You don't need to go to anywhere near that much trouble to see the differences. Start with this matrix: ```[[ 888445 887112 ] [ 887112 885781 ]] ``` Invert it. On an HP71 or HP48S, you will get: ```[[ 1279847.88527 -1281771.02151 ] [ -1281771.02151 1283697.0475 ]] ``` On an HP48G, HP49, HP49G+ or HP50G (in approximate mode), you will get: ```[[ 885436.50322 -886766.985569 ] [ -886766.985569 888099.46714 ]] ``` The HP48G and its successors have matrix arithmetic that was reworked by Paul McClellan to use 15 digit arithmetic for internal operations, rounding to 12 digits only occurring after everything is done. The HP71 and HP48S rounded internal operations to 12 digits at intermediate steps in the carrying out of matrix operations.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #6 Posted by GE on 25 Apr 2007, 6:49 a.m.,in response to message #5 by Rodger Rosenbaum On the Casio 9860, you get for the inverse : [ 888781 -887112] [-887112 888445] This answer seems to be exact. How do you explain this ???

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #7 Posted by Rodger Rosenbaum on 25 Apr 2007, 9:13 a.m.,in response to message #6 by GE Quote: On the Casio 9860, you get for the inverse : [ 888781 -887112] [-887112 888445] This answer seems to be exact. How do you explain this ??? Did you make a typo when you copied the 1,1 element of the inverse? I don't know anything about the Casio 9860, so I can only guess about why it gets the exact result. (The HP50 also gets the exact result when it's in exact mode.) What does the 9860 get for the determinant of the original matrix? What does it get for 885781 * 888445 - 7.8696E11 (do the calculation in one step if you can) ? What does it get for 1 / 7 - .142857 ? The HP48G, etc., usually uses the LU decomposition to determine the determinant of a matrix, but for a 2x2 matrix, it does it directly. Perhaps the 9860 treats a 2x2 matrix specially when inverting, and therefore gets an exact result.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #8 Posted by GE on 25 Apr 2007, 10:47 a.m.,in response to message #7 by Rodger Rosenbaum Hi Rodger, this machine is the non-CAS current top of the line Casio. It is the fastest non-CAS machine ever. See Xerces' speed tests in the Articles section. Yes there was a typo, the [1,1] element is 885781. It gets the determinant of the original matrix as exactly 1 (det M -1 gives Zero). 885781 * 888445 - 7.8696E11 returns 7700545. 1/7-.142857 returns 1.42857143e-07, and 1/7-.142857142857 returns 1.43e-13. Casios are notorious for rounding N+Epsilon to N (N integer) when epsilon is small, a bad behavior. You are probably right for the special case of dimension 2.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #9 Posted by Rodger Rosenbaum on 25 Apr 2007, 2:59 p.m.,in response to message #8 by GE Quote: 1/7-.142857 returns 1.42857143e-07, and 1/7-.142857142857 returns 1.43e-13. Casios are notorious for rounding N+Epsilon to N (N integer) when epsilon is small, a bad behavior. You are probably right for the special case of dimension 2. It looks like they're keeping 15 digits in displayed results, and therefore, using at least 15 digits for internal calculations.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #10 Posted by Egan Ford on 25 Apr 2007, 10:02 a.m.,in response to message #5 by Rodger Rosenbaum Quote: You don't need to go to anywhere near that much trouble to see the differences. The 71B code was intended to act as a memory check. That is why it defaults to the largest N that will fit into memory. I wrote this to check some used RAM I purchased, but when the residuals didn't match EMU71 I got concerned. I used the 50G as a tie breaker, but as demonstrated above that didn't help either. As for the memory check, the residuals did match after I used matching Math ROMs in my 71B and EMU71. Quote: The HP71 and HP48S rounded internal operations to 12 digits at intermediate steps in the carrying out of matrix operations. Thanks, that is helpful. I will use a 48SX for comparison. Now please tell me why there are two different Math ROMs for the 71B both with a version of 1A. Are there others? Thanks again.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #11 Posted by Rodger Rosenbaum on 25 Apr 2007, 2:55 p.m.,in response to message #10 by Egan Ford Would you try running the following little program with each of the MATHROM's you have, and posting the results? ```10 OPTION BASE 1 @ DIM C(5) @ COMPLEX R(4) 20 READ C 30 MAT R=PROOT(C) 40 MAT DISP R; 100 DATA 1,4,6,4,1 ``` You may have to type DELAY 9,0 before you run the program. Then you can press enter repeatedly to see each root. edited for grammar :-( Edited: 25 Apr 2007, 8:57 p.m. after one or more responses were posted

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #12 Posted by J-F Garnier on 25 Apr 2007, 3:32 p.m.,in response to message #11 by Rodger Rosenbaum Here are the results (Egan was so kind to send me a copy of "his" Math LEX that I named "math1x" for convenience): ```>lex math1x off @ lex math1a on >run (-1,0) (-.999999501113,0) (-1.00000024944,4.89897948557E-7) (-1.00000024944,-4.89897948557E-7) > >lex math1a off @ lex math1x on > run (-1,0) (-1,0) (-1,0) (-1,0) ``` Impressive! Rodger, it seems you know very well this Math Lex version, can you tell us more? J-F

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #13 Posted by Egan Ford on 25 Apr 2007, 5:05 p.m.,in response to message #12 by J-F Garnier I get the same results.

 Re: 50g/48GX vs. 71B matrix operations--getting slightly different resultsMessage #14 Posted by Rodger Rosenbaum on 25 Apr 2007, 8:54 p.m.,in response to message #12 by J-F Garnier That is a MATHROM that I modified about 20 years ago so that the internal 15-form arithmetic was properly rounded to even just like the 12-form arithmetic. I was curious as to just how much improvement could be had by doing this, and the PROOT function shows the improvement. The 4th order polynomial with equal roots is solved exactly, but not the 5th and higher! The overall improvement is fairly small except for some special cases. This LEX got out to the world by mistake! I didn't take care of the details of NAN's, INF's, etc., so this LEX shouldn't be relied on. I didn't change the version to 1x because I didn't ever expect to release it. It would be good to make that change so that is isn't confused with the real thing.

Go back to the main exhibit hall