Very Computer

a-((a/b)*b) = ?

It should be zero, of course.  The result of pentium calculation
is 256!!!
So, I checked that (I use fortran) and resutl is 256 on my
Pentium (90 MHz).  Too bad. But not so bad. I compiled the same
program with O2 optimization (f77 -O2 ...) Result this time is 0!
(as it should be).  
The results are independent of of single/double precision.

Anybody have any idea how is optimization doing this?

Thanks for bringing this up.

--
Paul Bash  

Quote:>Hot subject about Pentiums these days is recently discovered
>error in floating point unit (article in NYT, Thursday, Nov. 24).
>Example published there is as follows:
>a=4195835
>b=3145727

>a-((a/b)*b) = ?

>It should be zero, of course.  The result of pentium calculation
>is 256!!!
>So, I checked that (I use fortran) and resutl is 256 on my
>Pentium (90 MHz).  Too bad. But not so bad. I compiled the same
>program with O2 optimization (f77 -O2 ...) Result this time is 0!
>(as it should be).  
>The results are independent of of single/double precision.

>Anybody have any idea how is optimization doing this?

If you could figure out the correct result without doing any calculation,
a good optimizer could do the same thing, too.  If you have doubts,
look at the assembly code generated by the compiler:

    ues4:~/tmp 141> cat p5bug.c
    #include <stdio.h>
    main()
    {
        double a = 4195835, b = 3145727;
        printf("%f\n", a - (a / b) * b);
        return 0;
    }
    ues4:~/tmp 142> gcc -S p5bug.c
    ues4:~/tmp 143> cat p5bug.s
            .file   "p5bug.c"
    gcc2_compiled.:
    ___gnu_compiled_c:
    .text
    LC0:
            .ascii "%f\12\0"
            .align 4
    .globl _main
    _main:
            pushl %ebp
            movl %esp,%ebp
            subl $16,%esp
            call ___main
            movl $-1073741824,-8(%ebp)
            movl $1095762302,-4(%ebp)
            movl $-2147483648,-16(%ebp)
            movl $1095237631,-12(%ebp)
            fldl -8(%ebp)
            fdivl -16(%ebp)
            fmull -16(%ebp)
            fldl -8(%ebp)
            fsubp %st,%st(1)
            subl $8,%esp
            fstpl (%esp)
            pushl $LC0
            call _printf
            addl $12,%esp
            xorl %eax,%eax
            jmp L1
            .align 4,0x90
    L1:
            movl %ebp,%esp
            popl %ebp
            ret

Even without being an assembly language wizard, you can easily identify
the floating point instructions generated by the compiler.

And now, let's activate the first level of optimizations in gcc:

    ues4:~/tmp 144> gcc -S -O p5bug.c
    ues4:~/tmp 145> cat p5bug.s
            .file   "p5bug.c"
    gcc2_compiled.:
    ___gnu_compiled_c:
    .text
    LC0:
            .ascii "%f\12\0"
            .align 4
    .globl _main
    _main:
            pushl %ebp
            movl %esp,%ebp
            call ___main
            pushl $0
            pushl $0
            pushl $LC0
            call _printf
            xorl %eax,%eax
            movl %ebp,%esp
            popl %ebp
            ret

Poof!  No floating point instructions in the optimized version.  The
compiler has generated the equivalent of printf("%f\n", 0.0).  

Who said the Pentium chip is buggy?!  :-)

Dan
--
Dan Pop
CERN, CN Division

Mail:  CERN - PPE, Bat. 31 R-004, CH-1211 Geneve 23, Switzerland

: Who said the Pentium chip is buggy?!  :-)

The obvious fix!  Don't do any divisions!  :->

x-(x/y)*y

with 0. I think someone already explained this in comp.sys.intel.

Lech Borkowski

a-((a/b)*b) = a-(a*b/b) = a - a = 0

This is done without performing any divisions - it just uses the
algebraic laws.

Morten Joehnk

--
From the dictionary: recursive, adj: see recursive.

