From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me X-Spam-Level: X-Spam-Status: No, score=-1.3 required=5.0 tests=BAYES_00,INVALID_MSGID autolearn=no autolearn_force=no version=3.4.4 X-Google-Language: ENGLISH,ASCII-7-bit X-Google-Thread: 107079,c7637cfdf68e766 X-Google-Attributes: gid107079,public X-Google-Thread: 103376,c7637cfdf68e766 X-Google-Attributes: gid103376,public X-Google-Thread: f43e6,c7637cfdf68e766 X-Google-Attributes: gidf43e6,public X-Google-Thread: f8362,c7637cfdf68e766 X-Google-Attributes: gidf8362,public X-Google-Thread: 109d8a,c7637cfdf68e766 X-Google-Attributes: gid109d8a,public From: pmontgom@cwi.nl (Peter L. Montgomery) Subject: Re: floating point comparison Date: 1997/08/02 Message-ID: #1/1 X-Deja-AN: 261523236 Sender: news@cwi.nl (The Daily Dross) References: <01bc9dca$762b6000$e1e989cc@gort> <33E261F7.2EFB@mailgw.sanders.lockheed.com> Organization: CWI, Amsterdam Newsgroups: comp.lang.ada,sci.math.num-analysis,comp.software-eng,comp.theory,sci.math Date: 1997-08-02T00:00:00+00:00 List-Id: In article <33E261F7.2EFB@mailgw.sanders.lockheed.com> user writes: >Bob Binder wrote: >> I haven't followed the IEEE fp standard, but I can't imagine how >> (without resorting to some kind of BCD representation) any floating >> point scheme using a non-decimal radix gets around the fundamental fact >> that there are non-terminal fractions in this radix which are not >> non-terminal in a decimal radix. >> Try this: >> x, y float; >> x = 1.0; >> y = (x/10.0)*10.0; >> if x == y { >> // not likely } >> else { >> // more likely >> }; >Although I am reading this in comp.lang.ada, the suggestion to try a >simple example like this in C piqued my interest. > >On a Sun SPARC-5, with the built in compiler supplied by Sun, this >code (fixed for syntax) found (x == y) to be true. Using 3.0 and >7.0 (neither of which have any chance of yielding an exact >representation >after the division) also found (x == y) to be true. You lucked out due to IEEE rounding to nearest. Each nonzero floating point value is approximated by +- 1 times an integer mantissa M in [2^52, 2^53) times a power 2^e, where e is an integer. When approximating 1/10 = 2^(-56) * (2^56/10), we want M to be the integer nearest to 2^56/10. This is (2^56 + 4)/10. The computed value of 1.0/10.0 becomes ((2^56 + 4)/10) * 2^(-56). Now multiply back by 10.0. The exact product would be (2^56 + 4)*2^(-56), but the mantissa is required to be in [2^52, 2^53). The result is scaled to (2^52 + 0.25)*2^(-52). Now a rounding gives 2^52 * 2^(-52), or 1.0. The final result happens to match the original x, even though intermediate results were inexact. On the other hand, try (1.0/49.0)*49.0. The quotient 1/49 is approximated by ((2^58 - 23)/49) * 2^(-58). Now multiply back by 49. The normalized mantissa 2^53 - 23/32 rounds to 2^53 - 1, and the product is approximated by (2^53 - 1)*2^(-53). [Yes, I ignored denormals and other special operands above.] -- Peter L. Montgomery pmontgom@cwi.nl San Rafael, California A mathematician whose age has doubled since he last drove an automobile.