accuracy problem

E

Eric Wilhelm

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
 
E

Eric Wilhelm

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

Jürgen Exner

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
 
E

Eric Wilhelm

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
 
E

Eric J. Roode

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

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

Eric Wilhelm

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
 
J

Jürgen Exner

Eric said:
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

Jürgen Exner

Eric said:
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
 
E

Eric Wilhelm

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.

Yes, I've read perldoc -q 999 several times and had already ruled out
using Math::BigFloat due to the number of libraries which would need a
complete rewrite to carry the accuracy through. I had actually already
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
 
E

Eric J. Roode

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

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.


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

Bob Walton

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

Eric J. Roode

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

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

Tassilo v. Parseval

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
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

Forum statistics

Threads
473,744
Messages
2,569,483
Members
44,902
Latest member
Elena68X5

Latest Threads

Top