C As this below tiny program from Tim Coe, coe@vitsemi.com demonstrates
C the Pentium procesor sometimes returns wrong results when using real
C division. EPS-dependent Fortran version. This routine modified to
C use Cleve Moler's code to circumvent the Pentium Processor Divide
C bug.
C H. D. Knoble, hdk@psuvm.psu.edu, 11/27/94
C Penn State Center for Academic Computing
DOUBLE PRECISION X,Y,Z,FDIV
REAL EPS
C
C---Define SP EPS for the IBM PC. EPS is the smallest real number on
C this architecture such that 1+EPS>1 and 1-EPS<1.
EPS=5.9604640E-08
X = 4195835.0D0
Y = 3145727.0D0
C---Z Theoretically should be a computational (fuzzy) Zero.
C By substituting FDIV(X,Y) for (X / Y) below, the Pentium
C division bug can be avoided.
C Z = X - (X / Y) * Y
Z = X - FDIV(X,Y) * Y
WRITE(*,*) 'Residual Z=',Z
IF (DABS(Z).LT.EPS) THEN
WRITE(*,*) 'Since Z is essentially Zero, then this Pentium'
WRITE(*,*) 'Processor''s Divide bug has been circumvented.'
ELSE
WRITE(*,*) 'If Z is substantially not Zero, then this Pentium'
WRITE(*,*) 'Pentium Processor''s Divide bug has NOT been'
WRITE(*,*) 'circumvented.'
ENDIF
STOP
END
DOUBLE PRECISION FUNCTION FDIV(U,V)
C===Cleve Moler's Fix for Pentium Bug and comments.
C The idea is to use the hardware to divide u by v and produce the
C candidate quotient, z. This will be correct in the vast majority
C of cases, but it is necessary to check. The test uses the residual,
C r = u - v*z. Normally, r will be no larger than roundoff error in u.
C The constant EPS involved in the test is the distance from 1.0
C to the next larger floating point number. For double precision,
C EPS is 2**(-52). At first, we forgot to include the constant
C REALMIN in the test, but this term is required to deal correctly
C with denormal floating point numbers.
C
C If the residual fails the test, it indicates that the hardware bug
C has been encountered. When this occurs, we simply rescale the
C numerator and denominator by a factor of 3/4 and repeat the process.
C The scaling scrambles the bit patterns in u and v so that the second
C division almost certainly gives a satisfactory result.
C
C The factor 3/4 is important. It is the "simplest" factor which
C alters the bit patterns in the fractions of u and v. If the last
C bit in the fraction is zero, the scaling does not introduce any
C roundoff error. Furthermore, since 3/4 is less than 1, the
C scaling cannot overflow. And, if the operands are so small that
C scaling by 3/4 would underflow, the original division is done
C correctly and the scaling is not needed.
C
C The FDIV instruction can be used to correct itself! Here's the
C C-code converted to Fortran.
C
DOUBLE PRECISION X,Y,R,EPS,RMIN,RHO,U,V
DATA EPS /2.2204460492503131E-16/
DATA RMIN /2.225073858507202D-308/
DATA RHO /0.75/
X=U
Y=V
FDIV = X / Y
R = X - Y*FDIV
C---DO while (fabs(r) > EPS*fabs(x)+RMIN
IF (ABS(R).LE. EPS*ABS(X)+RMIN) RETURN
DO 1 I=1,9
X = RHO*X
Y = RHO*Y
FDIV = X/Y
R = X - Y*FDIV
IF (ABS(R).LE. EPS*ABS(X)+RMIN) RETURN
1 CONTINUE
STOP 9
END