[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## Re: what's the ULP in Octave

**From**: |
Kai Torben Ohlhus |

**Subject**: |
Re: what's the ULP in Octave |

**Date**: |
Wed, 2 Dec 2020 11:50:14 +0900 |

**User-agent**: |
Mozilla/5.0 (X11; Linux x86_64; rv:78.0) Gecko/20100101 Thunderbird/78.5.0 |

>> On Dec 1, 2020, at 7:55 AM, 康 珊 <kangshan0910@hotmail.com> 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