[Back to DATATYPE SWAG index]  [Back to Main SWAG index]  [Original]

{
> Still, I expect as a programmer for the CPU to have the right answer 100%
> of the time! My own errors are bad enough. :) It's still a bug and until
> it's 100% accurate and no less, it's not desirable to me and I would think
> others shall agree. Anything less is unreliable and unreliable doesn't
> cut it.

This is getting off-topic for the echo.  Please, if what you have to say has
nothing to do with programming in Pascal, don't say it here.

Here's something just posted to comp.lang.pascal on Usenet:

 From: Terje.Mathisen@hda.hydro.com (Terje Mathisen)
 Subject: FDIV fix for BP (fast & exact)
 Date: 1 Dec 1994 11:07:51 GMT

I have been working with Cleve Moler (from MathWorks) and Tim Coe about a sw
workaround for the Pentium FDIV bug.  Here is a BP version of the algorithm:

The FDIV_FIX unit defines a single function:

function fdiv(x: extended; y: extended): extended;

which will use about 25% more cycles than a single fdiv call, but always return
exact results for all double precision numbers, including the special cases of
Nan, Inf and denormal.
It also handles all extended precision numbers, giving a maximum error of a
single bit in the last position (1ulp).  This error can only be introduced if
the FDIV operations fails.
During the unit startup code, I test the FDIV instruction, and if if fails,
initialize a 16-byte table of critical nibble values.

If FDIV passes, the table is left empty, and no fixups will be done on
divisions.
=============== FDIV_FIX.PAS ===========================
}
{$N+,E-,G+}
Unit FDIV_FIX;

Interface

function fdiv(x, y : Extended): Extended;

const
  fdiv_risc_table : array [0..15] of byte = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  fdiv_scale : single = 0.9375;
  one_shl_63_l : Longint = $5f000000;
var
  one_shl_63 : single absolute one_shl_63_l;

Implementation

function fdiv(x, y : Extended): Extended;
Assembler;
asm
  fld [x]
@restart:
  mov bx,word ptr [y+6]
  fld [y]
  add bx,bx
   jnc @denormal
{$IFOPT G+}
  shr bx,4
{$ELSE}
  mov cl,4
  shr bx,cl
{$ENDIF}
  cmp bl,255
   jne @ok
{$IFOPT G+}
  shr bx,8
{$ELSE}
  mov bl,bh
  and bx,255
{$ENDIF}
  cmp byte ptr fdiv_risc_table[bx],bh
   jz @ok
  fld [fdiv_scale]
  fmul st(2),st
  fmulp st(1),st
   jmp @ok
@denormal:
  or bx,word ptr [y]
  or bx,word ptr [y+2]
  or bx,word ptr [y+4]
   jz @zero
  fld [one_shl_63]
  fmul st(2),st
  fmulp st(1),st
  fstp [y]
   jmp @restart
@zero:
@ok:
  fdivp st(1),st
end;

const
  a_bug : single = 4195835.0;
  b_bug : single = 3145727.0;

Procedure fdiv_init;
var
  r : double;
  i : Integer;
begin
  r := a_bug / b_bug;
  if a_bug - r * b_bug > 1.0 then begin
    i := 1;
    repeat
      fdiv_risc_table[i] := i;
      Inc(i,3);
    until i >= 16;
  end;
end;

begin
  fdiv_init;
end.

[Back to DATATYPE SWAG index]  [Back to Main SWAG index]  [Original]