HP PRIME not raising overflow/underflow may cause critical errors

[robve]

Apologies if this has already been reported as an issue (I checked with a quick web search but could not find a similar post).

I noticed that an HPPL program I wrote suffered from internal overflow and subsequent underflow. This produced garbage results by the program, which could easily be interpreted as legitimate output.

Consider the following statements in HOME or HPPL:

F := 3.1415926E499
G := 2E250
P := F/G^2

Because G^2 overflows, I would expect F/G^2 to underflow, but that is not what happens. If the overflow or underflow is not caught as a floating point exception, then at least the value should be +0 or NaN as a result of dividing F by +inf. But this is not the case:

F/G^2
0.31415926

I would expect this to be pi/40:


0.07853981...

This means that arithmetic can be broken at times:

H := G/2
F/G^2==1/4*F/H^2
0
F/G^2==F/G/G
0


which are both false (zero) according to the HP PRIME. Worse, the wrong values can still be compared to give the wrong true/false results, because they are not NaN. If they were NaN then all relational comparisons should be false as per IEEE 754.

The HP-71B with DEFAULT OFF produces an exception and with DEFAULT ON gives a warning. Unlike the HP-71B, it appears that common floating point standards (IEEE 754) are not followed by the HP PRIME, or am I wrong?

Has the HP PRIME any fp flags that can be set to raise internal fp overflow exceptions?

Why is there at least no warning message?

Why does IFERR in a HPPL program not catch overflow?

- Rob

[Hollerith]

This behaviour goes all the way back to the HP38G from which the Prime has evolved.

Overflows result in a value of MAXREAL or -MAXREAL depending on sign.
Underflows result in a value of 0.

Neither raise an error and so cannot be caught by IFERR.

Interestingly my HP38G and HP40GS correctly calculate 0^0 as 1. All my other calculators (hard and soft) raise an error.

[Arno K]

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

Overflows result in a value of MAXREAL or -MAXREAL depending on sign.
Underflows result in a value of 0.

Neither raise an error and so cannot be caught by IFERR.

Interestingly my HP38G and HP40GS correctly calculate 0^0 as 1. All my other calculators (hard and soft) raise an error.

I do not want to sound harsh, but 0^0 is undefined and not 1
Arno

[Hollerith]

Arno wrote

Quote:I do not want to sound harsh, but 0^0 is undefined and not 1

This has been argued over for many years (centuries?) with various proofs showing that it is 1 or undefined. I thought the current consensus among mathematicians was that 0^0 = 1.

Whatever is correct, HP has changed it's mind in the past.

[Gene]

Wolfram says it is undefined (today). FWIW.

[Stevetuc]

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

Overflows result in a value of MAXREAL or -MAXREAL depending on sign.
Underflows result in a value of 0.

Neither raise an error and so cannot be caught by IFERR.

Interestingly my HP38G and HP40GS correctly calculate 0^0 as 1. All my other calculators (hard and soft) raise an error.

limit(X^X,X,0) prime does say this is 1

