Re: what's the ULP in Octave

 >> On Dec 1, 2020, at 7:55 AM, 康 珊 <[email protected]> wrote:
 >>
 >> Hi,
 >> I wonder what's the ULP(unit of least precision) for those double 
operations such as pow in octave.


On 12/2/20 6:01 AM, Victor Norton via Help-octave wrote:
 > Try this
 >    octave> x = 1; eps = 1;
 >    octave> while(x + eps > x)
 >    > eps /= 10;
 >    > end
 >    octave> eps
 >    eps =    1.0000e-16
 > Is this the answer?
 >


There is one important change necessary.  Octave's double is binary not 
decimal [1,2], thus it should be "eps /= 2;".  Otherwise the program 
will give the correct answer.

However, there is an Octave function called "eps" [3] which gives the 
result as well.  No need to write own code, but for clarification it can 
be beneficial to run some code.

```````````````````````````````````````````````
x = 1;
e = 1;
i = 0;
while (x + e > x)
   printf ("1 + 2^(%3d) = %s\n", ...
           i, num2hex (x + e))
   e /= 2;
   i--;
endwhile
printf ("1 + 2^(%3d) = %s  <-- 'e' too small\n", ...
         i, num2hex (x + e))
```````````````````````````````````````````````

The output gives a hexadecimal representation of double numbers "1.0 + 
e" for shrinking values of "e":

1 + 2^(  0) = 4000000000000000
1 + 2^( -1) = 3ff8000000000000
1 + 2^( -2) = 3ff4000000000000
     ...
1 + 2^(-51) = 3ff0000000000002
1 + 2^(-52) = 3ff0000000000001
1 + 2^(-53) = 3ff0000000000000  <-- 'e' too small

Thus:

- 2^(-52) == eps() (approx. "2.2204e-16") is the smallest number that 
can be added to "1.0" to "change" this number (the ULP unit in the last 
place).

- 2^(-53) == eps()/2 (approx. "1.1102e-16") is the number that added to 
1.0 is too small to "change" this number.  The precision of double is 
not enough to store the exact result of this addition.

A final note: eps() is *relative* to 1.0.  To obtain eps relative to 
another double number, say 5.0, call "eps (5.0) == eps() * 4 == 2^(-50)"

HTH,
Kai


[1] https://en.wikipedia.org/wiki/Unit_in_the_last_place
[2] https://en.wikipedia.org/wiki/Double-precision_floating-point_format
[3] https://www.octave.org/doc/v6.1.0/XREFeps.html