HP-71B: 1^Inf Error

[J-F Garnier]

I only recently noticed that the HP-71B is raising an error when attempting to calculate 1^Inf.
This result is, at first glance, a bit surprising so I wanted to check on other machines.

The next high-end machines after the 71B were the 28C/S then the series 48 but none were supporting the infinite quantity any more, or other 71B features such as NaN, denormalized numbers and signed zero. Actually these quantities may still result from internal calculations but are not exposed to the user. Maybe because HP thought it was too complicate for most users, with too little benefit.

So we have to jump to the 49G (1999, so 15 years later) to find another HP calculator with direct support of the infinite quantity.
And indeed, the 49G also returns an undefined result ('?') for 1^Inf.
On the other hand, the contemporary TI-89 has no problem to calculate 1^Inf = 1 , but of course returns an undefined result for (-1)^Inf.

They may be good mathematical arguments to consider 1^Inf as undefined, but it's not what I'm mostly interested in.
I'm more curious of the reasons of the HP-71B development team: pure mathematical considerations or implementation constraints?
The choice they made at the time propagated to the next machines for several decades.

J-F

[anetzer]

(06-09-2023 08:17 AM)J-F Garnier Wrote:  This result is, at first glance, a bit surprising so I wanted to check on other machines.

Wolfram Alpha also declares the result undefined.

[Gilles]

(06-09-2023 08:58 AM)anetzer Wrote:  

(06-09-2023 08:17 AM)J-F Garnier Wrote:  This result is, at first glance, a bit surprising so I wanted to check on other machines.

Wolfram Alpha also declares the result undefined.

It's curious. But there must be an mathematical explanation.
With Wolfram Alpha :

1^inf -> Undefine
but
Lim (x->inf, 1^x) -> 1

EDIT : see https://en.wikipedia.org/wiki/Indeterminate_form and https://www.youtube.com/watch?v=I7LjiDDI7ZA

[anetzer]

(06-09-2023 10:06 AM)Gilles Wrote:  It's curious. But there must be an mathematical explanation.

L'Hopital's Rule:

https://web.auburn.edu/holmerr/1617/Text...screen.pdf

After you understood everything that is said there, example 31.7.1 illustrates that HP was perfectly right...

[Albert Chan]

The HP-71B is the first handheld to implement the proposed IEEE floating-point math standard.
Back then, IEEE treatment 1^Inf is not defined yet.

Current IEEE treatment of powers edge cases:

(−1.)^(anything, except  NaN) = 1.
(+1.)^(anything, *even* NaN) = 1.

https://en.cppreference.com/w/cpp/numeric/math/pow

---

But, HP71B IEEE implementation is not perfect ...

>E=1e-9 @ Z=0.
>acos((2,E))
(5.7735026919E-10,-1.31695789692)
>acos((2,Z))
(0,1.31695789692)      ! ???
>acos((2,-Z))
(0,1.31695789692)
>acos((2,-E))
(5.7735026919E-10,1.31695789692)

IEEE correctly implemented, it should be like this:

>>> from cmath import *
>>> E, Z = 1e-9, 0.
>>> acos(complex(2,E))
(5.7735026918962588e-10-1.3169578969248166j)
>>> acos(complex(2,Z))
-1.3169578969248166j
>>> acos(complex(2,-Z))
1.3169578969248166j
>>> acos(complex(2,-E))
(5.7735026918962588e-10+1.3169578969248166j)

[Werner]

Strange..
I can understand 1^Inf =1, but not
(-1)^Inf = 1 or
1^NaN = 1
'feels' wrong.
Werner

[J-F Garnier]

(06-09-2023 11:56 AM)Albert Chan Wrote:  But, HP71B IEEE implementation is not perfect ...

>E=1e-9 @ Z=0.
>acos((2,E))
(5.7735026919E-10,-1.31695789692)
>acos((2,Z))
(0,1.31695789692)      ! ???
>acos((2,-Z))
(0,1.31695789692)
>acos((2,-E))
(5.7735026919E-10,1.31695789692)

