# accuracy problem

Discussion in 'Perl Misc' started by Eric Wilhelm, Nov 15, 2003.

1. ### Eric WilhelmGuest

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, Nov 15, 2003

2. ### Eric WilhelmGuest

Re: accuracy (er, precision) problem

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

Eric Wilhelm, Nov 15, 2003

3. ### Jürgen ExnerGuest

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

Jürgen Exner, Nov 15, 2003
4. ### Eric WilhelmGuest

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 Wilhelm, Nov 16, 2003
5. ### Eric J. RoodeGuest

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Eric Wilhelm <> wrote in
news:

> 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+bmr2PkVwAoJAY
6LcSvscd+oJBSKee0cX71BNA
=U1By
-----END PGP SIGNATURE-----

Eric J. Roode, Nov 16, 2003
6. ### Eric WilhelmGuest

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 128), 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

Eric Wilhelm, Nov 16, 2003
7. ### Jürgen ExnerGuest

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, Nov 16, 2003
8. ### Jürgen ExnerGuest

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 128), 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

Jürgen Exner, Nov 16, 2003
9. ### Eric WilhelmGuest

Re: ack! problem

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 Wilhelm, Nov 16, 2003
10. ### Eric J. RoodeGuest

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Eric Wilhelm <> wrote in
news:

> 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-----

Eric J. Roode, Nov 16, 2003
11. ### Bob WaltonGuest

Purl Gurl wrote:

....

> Quite ironic in this perceived high tech era, we cannot
> design a computer which can do math anywhere near as well
> as the human mind or an official Number 2 pencil.
>

Well, actually, it is the computer software. Try Maxima -- it has a
wonderful "bfloat" package which will compute most algebraic and trig
expressions to a specified precision, which can easily be 1000 digits or
more. Example:

(C9) fpprec:100;
(D9) 100
(C10) 4*atan2(1,1),bfloat;
(D10)
3.1415926535897932384626433832795028841971693993751058209749445923078164#

06286208998628034825342117068B0
(C11)

It does what you might hope Math::BigFloat would do with trig functions.

Its free at: http://maxima.sourceforge.net/

> Purl Gurl

--
Bob Walton
Email: http://bwalton.com/cgi-bin/emailbob.pl

Bob Walton, Nov 16, 2003
12. ### Eric J. RoodeGuest

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Purl Gurl <> wrote in
news::

> I will have a look at this software. This will make a nice addition
> for my navigation system in my spacecraft. I found during my yearly
> visits to my homeworld, Alturus 5 out in the Gamma Sector, my usual
> navigated landing spot is a few miles short of my hometown. I really
> hate having to walk those few miles through a methane swamp. Phew!
>
> Personal computers today, in reality, do a pretty good job on
> math accuracy. I cannot think of circumstances requiring more
> accuracy than offered, save for navigating many thousands of
> light-year distances.

I don't think even that hypothetical application requires such accuracy.
I believe Isaac Asimov calculated that if you were to construct a circle
whose diameter was the diameter of the observable universe, and if you
were to compute its circumference using a value of pi that was "only"
accurate to 35 decimal places, your computation would be off by less than
the diameter of a hydrogen atom.

- --
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/AwUBP7gaumPeouIeTNHoEQIPigCffM90oTi1UNUjIqtxjsq1gMRHtxkAoKzt
1BgP+TtiDFFJgtk4fvOdMnlN
=XWrq
-----END PGP SIGNATURE-----

Eric J. Roode, Nov 17, 2003
13. ### Tassilo v. ParsevalGuest

Also sprach Purl Gurl:

> Personal computers today, in reality, do a pretty good job on
> math accuracy. I cannot think of circumstances requiring more
> accuracy than offered, save for navigating many thousands of
> light-year distances.

Actually (and in the fields of scientific computing) accuracy is a huge
issue. You often have to deal with problems such as inverting a
10000x10000 matrix or finding the eigenvalues of such a monstrosity. To
make that feasible at all, you throw many CPUs at it and therefore have
to distribute the load evenly across processors using techniques such
as openMP (will split the program for you more or less intelligently) or
MPI (where you have to come up yourself with parallalized algorithms).

When doing this in an unclever (but obvious) way, the computed result
wont even be close to the actual result, and with respect to this it's
quite irrelevant whether 'doubles' have a width of 8 bytes or 32 bytes.
It's still the programmer who has to do it right and there's no hope
that a machine could help him.

Tassilo
--
\$_=q#",}])!JAPH!qq(tsuJ[{@"tnirp}3..0}_\$;//::niam/s~=)]3[))_\$-3(rellac(=_\$({
pam{rekcahbus})(rekcah{lrePbus})(lreP{rehtonabus})!JAPH!qq(rehtona{tsuJbus#;
\$_=reverse,s+(?<=sub).+q#q!'"qq.\t\$&."'!#+sexisexiixesixeseg;y~\n~~dddd;eval

Tassilo v. Parseval, Nov 17, 2003