[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

[Hollerith]

There is a lot of information here

https://en.m.wikipedia.org/wiki/Zero_to_...er_of_zero

[Joe Horn]

(03-19-2021 12:27 PM)Hollerith Wrote:  Interestingly my HP38G and HP40GS correctly calculate 0^0 as 1. All my other calculators (hard and soft) raise an error.

The HP 48, and the HP 50g in Approximate mode (or with decimal points typed after the zeros) also return 1 for 0^0. The HP-71 issues a warning (depending on the current trap settings) but also returns 1.

[robve]

A simple textbook example: combinations, aka. nCr binomial coefficients:

Define nCr
X!/Y!*(X-Y)!

Now, this is NOT the most efficient way to compute nCr, but let's put that aside for now (I will comment on that later). We are just defining a classic textbook version of nCr, because we just want to try and apply a function we learned about in a textbook. Sadly, without exceptions being raised, we get garbage:

nCr(270,1)
1

At least an error message would have told us that we better not use this function for large n producing problematic results with low accuracy and overflow.

I am not being facetious. This can be critical. I am aware that most people will not write nCr this way, since we know it is not the most efficient way to compute nCr. The built-in COMB works fine, for example. We can also write some HPPL code to compute nCr the right way. The problem is that user-defined functions and HPPL code cannot be trusted when exceptions aren't raised when you expect them.

- Rob

[robve]

(03-19-2021 02:14 PM)Hollerith Wrote:  Arno wrote

Quote:I do not want to sound harsh, but 0^0 is undefined and not 1

This has been argued over for many years (centuries?) with various proofs showing that it is 1 or undefined. I thought the current consensus among mathematicians was that 0^0 = 1.

Whatever is correct, HP has changed it's mind in the past.

LOL

I was a teenager in the 80s and very lucky to have a Philips P2000T "home computer" and an Apple ][e. I got up at 7am in weekends to program my own BASIC concoctions and some games from magazines and later BASICODE broadcast on the radio. Can you believe it?

The P2000T Microsoft BASIC interpreter did this:

? 0↑0
1

Hmmm... this did not make much sense to me at first, because 0↑n=0 too and both cannot be true at the same time, right?

$$ \left. \begin{array}{l} 0^n = \underbrace{0\cdot 0\cdots 0}_\text{$n$ times} = 0 \cr x^0 = 1 \end{array} \right\} 0^0 = ? $$

This is not a matter of debate or a "current consensus among mathematicians". It is undefined. Just look at the expression above and you will never forget. Period. End of story.

- Rob

[Han]

(03-20-2021 02:28 AM)robve Wrote:  LOL

I was a teenager in the 80s and very lucky to have a Philips P2000T "home computer" and an Apple ][e. I got up at 7am in weekends to program my own BASIC concoctions and some games from magazines and later BASICODE broadcast on the radio. Can you believe it?

The P2000T Microsoft BASIC interpreter did this:

? 0↑0
1

Hmmm... this did not make much sense to me at first, because 0↑n=0 too and both cannot be true at the same time, right?

$$ \left. \begin{array}{l} 0^n = \underbrace{0\cdot 0\cdots 0}_\text{$n$ times} = 0 \cr x^0 = 1 \end{array} \right\} 0^0 = ? $$

This is not a matter of debate or a "current consensus among mathematicians". It is undefined. Just look at the expression above and you will never forget. Period. End of story.

- Rob

Context makes a huge difference in mathematics, though. And depending on the context, it may be 1, or 0, or undefined. It is defined in those contexts by choice in order to have a consistent system. Hence the "consensus." The wikipedia link gives plenty of well-reasoned examples of why they can be taken as any of these. Here is another resource:

https://www.maa.org/book/export/html/116806

EDIT: Also, your explanation is in fact not mathematically sound. For \( 0^n \) your assumption is that \( n \) is a positive integer. And what does it mean to write "zero" copies of the number 0 prior to multiplying? For \( x^0 \), is \(x \) assumed to be a real positive number? Apples and oranges. The former is a discrete problem while the latter is an analytical problem.

[robve]

(03-20-2021 06:05 AM)Han Wrote:  Context makes a huge difference in mathematics, though. And depending on the context, it may be 1, or 0, or undefined. It is defined in those contexts by choice in order to have a consistent system. Hence the "consensus." The wikipedia link gives plenty of well-reasoned examples of why they can be taken as any of these. Here is another resource:

LOL, I do not see any "consensus" in the resources you point to that \( 0^0=1 \) is a universal truth.

Perhaps you are an adherent of a branch of "constructive algebra" that allows you to make up your own rules as you go as long as the algebraic system you're working with is sound. That is perfectly fine. But do not impose your belief system on others when it is universally incorrect. That's the same as saying \( \pi=3 \). It may hold in your custom algebraic system (your "context"), but do not tell everyone that \( \pi=3 \) when they actually want to measure the circumference of a circle in Euclidian space.

