Velocity Reviews > Perl > accuracy problem

# accuracy problem

Eric Wilhelm
Guest
Posts: n/a

 11-15-2003
I'm running into some problems where 1+1 != 2.

I've tried compiling 5.8.0 and 5.9.0 with -Duselongdouble and
XS::APItest::have_long_double() returns true, but I still cannot
get the following code to perform correctly:

# test floating-point accuracy
@a = (40.23445,10.7784); # this will fail
# @a = (0,0); # this will succeed
\$r = 3/16;
\$pi = atan2(1,1) * 4;
\$ang = 12 * \$pi / 180;

@s = ( \$a[0] + \$r * cos(\$ang), \$a[1] + \$r * sin(\$ang) );
@f = ( \$a[0] - \$r * cos(\$ang), \$a[1] - \$r * sin(\$ang) );

\$dist = sqrt( (\$f[0] - \$s[0])**2 + (\$f[1] - \$s[1])**2);
\$d = 2 * \$r;
if(\$dist != \$d) {
print "failed: \$dist != \$d\n";
}
else {
print "okay: \$dist == \$d\n";
}
# end

this gives:
failed: 0.374999999999999998 != 0.375

unless @a is set to (0,0)

Recreating the process in C and declaring everything as double gives the
same results, but if you declare as long double it will come up with
matching numbers, thus the quest to get a Perl working with long doubles.

Anyone have similar experiences? Any way to fix this? Rounding is not
a solution, because the contents of @s and @f are the desired product,
and the distance between them needs to be 2 * \$r when computed with long
doubles (i.e. the inaccuracy apparently happens in the addition as
opposed to the distance calculation.)

Thanks,
Eric

Eric Wilhelm
Guest
Posts: n/a

 11-15-2003
On Sat, 15 Nov 2003 15:12:32 -0600, Eric Wilhelm wrote:

I just realized that I had the wrong terminology, though that doesn't
change the problem.

> I'm running into some problems where 1+1 != 2.
>
> I've tried compiling 5.8.0 and 5.9.0 with -Duselongdouble and
> XS::APItest::have_long_double() returns true, but I still cannot get the
> following code to perform correctly:

Forgot to mention that this in on an AMD MP2200 running Linux 2.4.21.

--Eric

Jürgen Exner
Guest
Posts: n/a

 11-15-2003
Eric Wilhelm wrote:
[...]
> this gives:
> failed: 0.374999999999999998 != 0.375
>
> Recreating the process in C and declaring everything as double gives
> the same results, but if you declare as long double it will come up
> with matching numbers, thus the quest to get a Perl working with long
> doubles.
>
> Anyone have similar experiences?

Everybody who tried to do computer numerics without knowing about computer
numerics.

> Any way to fix this?

Yes.
Is is mentioned in any computer numerics class thou shalt never compare
so-called real numbers (which are more aptly called floating point numbers)
for equality. Instead check if they are within an epsilon range of each
other.

Alternatively you can use a package for symbolic mathematic.

For further information see "perldoc -q 999" or your favourite book about
computer numerics

jue

Eric Wilhelm
Guest
Posts: n/a

 11-16-2003
On Sat, 15 Nov 2003 16:28:34 -0600, Jürgen Exner wrote:

>> Any way to fix this?

>
> Yes.
> Is is mentioned in any computer numerics class thou shalt never compare
> so-called real numbers (which are more aptly called floating point
> numbers) for equality. Instead check if they are within an epsilon range
> of each other.

Thanks, but I have no intention of comparing them directly. The problem
is that they will not be within the required epsilon range of each
other if the coordinates are not calculated with 128-bit precision.

I have long been aware of the advice to round numbers before trying to do
a comparison, etc, etc, but this is not a financial application and what
I have to make happen is that @ptB is halfway between @ptA and @ptC so
that an arc centered at @ptB can go through @ptA and @ptC in a piece of
external software which is using 'long double' for floating point numbers.

It seems from the INSTALL and README files in the perl source that perl
should be using 'long double' if configured as such, but the test does
not give the same results as C code which does use 'long double'
directly.

Thanks,
Eric

Eric J. Roode
Guest
Posts: n/a

 11-16-2003
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Eric Wilhelm <(E-Mail Removed)> wrote in
newsan.2003.11.15.15.13.11.640507.1618@something like.sbcglobalDOTnet:

