Floating point didacticism



#7

I have several postings here which indicate that some discussion of what one should expect from the (usually BCD) floating point arithmetic in calculators.

I'll start the discussion with an example to think about. In the V8N1 issue of Datafile, Wlodek Mier-Jedrzejowicz had an article about the "Savage" benchmark. The benchmark is called "Savage" not because it is vicious, but because Bill Savage presented it in Byte magazine.

It is, in HP71 Basic:

5 RADIANS

10 A=1

20 FOR I=1 TO 2499

30 A=TAN(ATN(EXP(LOG(SQR(A*A)))))+1

40 NEXT I

50 PRINT A

Upon first examination, one thinks that the calculator should get a result of 2500.00000000 if it did everything as it should. But this is not so. The HP71 actually gets 2499.99948647.

Should we be disappointed by this result? (All the other Saturn based machines, such as the HP48 get this result also.) What result should we expect, if not 2500.00000000, and why?


#8

The benchmark isn't necessarily telling you anything about the quality of the floating point arithmetic. The result is heavily dependent on the math library's implementation of a certain set of functions. In most function implementations, there are tradeoffs between speed, space and accuracy. Plus you need to know if the implementation is well-behaved and meets its maximum allowable error requirement over its entire range.

#9

Hi Rodger,
[...]
> In the V8N1 issue of Datafile, Wlodek Mier-Jedrzejowicz
> had an article about the "Savage" benchmark.

Thanks for mentioning that article - from January 1989, and a
very interesting set of results. For post-1989 machines, the
following "HP-Solver" code works on the HP19BII and 200LX (The
19BII needs colons instead of commas and the Greek capital
Sigma instead of SIGMA):

0*L(A,1)+SIGMA(I,1,N,1,
0*L(A,1+TAN(ATAN(EXP(LN(SQRT(g(A)*g(A))))))))
+L(B,g(A))-B

2499 N and solving for B gives:
2500.00005173 on the 19BII
2499.999999998308 on the HP200LX

It took 3.5 minutes on the 200LX and 6 minutes on the
19BII. The 19BII result surprised me as the machine is Saturn
based but the result is only .00005 too high whereas the other
Saturn based calcs give the 2499.99948647 which is .0005 too
low.

[...]
> Upon first examination, one thinks that the calculator
> should get a result of 2500.00000000 if it did everything
> as it should. But this is not so. The HP71 actually
> gets 2499.99948647. Should we be disappointed by this
> result? (All the other Saturn based machines, such as the
> HP48 get this result also.) What result should we expect,
> if not 2500.00000000, and why?

I was just pleased to make HP-Solver do the calculation<G>. On
my 200LX QB4 happily gave me 2477.244 in single precision but
when I put the same code in Power Basic, also in single
precision, it just displayed "2500". In DBL precision it gave
"2500.00000000001" and in extended precision "2500" again. And
the 200LX doesn't have a math co-processor. Power Basic seems
quite impressive!

Tony


#10

DataFile V8N8 has a much more elegant Savage Solver equation
than the one I gave. Wlodek (DataFile V8N8 page 25) gives this
beauty for the HP-27S which works fine on the 19BII and
100/200LX:

{T=0*SIGMA(I:2498:0:-1:L(X,1+TAN(ATAN(EXP(LN(SQRT(SQ(X))))))))+X}
You need to remember to start it off by storing 1 in X, and
then solve for T.
DataFile V8 has three articles by Wlodek on this "Savagery" :)
Very interesting.

Tony

#11

Hi,

With an HP86b (capricorn CPU), I got 24499.99942403

5 RAD
7 t=TIME
10 A=1
20 FOR I=1 TO 2499
30 A=TAN (ATN (EXP (LOG (SQR (A*A)))))+1
40 NEXT I
50 PRINT A
55 PRINT HMS$ (TIME -t)
60 END
255553
run
2499.99942403
00:01:22


Olivier

#12

Without working through the details, I might say that it's not necessarily that there's anything wrong with the arithmetic, but that non-linear functions will require greater intermediate precision to get what you seem to be after.