(03-20-2021 06:05 AM)Han Wrote:  EDIT: Also, your explanation is in fact not mathematically sound. For \( 0^n \) your assumption is that \( n \) is a positive integer.

Where do I state the condition that n is integer or x has constraints? You are saying that zero cannot be raised by a non-integer n and therefore the "trail of zeros" cannot have a "fraction" of zero? There are no constraints due to universality, i.e. constraints impose exclusions to create a custom algebraic system.

So thank you for actually helping to make my point that there are no constraints on n and x. Hence, the expression I've used is to illustrate in a very simple way that without constraints we both have \( 0^0=1 \) and \( 0^0=0 \), which is inconsistent. Therefore, \( 0^0 \) is undefined. An algebraic system falls apart when counter examples demonstrate that the axioms and rules of the system are inconsistent.

On this highlight, the great mathematician Hilbert was heavily criticized. It is fascinating history.

For all his successes, the nature of his proof created more trouble than Hilbert could have imagined. Although Kronecker had conceded, Hilbert would later respond to others' similar criticisms that "many different constructions are subsumed under one fundamental idea" — in other words (to quote Reid): "Through a proof of existence, Hilbert had been able to obtain a construction"; "the proof" (i.e. the symbols on the page) was "the object". https://en.wikipedia.org/wiki/David_Hilbert

Custom algebraic systems do exists. You can create any such system as you like, including infinitely many custom algebraic systems in which \( 0^0=1 \) as long as you prove that the system is sound. But the point is that these systems are customized to allow \( 0^0=1 \) to be sound. This deviates from our universal system of arithmetic, i.e. such custom algebraic system is not isomorphic to our universal system of arithmetic.

If you claim that \( 0^0=1 \) then it is up to you to prove that it is universally true in our common system of arithmetic. I'm eagerly awaiting...

- Rob

[Albert Chan]

(03-18-2021 08:35 PM)robve Wrote:  This means that arithmetic can be broken at times:

H := G/2
F/G^2==1/4*F/H^2
0
F/G^2==F/G/G
0

Even with IEEE-754, once hit by overflow/underflow, results are questionable.

lua> x = 1e300
lua> x/x/x, x/x^2
1e-300       0

We also have this weird x/x/x ≠ x/x^2 = x/x^2/2 = x/x^2/3 = x/x^3/4 ...

I do like IEEE-754 design that Inf does not turned back to finite value.
(not always, since 1/±Inf underflowed back to ±0.0)

Final result of Inf/NaN implied overflowed intermediate values.

(03-19-2021 04:12 PM)robve Wrote:  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

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:

"horner" is from the CAS side, thus calculating in binary math.
Since 1E312 > 2^1024, horner returned +Inf.

HOME> horner([1,0,-PI*PI], PI)     // to show horner is really doing binary math
6.25277607469ᴇ−13                     // 11/2^44 converted back to 12-digits decimal

(03-20-2021 02:23 PM)robve Wrote:  

(03-20-2021 06:05 AM)Han Wrote:  Context makes a huge difference in mathematics, though. And depending on the context, it may be 1, or 0, or undefined. It is defined in those contexts by choice in order to have a consistent system. Hence the "consensus." The wikipedia link gives plenty of well-reasoned examples of why they can be taken as any of these. Here is another resource:

LOL, I do not see any "consensus" in the resources you point to that \( 0^0=1 \) is a universal truth.

You might want to re-read Han's post again.
The quote actually support your claim that 0^0 is indeterminate.

n = exp(ln(n)/ln(x)*ln(x)) = x^(ln(n)/ln(x))

Taking the limit for any x, we get n back.

x->0: 0^(ln(n)/∞) = 0^0 = n
x->1: 1^(ln(n)/0) = 1^∞ = n

[robve]

(03-20-2021 04:48 PM)Albert Chan Wrote:  You might want to re-read Han's post again.
The quote actually support your claim that 0^0 is indeterminate.

I was primarily responding to the notion of "context" in the post and its relationship to the (non-)constraints on the expression that demonstrates indeterminism.