Note that the complex inverse trig functions are supported on the 71B only with the MATH version 2B that I wrote. You would get however, the same result on a HP-48, except you can't do acos((2,-0)) because there is no signed zero on the 48.

I remember that the HP-71B MATH ROM developers reported (in the HP Journal I guess), that they had trouble to extend the IEEE concepts to the complex domain. It was very new at the time. So they did their best and this was basically unchanged in the next 28C/S and series 48.
In this particular example, the signed zero is not managed properly in the complex domain. Maybe an additional reason for having abandoned it?

J-F

[Albert Chan]

(06-09-2023 01:03 PM)J-F Garnier Wrote:  In this particular example [acos(2+0j)], the signed zero is not managed properly in the complex domain.
Maybe an additional reason for having abandoned it?

I believe it is asin/acos sign bug in code.
Calculate acos(2+0j) by "hand", I am getting correct answer.

Implementing the complex arcsine and arccosine functions ... (page 326 for acos(z) algorithm),

lua> x, y = 2, 0
lua> signy = signbit(y) and -1 or 1
lua> r = sqrt((x+1)^2 + y^2)
lua> s = sqrt((x-1)^2 + y^2)
lua> a = (r+s)/2
lua> b = x/a
lua> I.new(acos(b), -signy * acosh(a))      -- = acos(2+0j)
(0-1.3169578969248166*I)

sign(y) bug is in cpr.a, which sometimes use -sign(x) instead. (docasc code comment)
Because asin(z) + acos(z) = pi/2, this bug affected both asin and acos.

Code:

*      with sign(y,x)=sign(y) if y<>0, -sgn(x) if y=0

>acos((2,0))
(0,1.31695789692)      ✘
>acos((-2,0))
(3.14159265359,-1.31695789692)      ✔

[Albert Chan]

To be more correct, sign(Im(acos(z))) = sign(-Im(z))

Here is Kahan's acos(z) implementation.
https://www.hpmuseum.org/forum/thread-13...#pid151453

lua> z = I.new(2,0)
lua> a = I.sqrt(1+z)
lua> b = I.sqrt(1-z)
lua> x2 = atan(b:real() / a:real()) * 2
lua> y2 = asinh(a:real()*b:imag() - a:imag()*b:real())
lua> I.new(x2, copysign(y2, -z:imag()))      -- = acos(2+0j)
(0-1.3169578969248166*I)

[J-F Garnier]

To go on from an historical perspective, here is a comparision of the answers of some machines of the end of last century to the seven indeterminate forms described here:

           ------ pre-CAS ------              ------ CAS-based ------
           HP-71B           HP-48SX            HP-49G           TI-89
           (1984)            (1990)               (1999)           (1998)
0/0        error              error                 undef            undef
Inf/Inf    error              n/a                   undef            undef
0xInf      error              n/a                   undef            undef
Inf-Inf    error              n/a                   undef            undef
0^0       1 (warning)     1 (no warning)    0^0=undef     1 with warning
                                                        but 0.^0=1
1^Inf      error              n/a                  undef            1 with warning
Inf^0      1 (warning)     n/a                  undef            1 with warning

The 1^Inf exception could have had a default value of 1 with warning on the 71B, but the designers decided otherwise.

Edit: To clarify my thought: the 71B can operate in DEFAULT OFF mode, where it issues an error for each exception, or in DEFAULT ON mode where it gives a default value (with warning) in certain cases as indicated above.

J-F

[Gilles]

0.99999999…9 ^infinity -> 0
1.00000000…1 ^infinity. -> infinity

[Guenter Schink]

Python on Prime:
1**float('inf') -> 1
-1**float('inf') -> -1
1**float('-inf') -> -1
-1**float('-inf') -> 1