> Anyone have similar experiences?

Oh yeah. Everyone who has ever used floating-point numbers has had this
problem.

> Any way to fix this?

Nope. It's a limitation of the precision of the hardware. It doesn't
matter if you use double-precision numbers, or quadruple, or whatever.
Some numbers simply are not representable in a finite mantissa/exponent
model.

> Rounding is not
> a solution, because the contents of @s and @f are the desired product,
> and the distance between them needs to be 2 * \$r when computed with long
> doubles (i.e. the inaccuracy apparently happens in the addition as
> opposed to the distance calculation.)

Rounding is not a solution? Not even, say, rounding to the 12th decimal
place? Surely nothing you're doing can require that level of precision.

The simple truth is, you cannot reliably use == or != with floating-point
numbers. You must do some sort of "is this close enough" comparison on
your own. You have to decide what constitutes "close enough".

- --
Eric
\$_ = reverse sort \$ /. r , qw p ekca lre uJ reh
ts p , map \$ _. \$ " , qw e p h tona e and print

-----BEGIN PGP SIGNATURE-----
Version: PGPfreeware 7.0.3 for non-commercial use <http://www.pgp.com>

iQA/AwUBP7bBMGPeouIeTNHoEQLf0wCg7A6oQOgdybi75AcRr9+bmr 2PkVwAoJAY
6LcSvscd+oJBSKee0cX71BNA
=U1By
-----END PGP SIGNATURE-----

Eric Wilhelm
Guest
Posts: n/a

 11-16-2003
On Sat, 15 Nov 2003 18:12:55 -0600, Eric J. Roode wrote:

> Rounding is not a solution? Not even, say, rounding to the 12th decimal
> place? Surely nothing you're doing can require that level of precision.
>
>

The issue is getting the calculations to resolve between cartesian
coordinates and trigonometric calculations. If I have an arc that has
been centered between two points and I calculate its endpoints, rounding
each coordinate (to any decimal place), they are not matching up to the
points that were averaged to get to the center point. While I am able to
sidestep the issue, it is causing problems in downstream programs which
are expecting higher precision.

> The simple truth is, you cannot reliably use == or != with
> floating-point numbers. You must do some sort of "is this close enough"
> comparison on your own. You have to decide what constitutes "close
> enough".

I have just realized that maybe I have what I was asking for:

0.375000000000000002 != 0.375 comes from the 5.8.0 custom build
0.375000000000004 != 0.375 comes from the 5.6.1 which has been installed

