HP PRIME not raising overflow/underflow may cause critical errors

[robve]

(03-19-2021 12:27 PM)Hollerith Wrote:  This behaviour goes all the way back to the HP38G from which the Prime has evolved.

That is disappointing. I was hoping there would be a flag to set to raise overflow exceptions.

Writing safe and reliable numerical subroutines in HPPL for the HP PRIME is nearly impossible without internal support for raising (and catching) overflow exceptions. It does not suffice to check that all parameter values are within range for an algorithm. Any line in the code where an expression may overflow must be guarded by an IF-THEN too.

Let's see if and how HP PRIME handles overflow in one of its built-in functions:

horner([1,-1,-1,-1,-1,-1,-1,-1,0],1e38)
1E304
horner([1,-1,-1,-1,-1,-1,-1,-1,0],1E39)
+Inf

Looks good. It is better to be safe than sorry to evaluate the polynomial for values that are too big. However, this is a bit restrictive, because evaluating this 8th degree polynomial with Horner's rule does not overflow for X=1E39:

1E39^8
1E312

Apparently horner performs a quick check on the magnitude of the list's leading coefficients to return +Inf. Quick checks like these are typically overly conservative. The result is actually 1E312:

Define poly
X^8-X^7-X^6-X^5-X^4-X^3-X^2-X
poly(1E38)
1E304
poly(1E39)
1E312

We expect the same with Horner's rule:

Define polyh
X*(-1+X*(-1+X*(-1+X*(-1+X*(-1+X*(-1+X*(-1+X)))))))
polyh(1E38)
1E304
polyh(1E39)
1E312

This checks out. With our polyh version we're not limited to X=1E38 and we can evaluate the polynomial all the way up to X=1E62:

polyh(1E62)
1E496

Now we can also do more interesting things like checking that for larger X the ratio polyh(X)/poly(X/10) is indeed very close to 1E8 in the limit:

polyh(1E21)/poly(1E20)
100000000
polyh(1E51)/poly(1E50)
100000000
polyh(1E61)/poly(1E60)
100000000

However

polyh(1E63)/poly(1E62)
9999.99999999

We get in trouble with division, which in this case internally evaluates 9.99999999999E499 (overflow) and then wrongly assumes Y/9.99999999999E499 is nonzero for any Y>=10

1/(1E250^2)
0
10/(1E250^2)
1E-499

The same, written slightly differently:

10*(1E-250^2)^-1
0

Oh really? Are you kidding me?

Not only are multiplication and division broken when overflow occurs in the denominator, subtraction (and therefore addition) are broken too:

Define polybig
X^8-9.9999E58*X

For any X>=1E59 we expect polybig(X)>polybig(X/10), but:

polybig(1E62)>polybig(1E61)
1
polybig(1E63)>polybig(1E62)
0

The second case overflows at X=1E63.

IEEE standards for floating point impose reasonable requirements regardless of binary coded or BCD representations. These requirements do not appear to be met. Strangely, the HP PRIME has Inf, so it seems natural to expect +Inf for positive overflow and -Inf for negative overflow. Note that NaN is required to reduce Inf-Inf to NaN, as an example. I will leave it aside for further discussion if +0 and -0 should be considered as well.

Now, I should say that I am honestly very impressed with the HP PRIME. It's capabilities are superior to other graphing calculators IMHO. However, this is a glaring flaw that cannot be overlooked or blamed on historical decisions.

- Rob