Consider the following (back to elementary scool stuff):

1) Assume a number N exists that will always yield the correct result on
   a pentium when used as a nominator or denominator in a division with any
   other number.

2) Then replace the expression
       c = a/b
   by
       c = (a*N) / (b*N)
   which of course is the same as
       c = (a/N) * (N/b)

This should give the correct results (excluding some small rounding errors)
unless the pentium can't multiply correctly either.

I know that replacing a division by two divisions and a multiplication may slow
things down a bit, but at least the result will be correct :->

Any suggestions for the value of N ?

/fch

: Any suggestions for the value of N ?

"Zero" comes to mind.   :-)

: : Any suggestions for the value of N ?

: "Zero" comes to mind.   :-)

How about Pi? :-)(-:

I tried using a C version of the test program with and
without optimization and it turns out that for any level
of optimization (-O1 to -O6) on this 90Mhz Pentium the
result is the correct 0.0000, without the -O option the
result is 256.000. This is of cource runing gcc 2.5.8 under
linux.

I also tried using my old TC version 2.0 optimized for speed,
but the result was wrong all the time. (I'm waiting for my
co-worker to try it under Visual C to se if this has any thing
to do with the OS or compiler). Maybe this is yet another
reason to make the switch from DOS to Linux?

(Pentium Power under DOS is kind of silly anyway.)

Hallvard P.

There may be plenty of reasons that one might wish to switch to Linux from
DOS, but the particular reason given above may not necessarily be one of
them.  The reason is that the optimization by gcc only hides the problem
(which is in the processor, not in gcc).  The problem will still occur in
instances in which gcc cannot optimize the division out of the emitted
object code.  

On the other hand, if the division were done correctly in the first place,
you wouldn't be able to tell the difference between the optimized and
unoptimized code by the numerical results, anyway.

I have heard that gcc may be modified soon to emit code that works around
the Pentium bug.  Most compilers may have to incorporate such a workaround
in the near future.  

Therefore, what may, in fact, be good reasons to switch to Linux are the
facts that gcc is likely to be one of the first compiler to incorporate a
work-around, and that most Linux code is distributed as source code that
can be recompiled to incorporate the work-around without waiting for a bug
fix from each publisher.

: > "Zero" comes to mind.   :-)

: Come on now, you can do better than that!
: Zero is not really a "viable" solution (sorry, I couldn't resist). :-)

But you've got to admit that it very neatly avoids the problem of "small
rounding errors!"

--------
"Slimy, yet satisfying!" - Simba, the lion, after eating some bugs

/fch

[...]

Quote:>Consider the following (back to elementary scool stuff):
>1) Assume a number N exists that will always yield the correct result on
>   a pentium when used as a nominator or denominator in a division with any
>   other number.
>2) Then replace the expression
>       c = a/b
>   by
>       c = (a*N) / (b*N)
>   which of course is the same as
>       c = (a/N) * (N/b)

Quote:>This should give the correct results (excluding some small rounding errors)
>unless the pentium can't multiply correctly either.
>I know that replacing a division by two divisions and a multiplication may slow
>things down a bit, but at least the result will be correct :->

Quote:>Any suggestions for the value of N ?

Quote:>/fch

All credit to MathWorks - this is a pretty quick response time to a bug
whih isn't even their fault...

ray