Günter

[EdS2]

Hmm that's all kinds of wrong.

[rprosperi]

(06-09-2023 12:21 PM)Werner Wrote:  Strange..
I can understand 1^Inf =1, but not
(-1)^Inf = 1 or
1^NaN = 1
'feels' wrong.
Werner

I think this proves that Infinity is even, right?

[Werner]

(06-09-2023 11:38 PM)rprosperi Wrote:  I think this proves that Infinity is even, right?

What does MOD(Inf, 2) return? ;-)
Werner

[C.Ret]

HI,

Another very interesting subject that trigger my curiosity, because I found the evaluation of 1^∞ very counter-intuitive.
So I tried to document myself a little bit. None of my math textbooks cover the subject. On the other hand, I found some video which explain well the concern of 1^∞.
One or two samples of tutorial video in english.

I understand better why my HP Prime displays this:


Indeed, as the designers of the HP-71B did, the designers of the HP Prime (its CAS in fact) may have deliberately chosen to return undef because they consider that entering 1.^∞ ( or even 1^∞ ) does not make it possible to determine whether it is an exact unit basis or not.
For other manufacturers, the interpretation may be different, so my Ti-92 II returns 1^∞=1 certainly considering that it is exactly a unit base.


So, the initial question from Jean-François remains:

(06-09-2023 08:17 AM)J-F Garnier Wrote:  I'm more curious of the reasons of the HP-71B development team: pure mathematical considerations or implementation constraints?

[Gilles]

(06-10-2023 10:50 AM)C.Ret Wrote:  HI,

I understand better why my HP Prime displays this:

Thanks C.Ret. I get (obviously) the same result with my 49g+ in exact mode. First time I use the 'lim' command with a list ;D
It proves that if we assert that 1^infinity=1, the whole algebra becomes incoherent.

Infinity is a strange thing. Someone told me one day that all the books in the world were already written in a sequence of decimals of PI. And not just all the books in the world already written, but all the books that can be written.

[rprosperi]

(06-10-2023 08:55 AM)Werner Wrote:  

(06-09-2023 11:38 PM)rprosperi Wrote:  I think this proves that Infinity is even, right?

What does MOD(Inf, 2) return? ;-)
Werner

Seems like an invalid argument to me... (ba dump bump...

[EdS2]

I notice IEEE floating point specifies the result of 0 × ∞ as NaN. If 1^x is computed as exp(log(1)*x) then that would be the same. (Why is it NaN? Perhaps because "The floating point zeros not only represent the real number zero. They also represent all real numbers that would round to something smaller than the smallest subnormal. That is why zero is signed. Even tiny numbers do have a sign if they are not actually zero.")

[Albert Chan]

(06-10-2023 01:47 PM)EdS2 Wrote:  I notice IEEE floating point specifies the result of 0 × ∞ as NaN.
If 1^x is computed as exp(log(1)*x) then that would be the same.

IEEE goal is not always mathematically correct.
It tried to avoid NaN if possible, with a default numerical value to (hopefully) simplify code.

https://www.open-std.org/jtc1/sc22/wg14/...eV5.10.pdf

page 182 Wrote:pow(x,0) is specified to return 1 for any x, because there are significant applications where 1 is
more useful than NaN. pow(f(t),g(t)) approaches 1 in all cases where f and g are analytic
functions and g(t) approaches zero. The result 1 better supports applications where the second
argument is integral. pow(NAN,0) returns 1.0 on the general principle that if a result is independent
of the numerical value of an argument, then that result is appropriate if that argument is a NaN.

(06-09-2023 11:38 PM)rprosperi Wrote:  I think this proves that Infinity is even, right?

Yes!

page 180 Wrote:The result of pow(–2,∞) is +∞, because all large positive floating-point values are even integers.

--> pow(-1, ∞) = pow(1, ∞) = 1

Again, this is not mathematically sound, but only for "convenience"