[Scilab-users] Strange floating point issue below %eps

classic Classic list List threaded Threaded
7 messages Options
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Strange floating point issue below %eps

Dear co-Scilabers,

I met a strange Scilab's bug this week-end. But today, i tried with Octave, Matlab2016 and R, and i found the same strange behavior. So, either i am missing something, or the bug affects all these languages in the same way. It is reported @ http://bugzilla.scilab.org/15276

In a few words, here it is:

The mantissa of any decimal number is recorded on 53 bits (numbered #0 to #52) + 1 bit for the sign.
This relative absolute accuracy sets the value of the epsilon-machine :
--> %eps == 2^0 / 2^52
 ans  =
  T
From here, it comes, as expected:

--> 2^52 + 1 == 2^52
 ans  =
  F

--> 2^53 - 1 == 2^53
 ans  =
  F

--> 2^53 + 1 == 2^53
   // (A)
 ans  =
  T

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!
But this occurs only when THIS bit #0 is set. For higher bits (hebelow the #1 one), we get back to the expected behavior:
--> (2^53 + 2^2) + 1 == (2^53 + 2^2)
 ans  =
  T

So, when the bit #0 is set and we add a value on the bit "#-1", the language somewhat behaves as if there was a "withholding" value on the bit #0, and seems to toogle it.
Is is a part of any IEEE floating point convention, or am i missing anything, or is it a true bug?
Again, R, Octave and Matlab behave exactly in the same way...

Looking forward to reading your thoughts,

Samuel





_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Antoine Monmayrant Antoine Monmayrant
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

Hello Samuel,

Is this some IEEE standard overflow wraparound?

For information, Julia (0.6) does not behave exactly in the same way.
Both 'print((2^53 + 2^2) + 1 == (2^53 + 2^2))' and  'print((2^53 + 2^1)
+ 1 == (2^53 + 2^1))' return false.

Antoine
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Antoine Monmayrant Antoine Monmayrant
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

In reply to this post by Samuel GOUGEON

Houps, my bad, julia does exactly the same when using Float64, not Int64:


println((2.0^53 + 2.0^2) + 1.0 == (2.0^53 + 2.0^2))

true

println((2.0^53 + 2.0^1) + 1.0 == (2.0^53 + 2.0^1))

false


It has to be some normal IEEE thing but I did not manage to find the proper reference for it.


Antoine


Le 18/09/2017 à 20:54, Samuel Gougeon a écrit :

Dear co-Scilabers,

I met a strange Scilab's bug this week-end. But today, i tried with Octave, Matlab2016 and R, and i found the same strange behavior. So, either i am missing something, or the bug affects all these languages in the same way. It is reported @ http://bugzilla.scilab.org/15276

In a few words, here it is:

The mantissa of any decimal number is recorded on 53 bits (numbered #0 to #52) + 1 bit for the sign.
This relative absolute accuracy sets the value of the epsilon-machine :
--> %eps == 2^0 / 2^52
 ans  =
  T
From here, it comes, as expected:

--> 2^52 + 1 == 2^52
 ans  =
  F

--> 2^53 - 1 == 2^53
 ans  =
  F

--> 2^53 + 1 == 2^53
   // (A)
 ans  =
  T

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!
But this occurs only when THIS bit #0 is set. For higher bits (hebelow the #1 one), we get back to the expected behavior:
--> (2^53 + 2^2) + 1 == (2^53 + 2^2)
 ans  =
  T

So, when the bit #0 is set and we add a value on the bit "#-1", the language somewhat behaves as if there was a "withholding" value on the bit #0, and seems to toogle it.
Is is a part of any IEEE floating point convention, or am i missing anything, or is it a true bug?
Again, R, Octave and Matlab behave exactly in the same way...

Looking forward to reading your thoughts,

Samuel






_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Antoine Monmayrant Antoine Monmayrant
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

Just digging a bit further in julia:

println("2.0^53+2.0^1   : $(bits(2.0^53+2.0^1))")
println("2.0^53+2.0^1+1 : $(bits(2.0^53+2.0^1+1))")
println("")
println("2.0^53+2.0^2   : $(bits(2.0^53+2.0^2))")
println("2.0^53+2.0^2+1 : $(bits(2.0^53+2.0^2+1))")

gives:

2.0^53+2.0^1   : 0100001101000000000000000000000000000000000000000000000000000001
2.0^53+2.0^1+1 : 0100001101000000000000000000000000000000000000000000000000000010

2.0^53+2.0^2   : 0100001101000000000000000000000000000000000000000000000000000010
2.0^53+2.0^2+1 : 0100001101000000000000000000000000000000000000000000000000000010

So the first two are different while the last two are bit identical.

I anyone can point to the specific bit in IEEE spec for this, I'll be happy!


Antoine

Le 19/09/2017 à 08:58, antoine monmayrant a écrit :

Houps, my bad, julia does exactly the same when using Float64, not Int64:


println((2.0^53 + 2.0^2) + 1.0 == (2.0^53 + 2.0^2))

true

println((2.0^53 + 2.0^1) + 1.0 == (2.0^53 + 2.0^1))

false


It has to be some normal IEEE thing but I did not manage to find the proper reference for it.


Antoine


Le 18/09/2017 à 20:54, Samuel Gougeon a écrit :

Dear co-Scilabers,

I met a strange Scilab's bug this week-end. But today, i tried with Octave, Matlab2016 and R, and i found the same strange behavior. So, either i am missing something, or the bug affects all these languages in the same way. It is reported @ http://bugzilla.scilab.org/15276

In a few words, here it is:

The mantissa of any decimal number is recorded on 53 bits (numbered #0 to #52) + 1 bit for the sign.
This relative absolute accuracy sets the value of the epsilon-machine :
--> %eps == 2^0 / 2^52
 ans  =
  T
From here, it comes, as expected:

--> 2^52 + 1 == 2^52
 ans  =
  F

--> 2^53 - 1 == 2^53
 ans  =
  F

--> 2^53 + 1 == 2^53
   // (A)
 ans  =
  T

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!
But this occurs only when THIS bit #0 is set. For higher bits (hebelow the #1 one), we get back to the expected behavior:
--> (2^53 + 2^2) + 1 == (2^53 + 2^2)
 ans  =
  T

So, when the bit #0 is set and we add a value on the bit "#-1", the language somewhat behaves as if there was a "withholding" value on the bit #0, and seems to toogle it.
Is is a part of any IEEE floating point convention, or am i missing anything, or is it a true bug?
Again, R, Octave and Matlab behave exactly in the same way...

Looking forward to reading your thoughts,

Samuel






_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Pinçon Bruno Pinçon Bruno
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

In reply to this post by Samuel GOUGEON
Le 18/09/2017 à 20:54, Samuel Gougeon a écrit :

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!


        Hi all,

   no issue here but simply the round to even rule. Every real number
  should be approximated by the nearest floating point number but in case of
  a number exactly between 2 successive floats the one with an even end
  digit win. You can see this rule as a mean to approximate those ambiguous cases,
  one in two down and one in two up, but has more subtle features in it such
  (it is explained in Knuth I).


  OK that
2^53 + 2^1 is exactly represented
  (it is a double float number). When computing :

          x = (2^53 + 2^1) + 1

  it is not a float but is exactly at the same distance between the two floats :

        
2^53 + 2^1   and 2^53 + 2^2  but this last one ends with a even digit (0)
  so :

         fl( (2^53 + 2^1) + 1) = 2^53 + 2^2

  hth
 Bruno


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Antoine Monmayrant Antoine Monmayrant
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

Thanks a lot Bruno!

I knew I had read something like that but could not remember what it was called.

Cheers,

Antoine

PS: Samuel, do you think you could close the bug by pointing to this "round to even" rule?


Le 19/09/2017 à 11:16, Pinçon Bruno a écrit :
Le 18/09/2017 à 20:54, Samuel Gougeon a écrit :

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!


        Hi all,

   no issue here but simply the round to even rule. Every real number
  should be approximated by the nearest floating point number but in case of
  a number exactly between 2 successive floats the one with an even end
  digit win. You can see this rule as a mean to approximate those ambiguous cases,
  one in two down and one in two up, but has more subtle features in it such
  (it is explained in Knuth I).


  OK that
2^53 + 2^1 is exactly represented
  (it is a double float number). When computing :

          x = (2^53 + 2^1) + 1

  it is not a float but is exactly at the same distance between the two floats :

        
2^53 + 2^1   and 2^53 + 2^2  but this last one ends with a even digit (0)
  so :

         fl( (2^53 + 2^1) + 1) = 2^53 + 2^2

  hth
 Bruno



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: Strange floating point issue below %eps

In reply to this post by Pinçon Bruno
Hello,

Le 19/09/2017 à 11:16, Pinçon Bruno a écrit :
Le 18/09/2017 à 20:54, Samuel Gougeon a écrit :

Now comes the issue:
In (A), the relative difference 1/2^53 is too small (< %eps) to be recorded and to change the number. OK.
Since 1 / (2^53 +2) is even smaller than 1 / (2^53), it should nor make a difference. Yet, it does:

--> (2^53 + 2^1) + 1 == (2^53 + 2^1)
 ans  =
  F

How is this possible ??!


        Hi all,

   no issue here but simply the round to even rule. Every real number
  should be approximated by the nearest floating point number but in case of
  a number exactly between 2 successive floats the one with an even end
  digit win. You can see this rule as a mean to approximate those ambiguous cases,
  one in two down and one in two up, but has more subtle features in it such
  (it is explained in Knuth I).


  OK that
2^53 + 2^1 is exactly represented
  (it is a double float number). When computing :

          x = (2^53 + 2^1) + 1

  it is not a float but is exactly at the same distance between the two floats :

        
2^53 + 2^1   and 2^53 + 2^2  but this last one ends with a even digit (0)
  so :

         fl( (2^53 + 2^1) + 1) = 2^53 + 2^2

  hth
 Bruno

Thank you Bruno for this explanation and confirmation!
Up to before, i was naively considering that all what is smaller than %eps
is ignored. But it is not the case. This rounding rule has some somewhat
puzzling side effects, since as initially illustrated, it leads to sometimes
take into account some relative quantities smaller than some other ones that
are ignored...

So, i guess that whether this algorithm has been chosen, it shall have more
advantages than this kind of pitfall.

Afterwards, i had a look to the Scilab %eps and number_properties help pages:
https://help.scilab.org/docs/6.0.0/en_US/percenteps.html
https://help.scilab.org/docs/6.0.0/en_US/number_properties.html
The IEEE 754 standard is quoted. Digging about it, for instance in
https://en.wikipedia.org/wiki/IEEE_754
https://fr.wikipedia.org/wiki/IEEE_754
is very instructive. There is a whole section about rounding methodS.

Le 19/09/2017 à 12:08, antoine monmayrant a écrit :
> PS: Samuel, do you think you could close the bug by pointing to this "round to even" rule?

I am keeping open and updating the report. IMO the posted %eps "exception" would deserve being documented as an example, either in the %eps or the number_properties page. Moreover, cross-references are missing between both pages. Finally, an external pointer to IEEE 754 would be welcome.
I will propose a revision.

Again thanks for your answers.
Best regards
Samuel


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users