-> Date: Thur Nov 24 22:00:20 EST 1994
-> Subject: Soft/Hardware Workaround for the FDIV bug
-> Newsgroups: comp.sys.intel
->
->         SOFTWARE/HARDWARE WORKAROUND FOR THE FDIV BUG
->
-> At the MathWorks, we have decided to issue a new release of
-> MATLAB which is "Pentium aware".  It incorporates a workaround
-> for the floating point division bug which restores full accuracy
-> without a serious degradation of efficiency.  And, we have
-> decided to post a description of our technique to the Internet
-> so that other people, including other commercial software
-> developers, can make use of it.
->
-> It is not easy to modify a large package like MATLAB to include
-> this change.  We would much prefer to have the compiler or the
-> operating system do the job for us, but this is not yet an
-> option.  The kernel of MATLAB is written in C.  We have replaced
-> several dozen instances of a '/' denoting a floating point
-> division by a function call.  We are now in the process of
-> confirming that we didn't introduce any source code errors
-> while doing this.
->
-> Our initial intention was to use the Pentium hardware FDIV
-> instruction to compute a candidate quotient, check its accuracy,
-> and then, if necessary, employ a standard Newton iterative
-> refinement technique to eliminate any inaccuracy produced
-> by the hardware.  But then we realized that an approach unique
-> to this situation was possible, and more effective.  The FDIV
-> instruction can be used to correct itself!  Here's the code:
->
->      #include <math.h>
->      #define EPS      2.2204460492503131e-16
->      #define REALMIN  2.2250738585072014e-308
->      #define RHO      0.75
->      
->      double fdiv(double x, double y)
->      {
->         int ok;
->         double r,z;
->         z = x/y;
->         r = x - y*z;
->         while (fabs(r) > EPS*fabs(x)+REALMIN) {
->            x = RHO*x;
->            y = RHO*y;
->            z = x/y;
->            r = x - y*z;
->         }
->         return(z);
->      }
->
-> The idea is to use the hardware to divide x by y and produce the
-> candidate quotient, z.  This will be correct in the vast majority
-> of cases, but it is necessary to check.  The test uses the residual,
-> r = x - y*z.  Normally, r will be no larger than roundoff error in x.
-> The constant EPS involved in the test is the distance from 1.0
-> to the next larger floating point number.  For double precision,
-> EPS is 2^(-52).  At first, we forgot to include the constant
-> REALMIN in the test, but this term is required to deal correctly
-> with denormal floating point numbers.
->
-> If the residual fails the test, it indicates that the hardware bug
-> has been encountered.  When this occurs, we simply rescale the
-> numerator and denominator by a factor of 3/4 and repeat the process.
-> The scaling scrambles the bit patterns in x and y so that the second
-> division almost certainly gives a satisfactory result.
->
-> The factor 3/4 is important.  It is the "simplest" factor which
-> alters the bit patterns in the fractions of x and y.  If the last
-> bit in the fraction is zero, the scaling does not introduce any
-> roundoff error.  Furthermore, since 3/4 is less than 1, the
-> scaling cannot overflow.  And, if the operands are so small that
-> scaling by 3/4 would underflow, the original division is done
-> correctly and the scaling is not needed.
->
-> How much does all this cost in execution speed?   We have not yet
-> done any serious timings outside the MATLAB context.  Since the
-> occurrence of the inaccurate results is so rare, any time spent
-> inside the while loop can be ignored.  The FDIV instruction itself
-> takes 30 or more cycles.  To this we must add the time for
-> computing the residual, the time for doing the test and, perhaps
-> the most significant, the time required by the function call.  We
-> guess that this might double the time, but the actual factor will
-> depend upon the surrounding environment.
->
-> If anybody actually does incorporate this approach in their
-> software, we'd like to hear about it.  In particular, we'd
-> appreciate hearing about any timing measurements from other
-> large packages.
->
-> I will be posting an article to the comp.soft-sys.matlab newsgroup
-> describing the Pentium aware MATLAB, and its availability, in
-> more detail.
->
-> Finally, although I have not yet met any of them in person,
-> I would like acknowledge and thank Thomas Nicely, Tim Coe and
-> Mike Carlton for their contributions to this enterprise.
->
->   -- Cleve Moler
->   Chairman and Chief Scientist
->   The MathWorks, Inc.

->

There's a standard opt called constant folding. What it does is
propogate constants through the program, and pre-execute any
expressions that are entirely constant. So that expression is just
being optimized away by the compiler.

The reason that it doesn't result in a bug is that the compiler
writers didn't waste time on figuring out whether or not there's an
FPU available - they just do floats by emulation (or whatever). So it
doesn't get run on the buggy pentium fpu.

        <MC>

--