Machine epsilon of the HP Prime
09-17-2022, 04:28 PM
Post: #1
 cdeaglejr Member Posts: 63 Joined: Jul 2022
Machine epsilon of the HP Prime
// 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;
09-17-2022, 05:28 PM
Post: #2
 rprosperi Super Moderator Posts: 6,111 Joined: Dec 2013
RE: Machine epsilon of the HP Prime
I suspect Home mode and CAS mode have different epsilons...

--Bob Prosperi
09-18-2022, 08:11 AM
Post: #3
 EdS2 Senior Member Posts: 564 Joined: Apr 2014
RE: Machine epsilon of the HP Prime
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.)
09-18-2022, 12:20 PM (This post was last modified: 09-18-2022 02:41 PM by Albert Chan.)
Post: #4
 Albert Chan Senior Member Posts: 2,496 Joined: Jul 2018
RE: Machine epsilon of the HP Prime
(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
09-19-2022, 06:31 AM
Post: #5
 EdS2 Senior Member Posts: 564 Joined: Apr 2014
RE: Machine epsilon of the HP Prime
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.
09-19-2022, 01:05 PM
Post: #6
 EdS2 Senior Member Posts: 564 Joined: Apr 2014
RE: Machine epsilon of the HP Prime
(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
09-20-2022, 11:24 AM
Post: #7
 Albert Chan Senior Member Posts: 2,496 Joined: Jul 2018
RE: Machine epsilon of the HP Prime
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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)