For example, suppose we only had 4-digit floating-point precision, and we square 2500, and take the log. Going to 6 digits, the answer is 15.6481, but we only have 4, so we'll round (correctly) to 15.65. The antilog of that number is 6261936, but we only have 4 digits, so we round, again correctly, to 6,262,000. We've done everything correctly, and yet half of our 4 digits are already wrong. (2500 squared is 6250000.) Taking the square root and correctly rounding, we get 2502. Actually the rounding in this one improved the accuracy, but the last digit is two counts off. If you do this several times, the error will compound.

Or take the cosine and arccosine functions. The cosine of .5625º still rounds to 1.000 if you only have 4 digits. When you take the arccos of that, you correctly get 0.000. Now you have 0.000 versus .5625º. If you take sine of .5735, you get .9999, but then taking the arccosine of .9999, you get .8103, so even the most significant digit is way off.

So does that mean we need 20 or 30 digits though?

I can't speak for others' fields of work; but for my own, generally the loss of precision means that somewhere along the way, it didn't matter that much. If you're taking the tangent of a number near ±90º, what appears to be a huge loss of accuracy in the tangent may translate to an insignificant error in the angle. Keeping 2-digit accuracy in the tangent on your real-life product may require impossibly tight control of the angle, and you probably approached the design entirely wrong.

For a given number of digits' precision, it is impossible to have the accuracy be the same number of digits for certain functions. It doesn't mean anything is wrong. The 10-12-digit precision of most calculators is far beyond what most of us need for real-life applications.

#13

Here is some results of the Savage benchmark with various vintage calculators :

Basic programmable (all Sharp):

PC-1262 : 2499.98731

PC-1403 : 2499.98731

PC-1475 single precision : 2499.987313

PC-1475 double precision : 2499.9999999999999925

PC-E500S single precision : 2500.000056

PC-E500S double precision : 2499.9999999999999999

Keystroke programmable :

HP-11C : 2499.970322

HP-41C : 2499.970322

HP-32SII : 2499.99948647

I don't know if this test is reliable about accuracy, but it gives interesting hints about the algorithm similarity between machines. For example, the seemingly different HP-41CX and HP-11C give the same result.
Performing this test is a lenghty process for these old ladies. It takes one hour and one minute for the HP-41CX to deal with it.

All the best

Pierre Brial

Reunion Island


#14

I wonder what Rodger was looking for when he started this thread :-) The Savage benchmark teaches one about the similarities and differences between various calculators - and it provided an interesting exercise in programming on calculators such as the HP-27S that are programmable only through the Solver. Given correct algorithms, the accuracy of the result depends on the number of digits of precision that calculator (or programming language) uses. I was more interested in the relative speeds of the calculators and hand computers when I wrote those articles. How fast would it be on a 300MHz PDA? ;-)

Wlodek


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP PRIME: derivative at a point Alasdair McAndrew 2 306 11-19-2013, 06:52 AM
Last Post: parisse
  HP28C/S ROM entry point lists update Christoph Giesselink 0 246 01-31-2013, 02:21 PM
Last Post: Christoph Giesselink
  HP-41 ROM Polling Point at FF5 Dan Grelinger 2 315 01-11-2013, 06:00 PM
Last Post: Monte Dalrymple
  33s Decimal Point Matt Agajanian 4 427 11-06-2012, 09:20 PM
Last Post: Matt Agajanian
  wp34S IEEE Floating-Point Conversions Jake Schwartz 12 760 06-08-2011, 10:46 AM
Last Post: Jake Schwartz
  HP 41CX Decimal Point restore Joseph Banki 5 409 04-25-2011, 12:38 PM
Last Post: C.Ret
  Decimal point before an operation Jim Yohe 5 488 01-30-2011, 06:32 AM
Last Post: Chris Randle (UK)
  decimal point budlor 3 322 11-07-2010, 11:42 PM
Last Post: Jim Yohe
  Interesting HP-32E with comma for decimal point Michael de Estrada 2 278 06-17-2010, 12:07 AM
Last Post: Vieira, Luiz C. (Brazil)
  HP-67 not displaying decimal point Etienne Victoria 3 324 03-21-2010, 03:23 PM
Last Post: Etienne Victoria

Forum Jump: