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