I guess the 'long double' is only 80 bits (somehow I cooked up the idea
that it was 12, but it is longer than the 64 bits that I was apparently
using before.

Now I have to setup a new installation with all of the libraries that were
being used (this problem came up in the middle of a lot of abstract and
complicated functions) and see if it solves the issue with the downstream
software.

Thanks,
Eric

Jürgen Exner
Guest
Posts: n/a

 11-16-2003
Eric Wilhelm wrote:
> On Sat, 15 Nov 2003 16:28:34 -0600, Jürgen Exner wrote:
>
>>> Any way to fix this?

>>
>> Yes.
>> Is is mentioned in any computer numerics class thou shalt never
>> compare so-called real numbers (which are more aptly called floating
>> point numbers) for equality. Instead check if they are within an
>> epsilon range of each other.

>
> Thanks, but I have no intention of comparing them directly. The
> problem is that they will not be within the required epsilon range of
> each other if the coordinates are not calculated with 128-bit
> precision.

You mean you require 38 (in words: thirtyeight!) digits of precision?
Are you doing some kind of astronomical navigation or so?

Well, ok, in that case you may want to re-read the fifth paragraph in
"perldoc -q 999"
Apparently you missed the advice given there when you read it the first
time.

> I have long been aware of the advice to round numbers before trying
> to do a comparison, etc, etc,

Ok, sorry, that wasn't clear from your first posting.
Again, please see "perldoc -q 999".

jue

Jürgen Exner
Guest
Posts: n/a

 11-16-2003
Eric Wilhelm wrote:
> On Sat, 15 Nov 2003 18:12:55 -0600, Eric J. Roode wrote:
>> The simple truth is, you cannot reliably use == or != with
>> floating-point numbers. You must do some sort of "is this close
>> enough" comparison on your own. You have to decide what constitutes
>> "close enough".

>
> I have just realized that maybe I have what I was asking for:
>
> 0.375000000000000002 != 0.375 comes from the 5.8.0 custom build
> 0.375000000000004 != 0.375 comes from the 5.6.1 which has been
> installed
>
> I guess the 'long double' is only 80 bits (somehow I cooked up the
> idea that it was 12, but it is longer than the 64 bits that I was
> apparently using before.
>
> Now I have to setup a new installation with all of the libraries that
> were being used (this problem came up in the middle of a lot of
> abstract and complicated functions) and see if it solves the issue
> with the downstream software.

Well, if the difference betwen those 64 bits and the 80 bits gets the
epsilon within the desired range: fine I guess.
But it's not the right solution. What about if someone runs your program on
a different computer and you get contradicting results?

Why not use Math::BigFloat as recommended in the FAQ?

jue

Eric Wilhelm
Guest
Posts: n/a

 11-16-2003
On Sat, 15 Nov 2003 20:59:45 -0600, Jürgen Exner wrote:

> Well, if the difference betwen those 64 bits and the 80 bits gets the
> epsilon within the desired range: fine I guess. But it's not the right
> solution.

You are absolutely right.

Ack! I just realized that the dxf file format was not storing more than 6
digits of accuracy. This was causing the problem in the downsteam app
which is what led me down the snipe hunt of trying to check the precision
inside of the perl library and the further adventures of building a new
perl interpreter.

Many thanks for setting my thinking straight on that problem (though of
course you couldn't have known all of the details from the posting of my
misguided unit test script.) Next time I will remember to question the
storage of the data on disk before pointing the finger at Perl.

using Math::BigFloat due to the number of libraries which would need a
punted the problem by skipping the arc entirely with the thought that it
could come back in another step if need be, but it was bugging me to no
end. Turns out that it all stems from a long-ago-forgotten bit of my own
code with the number "6" in it.

--Eric

Eric J. Roode
Guest
Posts: n/a

 11-16-2003
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Eric Wilhelm <(E-Mail Removed)> wrote in
newsan.2003.11.15.18.26.33.801363.1618@something like.sbcglobalDOTnet:

> On Sat, 15 Nov 2003 18:12:55 -0600, Eric J. Roode wrote:
>
>> Rounding is not a solution? Not even, say, rounding to the 12th
>> decimal place? Surely nothing you're doing can require that level of
>> precision.
>>
>>

> The issue is getting the calculations to resolve between cartesian
> coordinates and trigonometric calculations. If I have an arc that has
> been centered between two points and I calculate its endpoints,
> rounding each coordinate (to any decimal place), they are not matching
> up to the points that were averaged to get to the center point. While
> I am able to sidestep the issue, it is causing problems in downstream
> programs which are expecting higher precision.
>
>> The simple truth is, you cannot reliably use == or != with
>> floating-point numbers. You must do some sort of "is this close
>> enough" comparison on your own. You have to decide what constitutes
>> "close enough".

>
> I have just realized that maybe I have what I was asking for:
>
> 0.375000000000000002 != 0.375 comes from the 5.8.0 custom build
> 0.375000000000004 != 0.375 comes from the 5.6.1 which has been
> installed

If you don't consider two numbers that differ by 10^-12 percent to be
identical, then you are going about things the wrong way.

sub close_enough
{
my (\$a,\$b) = @_;
return abs(\$a - \$b) < 1e-10;
}

Or, if you want to be even more picky:

sub close_enough
{
my (\$a,\$b) = @_;
my \$delta = abs(\$a - \$b);
my \$base = (abs(\$a) + abs(\$b))/2;
return \$delta / \$base < 1e-10;
}

Forget about equality. Installing a more precise library will not solve
your problem; you're dealing with transcendental functions which simply
cannot be represented exactly in any number system.

- --
Eric
\$_ = reverse sort \$ /. r , qw p ekca lre uJ reh
ts p , map \$ _. \$ " , qw e p h tona e and print

-----BEGIN PGP SIGNATURE-----
Version: PGPfreeware 7.0.3 for non-commercial use <http://www.pgp.com>

iQA/AwUBP7cIEGPeouIeTNHoEQLPTACgtHvNl6hB9NphWra/sxYMEN812AUAoLCP
DhL6VJMkhcxiBYqZ7oMEL+Sv
=0A60
-----END PGP SIGNATURE-----