[erlang-questions] erlang float comparison
Angel J. Alvarez Miguel
Hi guys
I need to compare two floats something like 0.99999998... vs 1.000000
and we came across accuracy problems when testing dihedral angles on a molecule...
so we read http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
and wanted to implement something like:
if (*(int*)&f1 < *(int*)&f2)...
when ef1 and f2 are floats... (f1 about 0.9999899.. and f2=1.0
is think we should start with...
<<Padding:16,Myint:64>> = term_to_binary(1.9999976,[{minor_version,1}]).
and let Padding swallow the external format and tag while Myint gets the ieee754 float
But MyInt = 4611686017346523993 instead of 1073741822
these conversion should follow ieee754 and being lexicografic ordered
but ...
term_to_binary(1.99999976,[{minor_version,1}]).
<<131,70,63,255,255,255,191,147,83,89>>
term_to_binary(1.99999988,[{minor_version,1}]).
<<131,70,63,255,255,255,223,201,169,173>>
doesnt seem to be the same thatn you spect to see after reading that page...
..What im doing wrong....?
ieee754 layout....
|
+1.99999976 |
0x3FFFFFFE |
1073741822 |
|
+1.99999988 |
0x3FFFFFFF |
1073741823 |
|
+2.00000000 |
0x40000000 |
1073741824 |
|
+2.00000024 |
0x40000001 |
1073741825 |
|
+2.00000048 |
0x40000002 |
1073741826 |
Thanks!..
/Angel
Hynek Vychodil
_______________________________________________
erlang-questions mailing list
erlang-q...@erlang.org
http://erlang.org/mailman/listinfo/erlang-questions
Hynek Vychodil
Chief Scientists
GoodData
náměstí 28. října 1104/17, 602 00, Brno - Černá Pole
Office: +420 530 50 7704
E-mail: hy...@gooddata.com
Web: www.gooddata.com
Hynek Vychodil
Angel J. Alvarez Miguel
that's is the problem
sinosuke@coreshiba:~> erl
Erlang R15B01 (erts-5.9.1) [source] [64-bit] [smp:4:4] [async-threads:0]
[hipe] [kernel-poll:false]
[QuickCheck] [PropEr] [Erl0MQ]
Eshell V5.9.1 (abort with ^G)
1> 1.000000000000000002 > 1.0.
false
On Lunes, 14 de Mayo de 2012 09:50:57 usted escribió:
> 1> 0.99999998 < 1.000000.
> true
> 2> 1.99999976 < 1.99999988.
> true
>
> What's the problem?
>
> On Mon, May 14, 2012 at 9:41 AM, Angel J. Alvarez Miguel
<cl...@uah.es>wrote:
Hynek Vychodil
1> <<1.000000000000000002/float>>.
<<63,240,0,0,0,0,0,0>>
2> <<1.0/float>>.
<<63,240,0,0,0,0,0,0>>
3> 1.000000000000000002.
1.0
Anyway I think I know what you are about. You want work with numbers
in some interval as they are same. Something like
1> F = fun(X) when is_float(X) -> <<V:32, _/binary>> = <<X/float>>,
<<F/float>> = <<V:32, 0:32>>, <<T/float>> = <<V:32, 16#ffffffff:32>>,
{F, T} end.
#Fun<erl_eval.6.82930912>
2> F(1.0).
{1.0,1.0000009536743162}
3> F(0.99999999).
{0.9999995231628418,0.9999999999999999}
4> F(1.99999976).
{1.9999990463256836,1.9999999999999998}
So you would be able make approximate comparison:
5> Comp = fun(X, Y, N) when is_float(X), is_float(Y), is_integer(N),
N>16, N<64 -> <<XI:N,_/bits>> = <<X/float>>, <<YI:N,_/bits>> =
<<Y/float>>, if XI < YI -> lt; XI =:= YI -> eq; XI > YI -> gt end end.
#Fun<erl_eval.18.82930912>
6> Comp(1.99999976, 1.99999988, 32).
eq
7> Comp(1.99999976, 1.99999988, 48).
lt
lt
--
Hynek Vychodil
Chief Scientists
GoodData
náměstí 28. října 1104/17, 602 00, Brno - Černá Pole
Office: +420 530 50 7704
E-mail: hy...@gooddata.com
Web: www.gooddata.com
Bart van Deenen
> Well we just picked the wrong example
>
>
> that's is the problem
>
> sinosuke@coreshiba:~> erl
> Erlang R15B01 (erts-5.9.1) [source] [64-bit] [smp:4:4] [async-threads:0]
> [hipe] [kernel-poll:false]
>
> [QuickCheck] [PropEr] [Erl0MQ]
> Eshell V5.9.1 (abort with ^G)
> 1> 1.000000000000000002> 1.0.
> false
>
false
2> 1.000000000000000002 == 1.0.
true
3>
So I don't get the problem.
Angel J. Alvarez Miguel
That's what we need, i must confess that i have not earlier expeience with
binaries so needed a bit of thinking just to attack this problem but now im
starting to see it ...
Thanks Angel
Toby Thain
> Well we just picked the wrong example
>
>
> that's is the problem
>
> sinosuke@coreshiba:~> erl
> Erlang R15B01 (erts-5.9.1) [source] [64-bit] [smp:4:4] [async-threads:0]
> [hipe] [kernel-poll:false]
>
> [QuickCheck] [PropEr] [Erl0MQ]
> Eshell V5.9.1 (abort with ^G)
> 1> 1.000000000000000002> 1.0.
> false
>
suggested by others (tolerance won't help you with > or < or the example
above) or whether you need more (or arbitrary) precision.
--Toby
Angel J. Alvarez Miguel
> On 14/05/12 4:47 AM, Angel J. Alvarez Miguel wrote:
> > Well we just picked the wrong example
> >
> >
> > that's is the problem
> >
> > sinosuke@coreshiba:~> erl
> > Erlang R15B01 (erts-5.9.1) [source] [64-bit] [smp:4:4] [async-threads:0]
> > [hipe] [kernel-poll:false]
> >
> > [QuickCheck] [PropEr] [Erl0MQ]
> > Eshell V5.9.1 (abort with ^G)
> > 1> 1.000000000000000002> 1.0.
> > false
>
> It's not clear whether you need an approximate equality test as
> suggested by others (tolerance won't help you with > or < or the example
> above) or whether you need more (or arbitrary) precision.
Sorry folks, im just a bit puzzled for the strangeness of fp operations
I am talking about things like this:
8> 1.0000000000000001 > 1.0.
false
9> 1.0000000000000002 > 1.0.
true
And i wander that there is some method for avoiding such beasts..
Well ill be happy just if i can compare if Something > 1.0 and it just returns
false unless Something is really greater than 1.0 (whatever the accuracy is)
so i can relay on doing things like
if Angle > 1.0 the shutdown_orbital_laser_before_its_too_late().
because we found that:
math:acos(0.99999999999999991).
1.4901161193847656e-8
math:acos(0.99999999999999992).
1.4901161193847656e-8
math:acos(0.99999999999999993).
1.4901161193847656e-8
math:acos(0.99999999999999999).
0.0
math:acos(1.0000000000000001).
0.0
math:acos(1.0000000000000002).
** exception error: bad argument in an arithmetic expression
in function math:acos/1
called as math:acos(1.0000000000000002)
this last case is waht we wanted to avoid, specifically
/Angel
Richard O'Keefe
On 14/05/2012, at 7:41 PM, Angel J. Alvarez Miguel wrote:
>
> Hi guys
>
>
> I need to compare two floats something like 0.99999998... vs 1.000000
> and we came across accuracy problems when testing dihedral angles on a molecule...
>
>
> so we read http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
As the author has said elsewhere,
Years ago I wrote an article about how to do epsilon
floating-point comparisons by using integer comparisons.
That article has been quite popular (it is frequently
cited, and the code samples have been used by a number
of companies) and this worries me a bit, because the
article has some flaws. I’m not going to link to the
article because I want to replace it, not send people
looking for it.
The crucial questions are "what exactly do you want to accomplish,
*semantically*, and what is the problem you are *really* trying
to solve by doing that?"
IEEE formats look simple until you discover that there are, in
reality, endianness issues and worse. For my Smalltalk system I
had to implement a number of functions that I couldn't rely on
being present in C; the code is unpleasantly riddled with
machine-dependent #ifdefs. The way things are stored in memory
and the way they are stored in registers aren't always related
the way you expect.
> is think we should start with...
>
> <<Padding:16,Myint:64>> = term_to_binary(1.9999976,[{minor_version,1}]).
speaks of "THE 64-bit IEEE format", but of course there is no
such animal. The IEEE format is a register format; both floats
and an integers have endianness issues, and it looks as though
here they do not match.
Richard O'Keefe
On 14/05/2012, at 8:47 PM, Angel J. Alvarez Miguel wrote:
> 1.000000000000000002 > 1.0.
That's because in 64-bit arithmetic, 1.000000000000000002 *is* 1.0.
Richard O'Keefe
On 15/05/2012, at 1:53 AM, Angel J. Alvarez Miguel wrote:
>
> Sorry folks, im just a bit puzzled for the strangeness of fp operations
>
> I am talking about things like this:
>
> 8> 1.0000000000000001 > 1.0.
> false
> 9> 1.0000000000000002 > 1.0.
> true
>
> And i wander that there is some method for avoiding such beasts..
Floating point numbers are not mathematical real numbers.
They have a bounded, indeed a fixed, precision.
Single-precision IEEE floats have about 7 digits;
double-precision ones have about 16 digits.
Considering only strictly positive numbers, and being
rather vaguer about the details than I'd like, in order
to keep this simple,
double precision can distinguish x from y
only when |x-y|/max(|x|,|y|) is greater than about 10**-16.
When you enter a decimal number like 1.0000000000000001
it is rounded to the nearest number that can be represented
as a float, and that happens to be 1.0.
>
> Well ill be happy just if i can compare if Something > 1.0 and it just returns
> false unless Something is really greater than 1.0 (whatever the accuracy is)
The problem is that the Something you *type* (or try to calculate)
is not the Something' you *get*; the very best you have any right
to hope for is that Something' is Something rounded to the nearest
representable number. For any single
+ - * / or sqrt, that's what you get. For anything more complicated,
Good Luck (you're going to need it)!
Oh yeah, one other thing: decimal->binary conversion and binary->
decimal conversion are quite hard to get right. I believe Scheme
was the first programming language to require that the results be
as good as possible. C does not. In particular, you should not
be surprised if scanf(printf(x)) != x; since other languages often
use the C libraries for binary<->decimal conversion, you can expect
similar problems in any language that does not explicitly guarantee
best-possible binary<->decimal conversion (and life being what it
is, don't be too surprised to find problems in implementations of
languages that _do_ make such guarantees; implementors may be
unaware of the difficulties).
You know, I have to wonder what the people who fund the development
of autonomous weapons think they are playing at, sometimes...
Toby Thain
> On Lunes, 14 de Mayo de 2012 15:11:52 usted escribió:
>> On 14/05/12 4:47 AM, Angel J. Alvarez Miguel wrote:
>>> Well we just picked the wrong example
>>>
>>>
>>> that's is the problem
>>>
>>> sinosuke@coreshiba:~> erl
>>> Erlang R15B01 (erts-5.9.1) [source] [64-bit] [smp:4:4] [async-threads:0]
>>> [hipe] [kernel-poll:false]
>>>
>>> [QuickCheck] [PropEr] [Erl0MQ]
>>> Eshell V5.9.1 (abort with ^G)
>>> 1> 1.000000000000000002> 1.0.
>>> false
>>
>> It's not clear whether you need an approximate equality test as
>> suggested by others (tolerance won't help you with> or< or the example
>> above) or whether you need more (or arbitrary) precision.
>
> Sorry folks, im just a bit puzzled for the strangeness of fp operations
>
> I am talking about things like this:
>
> 8> 1.0000000000000001> 1.0.
> false
> 9> 1.0000000000000002> 1.0.
> true
As others have pointed out, this is not strange at all. These are the
same values, once precision has been lost by converting the decimal
string to floating point (which has limited precision, see e.g.[1]).
>
> And i wander that there is some method for avoiding such beasts..
>
> Well ill be happy just if i can compare if Something> 1.0 and it just returns
> false unless Something is really greater than 1.0 (whatever the accuracy is)
> so i can relay on doing things like
>
> if Angle> 1.0 the shutdown_orbital_laser_before_its_too_late().
>
To PROPERLY make these comparisons, you need more precision (e.g.
BigDecimal, or fixed point/rational which would be simple with Erlang's
bignum integers).
There are many other issues with floating point, but to deal with the
above < and > comparisons, you simply need sufficient precision.
--Toby
[1] http://en.wikipedia.org/wiki/IEEE_754-1985
Richard O'Keefe
molecules. I find myself wondering what the actual problem is that
the OP is trying to solve, and how much difference between two angles
is *physically* reasonable.
In one document I just found on the web, the H-O-H bond angle is
given as 107.5 degrees in one place and 104.5 degrees in another.
Another document calls it 105 degrees.
My 1st-year chemistry days are long past, but web pages I've checked
today seem not to quote any bond angles to better than 1 part in 1000
precision. And since angles are constrained to [0,one circle), it
looks to me very much as if we have a case for representing bond angles
and dihedral angles as *fixed point* numbers, or at any rate for
rounding them to multiples of say (one circle)/65536 and comparing
the angles as integers.
Angel J. Alvarez Miguel
Thanks for all the comments and all those details. after reading the articles,
that Richard pointed out, and think a bit about what we need i endended with two
solutions that for the most part allow us to get it fixed for the moment:
The problem can be stated as:
Result = case Myfloat > 1.0 of
true -> 1.0;
false MyFloat
end
...
...
math:acos(Result).
after that we came across some values that fail to pass the test > 1.0 but are seen
as greater to the math:acos function rendering a bad arithmetic expression error so finally we manage to do
Result = case trunc(MyFloat) < 1 of
true -> MyFloat;
false -> 1.0
end
...
...
math:acos(Result).
that unless anyone see some hidden mosnter that will bite us the... i think it solves our problem.
"...My 1st-year chemistry days are long past, but web pages I've checked
precision. And since angles are constrained to [0,one circle), it
looks to me very much as if we have a case for representing bond angles
the angles as integers."
an earlier python tool he had developed wher he managed to fix this in a diferent fashion.
So i was working on the OTP part while we found this "feature" but provided that we have enough
time to check and surely we will refactor some code to be carried ogen_fsm,'s maybe...
/Angel
-
Hynek Vychodil
Result = if Myfloat < 1.0 -> Myfloat; true -> 1.0 end,
but little bit more efficient.
--
Hynek Vychodil
Chief Scientists
GoodData
náměstí 28. října 1104/17, 602 00, Brno - Černá Pole
Office: +420 530 50 7704
E-mail: hy...@gooddata.com
Web: www.gooddata.com
Pierpaolo Bernardi
> The problem can be stated as:
>
> Result = case Myfloat > 1.0 of
> true -> 1.0;
> false MyFloat
> end
> ...
> ...
> math:acos(Result).
>
> after that we came across some values that fail to pass the test > 1.0 but are seen
> as greater to the math:acos function rendering a bad arithmetic expression error
If true, this indicates a serious bug. You should report such cases.
I am somewhat skeptical that such a gross bug would lie undetected so
far in erlang, so please supply an actual reproducible example.
Cheers,
P.
Pierpaolo Bernardi
<olop...@gmail.com> wrote:
> On Tue, May 15, 2012 at 11:04 AM, Angel J. Alvarez Miguel <cl...@uah.es> wrote:
>
>> The problem can be stated as:
>>
>> Result = case Myfloat > 1.0 of
>> true -> 1.0;
>> false MyFloat
>> end
>> ...
>> ...
>> math:acos(Result).
>>
>> after that we came across some values that fail to pass the test > 1.0 but are seen
>> as greater to the math:acos function rendering a bad arithmetic expression error
>
> If true, this indicates a serious bug. You should report such cases.
Assuming is not MyFloat < -1, that is.
Cheers
Angel J. Alvarez Miguel
> On Tue, May 15, 2012 at 11:04 AM, Angel J. Alvarez Miguel <cl...@uah.es>
wrote:
> > The problem can be stated as:
> >
> > Result = case Myfloat > 1.0 of
> > true -> 1.0;
> > false MyFloat
> > end
> > ...
> > ...
> > math:acos(Result).
> >
> > after that we came across some values that fail to pass the test > 1.0
> > but are seen as greater to the math:acos function rendering a bad
> > arithmetic expression error
>
> If true, this indicates a serious bug. You should report such cases.
>
> I am somewhat skeptical that such a gross bug would lie undetected so
> far in erlang, so please supply an actual reproducible example.
>
> Cheers,
> P.
It starts smelling a bit, like i messed something up...
Richard O'Keefe
On 15/05/2012, at 9:04 PM, Angel J. Alvarez Miguel wrote:
>
> The problem can be stated as:
>
> Result = case Myfloat > 1.0 of
> true -> 1.0;
> false MyFloat
> end
> ...
> ...
> math:acos(Result).
>
> after that we came across some values that fail to pass the test > 1.0 but are seen
> as greater to the math:acos function rendering a bad arithmetic expression error
The obvious question is *what are these values*?
I've just written a test loop searching for such
values and so far haven't found any.
Is there any difference between emulated and native code?
Why aren't you worried about Myfloat < -1.0?
By the way, this seems a rather strange way to write
Result = if Myfloat > 1.0 -> 1.0 ; true -> Myfloat end
Do you know what trig library your copy of Erlang was linked
with? Here's the relevant code from fdlibm:
hx = __HI(x);
ix = hx & 0x7fffffff;
if (ix >= 0x3ff00000) { /* |x| >= 1 */
if (((ix - 0x3ff00000) | __LO(x)) == 0) { /* |x|==1 */
if (hx > 0)
return 0.0; /* acos(1) = 0 */
else
return pi + 2.0 * pio2_lo; /* acos(-1)= pi */
}
return (x - x) / (x - x); /* acos(|x|>1) is NaN */
}
There's no way _that_ could claim an IEEE double was greater
than 1.0 if it wasn't.
Angel J. Alvarez Miguel
While presenting this "case" these days we also have changed some other code so
(right now) we are.. ah... unable to reproduce such scenario properly, im rechecking
the whole thing (i said we dont have any versioning system?)
Also, asking my workmate about this, results also in some vague memories about acos
failing on a number around 1.0-e8 still we can detail exactly what the number was.. d'oh
So we stated to log all runs after this fiasco in order to honor with real data our (future) claims...
a few quick checks confirm that you all are right while by all means i must be wrong...
so I must state that there isnt erlang fault nor ieee floats...
...in the end is probably my fault, we messed up something, i agree it is rather dificult to find
such values passing the test but failing the acos fun.
But after all and being unsure about floating point behaviour (its probably the first time I try to work
on such numerical thing after many years of boring sysadmin shell scripting...) that i wrote this
on the list....
Sorry for bothering you all, nevertheless, such coverage is one of the things i really appreciate,
every time the list machinery starts growling, we (me) learn a lot of pretty interesting things..
/Angel
-
Richard O'Keefe
On 17/05/2012, at 1:20 AM, Angel J. Alvarez Miguel wrote:
> Also, asking my workmate about this, results also in some vague memories about acos
> failing on a number around 1.0-e8 still we can detail exactly what the number was.. d'oh
then acos(epsilon) approx= pi/2 - epsilon. On the other hand, since pi
is not machine-representable, this might well *look* like a hard case.
One of the features I love about Ada is that the trig functions have
an optional argument saying how much a full circle is, so that you
can express angles like pi/2 precisely as "0.25 of a full circle.
Solaris has a full set of {sin,cos,tan,asin,acos,atan,atan2,sincos}pi()
functions with an implicit multiplication or division as appropriate by
an infinitely precise pi, which is nearly as good.
When it comes to inverse trig functions, I am leery of anything but
two-argument arc-tangent. I beg you to examine every use of
math:acos/1, math:asin/1, and math:atan/1 and very carefully verify
that it is OK in that particular case to get only half the circle
back and that the half you get is in fact the half you want and to
DOCUMENT in the code what you expect and to test that you were right.
I think it's fair to say that floating point is not only weirder than
we imagine, it's weirder than most of us *can* imagine.