(03-20-2021 04:48 PM)Albert Chan Wrote:  n = exp(ln(n)/ln(x)*ln(x)) = x^(ln(n)/ln(x))

Taking the limit for any x, we get n back.

x->0: 0^(ln(n)/∞) = 0^0 = n
x->1: 1^(ln(n)/0) = 1^∞ = n

Sure, but taking x to the limit close to 0 is not the same as hitting 0.

With respect to the "power operator" X^Y, I find that there are some fascinating historical anecdotes to implementations in programming languages. If you take programming languages as "context" rather than universal truth then yes, sometimes we get 0^0=1.

Let's start by taking an example in which the power operator does not work, but why? On some the earliest Basic machines, the following code never printed "nine":

X=3
IF X^2=9 THEN PRINT "nine"

Take a close look. Why would "nine" not be printed if 3^2 is clearly 9.

The user manual of these early machines, like the TRS-80 PC-2, warned about the accuracy of the power operator/exponentiation. As we all know, the reason is that the power operator is internally implemented with EXP and LOG:

X^Y = EXP(Y*LOG(X))

EXP and LOG are accurate only up to (but not necessarily including) the last digit of the mantissa, because these functions are typically implemented using their Taylor series expansions. Therefore 3^2 does not compute 9 exactly.

Right away we also realize that X cannot be zero or negative. Allowing X to be non-positive requires adding some conditions and additional code, which is non-trivial, e.g. (-27)^(1/3)=-3 but (-2)^(1/2) is complex.

The engineers of the early home computers had very little ROM space to work with, so any additional code should be kept very small. One simple adddition to improve the speed and accuracy is to use repeated multiplication for a small range of Y, which is faster than EXP and LOG, (or even better, use Exponentiation by Squaring, aka. the "Indian power algorithm", but with lots more code):

if (Y >= 0 && Y <= SOME_MAX && Y == (int)Y) {
  result = 1;
  while (Y-- > 0)
    result *= X
  return result;
}
if (X > 0)
  return exp(Y*log(X));
ERROR!

Note that if Y==0 the result=1 is returned for any X. Therefore, I can imagine that hilarity ensues among mathematicians. However, all that the engineers did is to follow the software engineering requirements specification, which must have been something like this:

implement X^Y for the valid range of real values X and Y

The valid range of exponentiation excludes (X,Y)=(0,0). Therefore, they will argue that the implementation is correct. It is indeed correct as long as programmers do not use 0^0 and then rely somehow on the result being 1.

In similar vein, let's implement factorial N! for integer N:

k = N;
r = 1;
while (k > 0) {
  r *= k;
  k -= 1;
}
return r;

Without checking the range of the factorial function, (-1)! gives 1 which is not correct. It is a "feature" (ahem, a bug) of our program. It expects a valid range for the factorial.

In Axiomatic Semantics of program correctness proofs, we can state this formally with a precondition to the algorithm:

PRE: N >= 0
POST: r = N!

We can then prove the correctness of our algorithm, using the loop invariant \( I \equiv (r\,k! = N! \wedge k \ge 0) \):


\( \{ N\ge 0 \} \)
k = N;
\( \{ k! = N! \wedge k\ge 0 \} \)
r = 1;
\( \{ r\,k! = N! \wedge k\ge 0 \} \)
while (k > 0) {
  \( \{ r\,k! = N! \wedge k\ge 0 \wedge k>0 \} \Rightarrow \{ r\,k\,(k-1)! = N! \wedge k-1 \ge 0 \} \)
  r *= k;
  \( \{ r\,(k-1)! = N! \wedge k-1\ge 0 \} \)
  k -= 1;
  \( \{ r\,k! = N! \wedge k\ge 0 \} \)
}
\( \{ r\,k! = N! \wedge k\ge 0 \wedge k\le 0 \} \Rightarrow \{ r = N! \wedge k=0 \} \Rightarrow \{ r = N! \} \)
return r;

Our program works!

But mind the precondition, please! Otherwise, expect undefined behavior.

- Rob