accuracy problem

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

  1. Eric Wilhelm

    Eric Wilhelm Guest

    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
    #1
    1. Advertising

  2. Eric Wilhelm

    Eric Wilhelm Guest

    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
    #2
    1. Advertising

  3. 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
    #3
  4. Eric Wilhelm

    Eric Wilhelm Guest

    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
    #4
  5. -----BEGIN PGP SIGNED MESSAGE-----
    Hash: SHA1

    Eric Wilhelm <> wrote in
    news:p:

    > 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
    #5
  6. Eric Wilhelm

    Eric Wilhelm Guest

    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
    #6
  7. 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
    #7
  8. 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
    #8
  9. Eric Wilhelm

    Eric Wilhelm Guest

    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.

    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
    Eric Wilhelm, Nov 16, 2003
    #9
  10. -----BEGIN PGP SIGNED MESSAGE-----
    Hash: SHA1

    Eric Wilhelm <> wrote in
    news:p:

    > 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
    #10
  11. Eric Wilhelm

    Bob Walton Guest

    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
    #11
  12. -----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
    #12
  13. 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
    #13
    1. Advertising

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

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. =?Utf-8?B?Q2hhcmxlc0E=?=

    Accuracy and CSS

    =?Utf-8?B?Q2hhcmxlc0E=?=, Jan 16, 2006, in forum: ASP .Net
    Replies:
    5
    Views:
    658
  2. Hatzigiannakis Nikos

    Problem with accuracy in DEV C++

    Hatzigiannakis Nikos, Aug 5, 2008, in forum: C Programming
    Replies:
    10
    Views:
    499
    Dik T. Winter
    Aug 6, 2008
  3. Debashish Saha

    accuracy problem in calculation

    Debashish Saha, Nov 8, 2012, in forum: Python
    Replies:
    0
    Views:
    135
    Debashish Saha
    Nov 8, 2012
  4. Chris Angelico

    Re: accuracy problem in calculation

    Chris Angelico, Nov 8, 2012, in forum: Python
    Replies:
    1
    Views:
    151
    Grant Edwards
    Nov 8, 2012
  5. Dave Angel

    Re: accuracy problem in calculation

    Dave Angel, Nov 8, 2012, in forum: Python
    Replies:
    0
    Views:
    156
    Dave Angel
    Nov 8, 2012
Loading...

Share This Page