Machine epsilon of the HP Prime

[cdeaglejr]

// program mach_eps HP Prime

// September 17, 2022

// machine epsilon of the HP Prime

// From Wikipedia;

// "Machine epsilon or machine precision is an upper bound
// on the relative approximation error due to rounding in
// floating point arithmetic. This value characterizes
// computer arithmetic in the field of numerical analysis,
// and by extension in the subject of computational science."

/////////////////////////////////////////////////////////////

export mach_eps()

begin

local epsmach;

local one := 1.0;

local two := 2.0;

local u := 1.0;

repeat

u := u / two;

until ((one + u) = one);

epsmach := u;

print();

print("machine epsilon of the HP Prime");

print("\neps = " + epsmach);

print("\n1.0 + eps = " + (1.0 + epsmach));

end;

[rprosperi]

I suspect Home mode and CAS mode have different epsilons...

[EdS2]

I fear there might be a couple of problems with this approach...

I should say that I'm taking machine epsilon to be the same as a unit in the last place (an ULP). If that's wrong, please tell me!

First, I would have thought the correct value of epsilon is not the final u, but the one before that.

Second, on a binary machine with correct rounding, we can add half an ULP and find that it makes a difference, because half an ULP may round upwards.

So, in the case of BBC Basic with a 32 bit mantissa, the program as written returns 2^-33 as the epsilon, whereas I think the right answer is 2^-31.

(Actually, I'm thinking now that even a decimal machine with correct rounding should round up when half an ULP is added... but then again, in this sense I think I'd have to say the majority of HP machines don't tend to round this way. With a decimal machine, halving the candidate epsilon each time might be asking for trouble: at the very least, the final epsilon found will be close to a power of two, but it ought to be an exact power of ten.)

[Albert Chan]

(09-18-2022 08:11 AM)EdS2 Wrote:  BBC Basic with a 32 bit mantissa, the program as written returns 2^-33 as the epsilon,
whereas I think the right answer is 2^-31.

Your "right" answer was the variant definition, machine epsilon = ULP of 1.

1. = 2^31 ULP           → ε = ULP = 2^-31

OP use formal definition, u = ε/2 = 2^-32

Round-to-nearest, halfway-to-even, 1. + u == 1., for any even base

>>> from gmpy2 import *
>>> get_context().precision = 32 # binary32
>>> one = mpfr(1)
>>> one + 2**-31 # 1 + ε
mpfr('1.0000000005',32)
>>> one + 2**-32 # 1 + u
mpfr('1.0',32)

10-digits decimal, correct rounding:

C:\> rpn =10 1e-9 1+
1.000000001
C:\> rpn =10 1e-9 2/ 1+
1

[EdS2]

Ah, thanks, so I may not be fully compliant, but I am a recognised variant!

Quote:This alternative definition is much more widespread outside academia: machine epsilon is the difference between 1 and the next larger floating point number.

This, also from the Wikipedia page, seems to me to cast doubt on this approach:

Quote:The following simple algorithm can be used to approximate the machine epsilon, to within a factor of two (one order of magnitude) of its true value, using a linear search.

(Emphasis added)

[EdS2]

(09-17-2022 05:28 PM)rprosperi Wrote:  I suspect Home mode and CAS mode have different epsilons...

Indeed, I think they should! See for example Joe Horn's explanation:

(10-23-2015 05:41 AM)Joe Horn Wrote:  

(10-23-2015 02:35 AM)luc4as Wrote:  When I try to calculate 3.6^2-4*3.24 in CAS mode, the calculator doesn't return 0, which is the correct answer. Home mode calculates it just fine.

Short answer: It's wrong in CAS for the same reason that .5-.3-.2 doesn't return 0, namely, garbage in garbage out. You THINK you put 3.6 and 3.24 into CAS, but you didn't.

Long answer: It's because Home uses BCD (which can represent 3.6 and 3.24 exactly) and CAS uses binary floating point (which can't). When you type 3.6 in Home, it's EXACTLY 3.6, but when you type 3.6 in CAS, it actually generates this value:

3.599999999999994315658113919198513031005859375 (exactly)

... because that's the largest number less than or equal to 3.6 which is representable with a 48-bit mantissa. For what it's worth, it's stored internally in hex scientific notation as 1.CCCCCCCCCCCCp+1 where "p" means "times 2 to the power of".

...

"If Home and CAS behaved exactly the same, there would be no need for both to exist." -- Joe Horn

[Albert Chan]

For HP Prime, formal definition (maximum relative error) and variant definition (ulp of 1.) are the same.
CAS round number to 53 bits double (round-to-nearest), then chop to 48 bits,

u = maximum relative error = (1 ULP) / (2^47 ULP) = 2^-47

1 = 2^47 ULP, unaffected by rounding-mode.
This probably made variant defintion, ε = ulp of 1., more popular.