`There is an error here which shouldn't be propagated:	if (fabs(a-b) < 1.0e-38)		..."Machine epsilon" for doubles (namely, the difference between 1.0 andthe next largest number) is on the order of 2e-16.  This is a roughestimate of the granularity of round-off error; and in fact 1.0 / 0.2and 5.0 can't possibly be as close as 1.0e-38, unless they're exactlyequal.There are four epsilon-ish things to be aware of:*	Difference between 0.0 and next float  above: ~= 1.4e-45*	Difference between 0.0 and next double above: ~= 4.9e-324*	Difference between 1.0 and next float  above: ~= 1.2e-7*	Difference between 1.0 and next double above: ~= 2.2e-16The first two are more useful for things like detecting underflow; thelast two (some numerical folks suggest using their square roots) aremore useful for implementing an "approximately equals".----------------------------------------------------------------The poster was incorrect in expecting 1.0 / 0.2 to be exactly equal toanything, as was explained to him.  But the problem doesn't have to dowith whether a number is transcendental, or irrational, or rational:the number must be rational *and* must have a mantissa whosedenominator is a power of two *and* that power of two must be less thanor equal to 23 (for single) or 52 (for double).  And of course 1/5 is2^-3 * 8/5, of which the mantissa has denominator 5, which isn't a powerof two.So we all should know not to expect floating-point numbers to beexactly equal to anything; that's been established.  However, anothermore basic question was not answered; curiosity (if nothing else)demands an answer.  Namely, it's OK to say we can't expect 1.0/0.2 ==5.0.  But why is the result of (what is apparently) the samecomputation *sometimes* the same, and *sometimes* different? That's thequestion.And I think it's fair for the poster to want to know why.If you disassemble the sample program, you'll see that withoutoptimization, 1.0 is divided by 0.2 at *run* time, and compared with5.0; with optimization, the division is done, and the "<" and"==" comparisons are done, at *compile* time.  OK, but: If we're notcross-compiling (most people don't), then the compiler creating a.outis running on perhaps the same box as a.out is!  Why does gcc, foldingthe constant in the optimized a.out, get a different answer for 1.0/0.2than the unoptimized a.out gets for 1.0/0.2?Not only that, without optimization:	if (1/h < 5.0)		...gives a different answer inside a.out than	x = 1/h;	if (x < 5.0)		...The key is that Pentiums (Pentia?) have 80-bit floating-point numbersin the FPU.  Without optimization, at compile time, gcc represents 5.0as 0x4014000000000000.  0.2 is 0x3fc999999999999a.  These are both64-bit doubles -- 1 sign bit, 11 exponent bits, & 52 explicit mantissabits (and 1 implicit leading mantissa bit, not stored in memory.)In the case "if (1/h < 5.0)", at run time, 1.0 is loaded into the FPUusing fld1; then "fdivl {address of 0.2 in memory}".  The result is the*80-bit* number 0x40019ffffffffffffd80.  The 64-bit number 5.0(0x4014000000000000) is loaded into the FPU to become the 80-bit number0x4001a000000000000000.  Then, these two 80-bit numbers are compared inthe FPU; they're of course not the same.What's different in the case "x = 1/h; if (x < 5.0) ..." is that both80-bit numbers are stored from the FPU to memory as 64-bit (roundingoff the mantissa bits which differ), at which point they're both0x4014000000000000, then loaded *back* into the FPU where they'reboth 0x4001a000000000000000.This isn't an FPU bug, by any stretch of the imagination, nor is it acompiler bug.  But it's a subtle difference between the Pentium's FPUand other FPUs, of which it may occasionally be useful to be aware.-----Original Message-----From: Richard B. Johnson [mailto:root@chaos.analogic.com]Sent: Thursday, April 25, 2002 7:23 AMTo: rpmCc: Jesse Pollard; Nikita@Namesys.COM; Andrey Ulanov;linux-kernel@vger.kernel.orgSubject: Re: FPU, i386On Thu, 25 Apr 2002, rpm wrote:> On Wednesday 17 April 2002 08:10 pm, Jesse Pollard wrote:> > ---------  Received message begins Here  ---------> >> > > if (int(1/h * 100) == int(5.0 * 100))> >> > will give a "proper" result within two decimal places. This is still> > limited since there are irrational numbers within that range that COULD> > still come out with a wrong answer, but is much less likely to occur.> >> > Exact match of floating point is not possible - 1/h is eleveated to a> > float.> >> > If your 1/h was actually num/h, and num computed by summing .01 100times> > I suspect the result would also be "wrong".> >> > why is exact match of floating point not possible ?Because many (read most) numbers are not exactly representablein floating-point. The purpose of floating-point it to representreal numbers with a large dynamic range. The trade-off is thatfew such internal representations are exact.As a simple example, 0.33333333333.....  can't be represented exactlyeven with paper-and-pencil. However, as the ratio of two integersit can be represented exactly, i.e., 1/3 . Both 1 and 3 mustbe integers to represent this ratio exactly.All real numbers (except trancendentials) can represented exactlyas the ratio of two integers but floating-point uses only onevalue, not two integers, to represent the value. So, an exactrepresentation of a real number, when using a single variablein a general-purpose way, is, for all practical purposes, notpossible. Instead, we get very close.When it comes to '==' close is not equal. There are macros in<math.h> that can be used for most floating-point logic. Youshould check them out. If we wanted to check for '==' we reallyneed to do something like this:    double a, b;    some_loop() {       if(fabs(a-b) < 1.0e-38)           break;     }Where we get the absolute value of the difference between twoFP variables and compare against some very small number.To use the math macros, the comparison should be something like:        if (isless(fabs(a-b), 1.0e-38))             break;Cheers,Dick JohnsonPenguin : Linux version 2.4.18 on an i686 machine (797.90 BogoMips).                 Windows-2000/Professional isn't.-To unsubscribe from this list: send the line "unsubscribe linux-kernel" inthe body of a message to majordomo@vger.kernel.orgMore majordomo info at  http://vger.kernel.org/majordomo-info.htmlPlease read the FAQ at  http://www.tux.org/lkml/-To unsubscribe from this list: send the line "unsubscribe linux-kernel" inthe body of a message to majordomo@vger.kernel.orgMore majordomo info at  http://vger.kernel.org/majordomo-info.htmlPlease read the FAQ at  http://www.tux.org/lkml/`