{$ifdef LINUX}
{$include ../defines.pas}
{$else}
{$include ..\defines.pas}
{$endif}
unit fpu_system;

interface

uses config, fpu_types, cpu;

procedure add_two_Xsig (dest: PXsig; const x2: PXsig; exp: ps32);
procedure negate_Xsig (x: PXsig);
function XSIG_LL (x: Xsig): pu64;
function LL_MSW (x: Pointer): pu64;
procedure add_Xsig_Xsig (dest: PXsig; const x2: PXsig);
function mul_32_32 (const arg1: u32; const arg2: u32): u32;
function access_limit: pu_char;
{$include status_w.pas}

const
  ERR_MARGIN = $20;

 { Keep the order TAG_Valid, TAG_Zero, TW_Denormal }
 { The following fold to 2 (Special; in the Tag Word }

const
  FWAIT_OPCODE     = $9b;
  OP_SIZE_PREFIX   = $66;
  ADDR_SIZE_PREFIX = $67;
  PREFIX_CS        = $2e;
  PREFIX_DS        = $3e;
  PREFIX_ES        = $26;
  PREFIX_SS        = $36;
  PREFIX_FS        = $64;
  PREFIX_GS        = $65;
  PREFIX_REPE      = $f3;
  PREFIX_REPNE     = $f2;
  PREFIX_LOCK      = $f0;
  PREFIX_CS_       = 1;
  PREFIX_DS_       = 2;
  PREFIX_ES_       = 3;
  PREFIX_FS_       = 4;
  PREFIX_GS_       = 5;
  PREFIX_SS_       = 6;
  PREFIX_DEFAULT   = 7;

type
  FUNC     = procedure;
  FUNC_ST0 = procedure (st0_ptr: PFPU_REG; st0_tag: u_char);

const
  PROTECTED_ = 4;

const
  SIXTEEN = 1;         // We rely upon this being 1 (true)

const
  VM86 = SIXTEEN;

const
  PM16 = (SIXTEEN or PROTECTED_);

const
  SEG32 = PROTECTED_;

function FPU_normalize_nuo (x: PFPU_REG; bias: integer): intg;
function FPU_u_sub (arg1: PFPU_REG; arg2: PFPU_REG; dest: PFPU_REG;
  control_w: u16; sign: u_char; expa: s32; expb: s32): intg;
function FPU_u_mul (a: PFPU_REG; b: PFPU_REG; c: PFPU_REG; cw: u16;
  sign: u_char; expon: s32): intg;
function FPU_u_div (a: PFPU_REG; b: PFPU_REG; dest: PFPU_REG;
  control_w: u16; sign: u_char): intg;
function FPU_u_add (arg1: PFPU_REG; arg2: PFPU_REG; answ: PFPU_REG;
  control_w: u16; sign: u_char; expa: s32; expb: s32): intg;
function wm_sqrt (n: PFPU_REG; dummy1: intg; dummy2: intg; control_w: u16;
  sign: u_char): intg;
function FPU_shrx (arg1: pointer; arg2: u32): intg;
function FPU_shrxs (arg1: pointer; arg2: u32): intg;
function FPU_div_small (x: pu64; y: u32): intg;

function no_ip_update: u_char;

function exponent (x: Pfpu_reg): longword;
function top: u32;
function register_base: pu_char;
procedure poppop;
procedure FPU_illegal;
procedure pop_0 ();
function status_word: intg;
procedure push;

var
  data_sizes_16:  array[0..32] of u_char;
  I387:           i387_t;
  fpu_cpu_ptr:    PBX_CPU_C;
  fpu_iptr:       PBxInstruction_tag;
  CONST_1:        FPU_REG;
  CONST_2:        FPU_REG;
  CONST_HALF:     FPU_REG;
  CONST_L2T:      FPU_REG;
  CONST_L2E:      FPU_REG;
  CONST_PI:       FPU_REG;
  CONST_PI2:      FPU_REG;
  CONST_PI2extra: FPU_REG;
  CONST_PI4:      FPU_REG;
  CONST_LG2:      FPU_REG;
  CONST_LN2:      FPU_REG;
  CONST_Z:        FPU_REG;
  CONST_PINF:     FPU_REG;
  CONST_INF:      FPU_REG;
  CONST_MINF:     FPU_REG;
  CONST_QNaN:     FPU_REG;

  fp_fpu: integer;

procedure reg_copy (x: PFPU_REG; y: PFPU_REG);
function fpu_register (index: integer): PFPU_REG;
procedure setcc (cc: u32);

implementation

uses Service_fpu, Math, reg_round, errors, fpu_tags, wmfpuemu_Glue;

procedure reg_copy (x: PFPU_REG; y: PFPU_REG);
  var
    src, dest: PChar;
    I: integer;
  begin
    y^.exp := x^.exp;
  {src:=pchar(x);
  dest:=pchar(y);
  if Cardinal(Dest) > Cardinal(Src) then
    for I := sizeof(int64)-1 downto 0 do
      dest[I] := src[I]
  else
    for I := 0 to sizeof(int64)-1 do
      dest[I] := src[I];}
    significand(y)^ := significand(x)^;
    //PU64(@y^.SIGH)^:=PU64(@X^.SIGH)^;
  end;

function FPU_normalize_nuo (x: PFPU_REG; bias: integer): intg;
  begin

    if ((x^.sigh and $8000) = 0) then
      begin
      if (x^.sigh = 0) then
        begin
        if (x^.sigl = 0) then
          begin
          x^.exp := EXP_UNDER;
          Result := TAG_Zero;
          Exit;
          end;
        x^.sigh := x^.sigl;
        x^.sigl := 0;
        x^.exp  := x^.exp - 32;
        end;
      while ((x^.sigh and $8000) = 0) do
        begin
        x^.sigh := x^.sigh shl 1;
        if (x^.sigl and $8000) <> 0 then
          x^.sigh := x^.sigh or 1;
        x^.sigl := x^.sigl shl 1;
        Dec(x^.exp);
        end;
      end;

    x^.exp := x^.exp + bias;

    Result := TAG_Valid;
  end;

function FPU_u_sub (arg1: PFPU_REG; arg2: PFPU_REG; dest: PFPU_REG;
  control_w: u16; sign: u_char; expa: s32; expb: s32): intg;
  var
    shifted, answ: FPU_REG;
    extent: u32;
    ediff, ed2, borrow: intg;
  begin
    ediff := expa - expb;
{$ifdef PARANOID}
  if ( ediff < 0 )
    begin
      EXCEPTION(EX_INTERNAL|$206);
      return -1;
    end;
{$endif}

    answ.exp := expa;

{$ifdef PARANOID}
  if ( !(arg1^.sigh  and $8000) or !(arg2^.sigh  and $8000) )
    begin
      EXCEPTION(EX_INTERNAL|$209);
      return -1;
    end;
{$endif}

    if (ediff = 0) then
      begin
      shifted.sigl := arg2^.sigl;
      shifted.sigh := arg2^.sigh;
      extent := 0;
      end
    else
      if (ediff < 32) then
        begin
        ed2 := 32 - ediff;
        extent := arg2^.sigl shl ed2;
        shifted.sigl := arg2^.sigl shr ediff;
        shifted.sigl := shifted.sigl or (arg2^.sigh shl ed2);
        shifted.sigh := arg2^.sigh shr ediff;
        end
      else
        if (ediff < 64) then
          begin
          Dec(ediff, 32);
          if (ediff) = 0 then
            begin
            extent := arg2^.sigl;
            shifted.sigl := arg2^.sigh;
            shifted.sigh := 0;
            end
          else
            begin
            ed2 := 32 - ediff;
            extent := arg2^.sigl shr ediff;
            extent := extent or (arg2^.sigh shl ed2);
            if (arg2^.sigl shl ed2) <> 0 then
              extent := extent or 1;
            shifted.sigl := arg2^.sigh shr ediff;
            shifted.sigh := 0;
            end;
          end
        else
          begin
          Dec(ediff, 64);
          if (ediff) = 0 then
            begin
            extent := arg2^.sigh;
            if (arg2^.sigl) <> 0 then
              extent := extent or 1;
            shifted.sigl := 0;
            shifted.sigh := 0;
            end
          else
            begin
            if (ediff < 32) then
              begin
              extent := arg2^.sigh shr ediff;
              if ((arg2^.sigl <> 0) or ((arg2^.sigh shl (32 - ediff)) <> 0)) then
                extent := extent or 1;
              end
            else
              extent := 1;
            shifted.sigl := 0;
            shifted.sigh := 0;
            end;
          end;

    extent := -extent;
    borrow := extent;
    answ.sigl := arg1^.sigl - shifted.sigl;
    if (answ.sigl > arg1^.sigl) then
      begin
      if (borrow) <> 0 then
        Dec(answ.sigl);
      borrow := 1;
      end
    else
      if (borrow) <> 0 then
        begin
        Dec(answ.sigl);
        if (answ.sigl <> $ffffffff) then
          borrow := 0;
        end;
    answ.sigh := arg1^.sigh - shifted.sigh;
    if (answ.sigh > arg1^.sigh) then
      begin
      if (borrow) <> 0 then
        Dec(answ.sigh);
      borrow := 1;
      end
    else
      if (borrow) <> 0 then
        begin
        Dec(answ.sigh);
        if (answ.sigh <> $ffffffff) then
          borrow := 0;
        end;

{$ifdef PARANOID}
  if ( borrow )
    begin
      (* This can only occur if the code is bugged *)
      EXCEPTION(EX_INTERNAL|$212);
      return -1;
    end;
{$endif}

    if (answ.sigh and $80000000) <> 0 then
      begin
      (*
  The simpler '*dest := answ' is broken in gcc
      *)
      dest^.exp := answ.exp;
      dest^.sigh := answ.sigh;
      dest^.sigl := answ.sigl;
      Result := FPU_round(dest, extent, 0, control_w, sign);
      exit;
      end;

    if (answ.sigh = 0) then
      begin
      if (answ.sigl) <> 0 then
        begin
        answ.sigh := answ.sigl;
        answ.sigl := extent;
        extent := 0;
        Dec(answ.exp, 32);
        end
      else
        if (extent) <> 0 then
          begin
(*
 *   A rare case, the only one which is non-zero if we got here
 *         is:           1000000 .... 0000
 *                      -0111111 .... 1111 1
 *                       -------------------- 
 *                       0000000 .... 0000 1 
 *)
          if (extent <> $8000) then
            begin
            (* This can only occur if the code is bugged *)
            //EXCEPTION(EX_INTERNAL|$210); manca
            Result := -1;
            exit;
            end;
          dest^.sigh := extent;
          extent := 0;
          dest^.sigl := extent;
          Dec(dest^.exp, 64);
          Result := FPU_round(dest, extent, 0, control_w, sign);
          end
        else
          begin
          dest^.exp := 0;
          dest^.sigl := 0;
          dest^.sigh := dest^.sigl;
          Result := TAG_Zero;
          Exit;
          end;
      end;

    while ((answ.sigh and $80000000) = 0) do
      begin
      answ.sigh := answ.sigh shl 1;
      if (answ.sigl and $80000000) <> 0 then
        answ.sigh := answ.sigh or 1;
      answ.sigl := answ.sigl shl 1;
      if (extent and $80000000) <> 0 then
        //!!! dacci una occhiata al codice con l'originale
        answ.sigl := answ.sigl or 1;
      extent := extent shl 1;
      Dec(answ.exp);
      end;

    dest^.exp  := answ.exp;
    dest^.sigh := answ.sigh;
    dest^.sigl := answ.sigl;

    Result := FPU_round(dest, extent, 0, control_w, sign);

  end;

function FPU_u_mul (a: PFPU_REG; b: PFPU_REG; c: PFPU_REG; cw: u16;
  sign: u_char; expon: s32): intg;
  var
    mu, ml, mi: u64;
    lh, ll, th, tl: u32;
  begin

{$ifdef PARANOID}
  if ( ! (a^.sigh  and $8000) or ! (b^.sigh  and $8000) )
    begin
      EXCEPTION(EX_INTERNAL|$205);
    end;
{$endif}

    ml := a^.sigl;
    ml := ml * b^.sigl;
    ll := ml;
    lh := ml shr 32;

    mu := a^.sigh;
    mu := mu * b^.sigh;

    mi := a^.sigh;
    mi := mi * b^.sigl;
    tl := mi;
    th := mi shr 32;
    lh := lh + tl;
    if (tl > lh) then
      Inc(mu);
    Inc(mu, th);

    mi := a^.sigl;
    mi := mi * b^.sigh;
    tl := mi;
    th := mi shr 32;
    lh := lh + tl;
    if (tl > lh) then
      Inc(mu);
    mu := mu + th;

    ml := lh;
    ml := ml shl 32;
    ml := ml + ll;

    Dec(expon, EXP_BIAS - 1);
    if (expon <= EXP_WAY_UNDER) then
      expon := EXP_WAY_UNDER;

    c^.exp := expon;

    if ((mu and $800000000000) = 0) then
      begin
      mu := mu shl 1;
      if (ml and $800000000000) <> 0 then
        mu := mu or 1;
      ml := ml shl 1;
      Dec(c^.exp);
      end;

    ll := ml;
    lh := ml shr 32;

    if (ll) <> 0 then
      lh := lh or 1;

    c^.sigl := mu;
    c^.sigh := mu shr 32;

    Result := FPU_round(c, lh, 0, cw, sign);

  end;

function FPU_u_div (a: PFPU_REG; b: PFPU_REG; dest: PFPU_REG;
  control_w: u16; sign: u_char): intg;
  var
    exp:  s32;
    divr32, rem, rat1, rat2, work32, accum3, prodh: u32;
    work64, divr64, prod64, accum64: u64;
    ovfl: u8;
  begin

    exp := s32(a^.exp) - s32(b^.exp);

    if (exp < EXP_WAY_UNDER) then
      exp := EXP_WAY_UNDER;

    dest^.exp := exp;
{$ifdef PARANOID}
  if ( !(b^.sigh  and $8000) )
    begin
      EXCEPTION(EX_INTERNAL|$202);
    end;
{$endif}

    work64 := significand(a)^;

  (* We can save a lot of time if the divisor has all its lowest
     32 bits equal to zero. *)
    if (b^.sigl = 0) then
      begin
      divr32 := b^.sigh;
      ovfl := u8(a^.sigh >= divr32);
      rat1 := work64 div divr32;
      rem  := work64 mod divr32;
      work64 := rem;
      work64 := work64 shl 32;
      rat2 := work64 div divr32;
      rem  := work64 mod divr32;

      work64 := rem;
      work64 := work64 shl 32;
      rem := work64 div divr32;

      if (ovfl) <> 0 then
        begin
        rem := rem shr 1;
        if (rat2 and 1) <> 0 then
          rem := rem or $8000;
        rat2 := rat2 shr 1;
        if (rat1 and 1) <> 0 then
          rat2 := rat2 or $8000;
        rat1 := rat1 shr 1;
        rat1 := rat1 or $8000;
        Inc(dest^.exp);
        end;
      dest^.sigh := rat1;
      dest^.sigl := rat2;

      Dec(dest^.exp);
      Result := FPU_round(dest, rem, 0, control_w, sign);
      Exit;
      end;

    (* This may take a little time... *)

    accum64 := work64;
    divr64 := significand(Pfpu_reg(b))^;
    ovfl := u8(accum64 >= divr64);
    if (ovfl) <> 0 then
      Dec(accum64, divr64);
    divr32 := b^.sigh + 1;

    if (divr32 <> 0) then
      begin
      rat1 := accum64 div divr32;
      end
    else
      rat1 := accum64 shr 32;
    prod64 := rat1 * u64(b^.sigh);

    Dec(accum64, prod64);
    prod64 := rat1 * u64(b^.sigl);
    accum3 := prod64;
    if (accum3) <> 0 then
      begin
      accum3 := -accum3;
      Dec(accum64);
      end;
    prodh := prod64 shr 32;
    Dec(accum64, prodh);

    work32 := accum64 shr 32;
    if (work32) <> 0 then
      begin
{$ifdef PARANOID}
      if ( work32 !:= 1 )
	begin
	  EXCEPTION(EX_INTERNAL|$203);
	end;
{$endif}

      (* Need to subtract the divisor once more. *)
      work32 := accum3;
      accum3 := work32 - b^.sigl;
      if (accum3 > work32) then
        Dec(accum64);
      Inc(rat1);
      Dec(accum64, b^.sigh);

{$ifdef PARANOID}
      if ( (accum64 shr 32) )
	begin
	  EXCEPTION(EX_INTERNAL|$203);
	end;
{$endif}
      end;

  (* Now we essentially repeat what we have just done, but shifted
     32 bits. *)

    accum64 := accum64 shl 32;
    accum64 := accum64 or accum3;
    if (accum64 >= divr64) then
      begin
      Dec(accum64, divr64);
      Inc(rat1);
      end;
    if (divr32 <> 0) then
      begin
      rat2 := accum64 div divr32;
      end
    else
      rat2 := accum64 shr 32;
    prod64 := rat2 * u64(b^.sigh);

    Dec(accum64, prod64);
    prod64 := rat2 * u64(b^.sigl);
    accum3 := prod64;
    if (accum3) <> 0 then
      begin
      accum3 := -accum3;
      Dec(accum64);
      end;
    prodh := prod64 shr 32;
    Dec(accum64, prodh);

    work32 := accum64 shr 32;
    if (work32) <> 0 then
      begin
{$ifdef PARANOID}
      if ( work32 !:= 1 )
	begin
	  EXCEPTION(EX_INTERNAL|$203);
	end;
{$endif}

      (* Need to subtract the divisor once more. *)
      work32 := accum3;
      accum3 := work32 - b^.sigl;
      if (accum3 > work32) then
        Dec(accum64);
      Inc(rat2);
      if (rat2 = 0) then
        Inc(rat1);
      Dec(accum64, b^.sigh);

{$ifdef PARANOID}
      if ( (accum64 shr 32) )
	begin
	  EXCEPTION(EX_INTERNAL|$203);
	end;
{$endif}
      end;

    (* Tidy up the remainder *)

    accum64 := accum64 shl 32;
    accum64 := accum64 or accum3;
    if (accum64 >= divr64) then
      begin
      Dec(accum64, divr64);
      Inc(rat2);
      if (rat2 = 0) then
        begin
        Inc(rat1);
{$ifdef PARANOID}
	  (* No overflow should be possible here *)
	  if ( rat1 = 0 )
	    begin
	      EXCEPTION(EX_INTERNAL|$203);
	    end;
{$endif}
        end;

      end;

  (* The basic division is done, now we must be careful with the
     remainder. *)

    if (ovfl) <> 0 then
      begin
      if (rat2 and 1) <> 0 then
        rem := $8000
      else
        rem := 0;
      rat2 := rat2 shr 1;
      if (rat1 and 1) <> 0 then
        rat2 := rat2 or $8000;
      rat1 := rat1 shr 1;
      rat1 := rat1 or $8000;

      if (accum64) <> 0 then
        rem := rem or $ff0000;

      Inc(dest^.exp);
      end
    else
      begin
      (* Now we just need to know how large the remainder is
   relative to half the divisor. *)
      if (accum64 = 0) then
        rem := 0
      else
        begin
        accum3 := accum64 shr 32;
        if (accum3 and $8000) <> 0 then
          begin
          (* The remainder is definitely larger than 1/2 divisor. *)
          rem := $ff000000;
          end
        else
          begin
          accum64 := accum64 shl 1;
          if (accum64 >= divr64) then
            begin
            Dec(accum64, divr64);
            if (accum64 = 0) then
              rem := $8000
            else
              rem := $ff000000;
            end
          else
            rem := $7f000000;
          end;
        end;
      end;

    dest^.sigh := rat1;
    dest^.sigl := rat2;

    Dec(dest^.exp);
    Result := FPU_round(dest, rem, 0, control_w, sign);

  end;

function FPU_u_add (arg1: PFPU_REG; arg2: PFPU_REG; answ: PFPU_REG;
  control_w: u16; sign: u_char; expa: s32; expb: s32): intg;
  var
    rtmp: PFPU_REG;
    shifted: FPU_REG;
    extent: u32;
    ediff, ed2, eflag, ovfl, carry: intg;
  begin
    ediff  := expa - expb;
    extent := 0;
    if (ediff < 0) then
      begin
      ediff := -ediff;
      rtmp  := arg1;
      arg1  := arg2;
      arg2  := rtmp;
      expa  := expb;
      end;

    (* Now we have exponent of arg1 >= exponent of arg2 *)

    answ^.exp := expa;

{$ifdef PARANOID}
  if ( !(arg1^.sigh  and $8000) or !(arg2^.sigh  and $8000) )
    begin
      EXCEPTION(EX_INTERNAL|$201);
      return -1;
    end;
{$endif}

    if (ediff = 0) then
      begin
      extent := 0;
      shifted.sigl := arg2^.sigl;
      shifted.sigh := arg2^.sigh;
      end
    else
      if (ediff < 32) then
        begin
        ed2 := 32 - ediff;
        extent := arg2^.sigl shl ed2;
        shifted.sigl := arg2^.sigl shr ediff;
        shifted.sigl := shifted.sigl or (arg2^.sigh shl ed2);
        shifted.sigh := arg2^.sigh shr ediff;
        end
      else
        if (ediff < 64) then
          begin
          Dec(ediff, 32);
          if (ediff) = 0 then
            begin
            eflag  := 0;
            extent := arg2^.sigl;
            shifted.sigl := arg2^.sigh;
            end
          else
            begin
            ed2 := 32 - ediff;
            eflag := arg2^.sigl;
            if (eflag) <> 0 then
              extent := extent or 1;
            extent := arg2^.sigl shr ediff;
            extent := extent or (arg2^.sigh shl ed2);
            shifted.sigl := arg2^.sigh shr ediff;
            end;
          shifted.sigh := 0;
          end
        else
          begin
          Dec(ediff, 64);
          if (ediff) = 0 then
            begin
            eflag  := arg2^.sigl;
            extent := arg2^.sigh;
            end
          else
            begin
            ed2 := 64 - ediff;
            eflag := arg2^.sigl or arg2^.sigh;
            extent := arg2^.sigh shr ediff;
            end;
          shifted.sigl := 0;
          shifted.sigh := 0;
          if (eflag) <> 0 then
            extent := extent or 1;
          end;

    answ^.sigh := arg1^.sigh + shifted.sigh;
    ovfl := u8(shifted.sigh > answ^.sigh);
    answ^.sigl := arg1^.sigl + shifted.sigl;
    if (shifted.sigl > answ^.sigl) then
      begin
      Inc(answ^.sigh);
      if (answ^.sigh = 0) then
        ovfl := 1;
      end;
    if (ovfl) <> 0 then
      begin
      carry  := extent and 1;
      extent := extent shr 1;
      extent := extent or carry;
      if (answ^.sigl and 1) <> 0 then
        extent := extent or $8000;
      answ^.sigl := answ^.sigl shr 1;
      if (answ^.sigh and 1) <> 0 then
        answ^.sigl := answ^.sigl or $8000;
      answ^.sigh := answ^.sigh shr 1;
      answ^.sigh := answ^.sigh or $8000;
      Inc(answ^.exp);
      end;

    Result := FPU_round(answ, extent, 0, control_w, sign);
  end;

function wm_sqrt (n: PFPU_REG; dummy1: intg; dummy2: intg; control_w: u16;
  sign: u_char): intg;
  var
    nn, guess, halfn, lowr, mid, upr, diff, uwork: u64;
    work: s64;
    ne, guess32, work32, diff32, mid32: u32;
    shifted: intg;
  begin

    nn := significand(Pfpu_reg(n))^;
    ne := 0;
    if (exponent16(Pfpu_reg(n)) = EXP_BIAS) then
      begin
      (* Shift the argument right one position. *)
      if (nn and 1) <> 0 then
        ne := $8000;
      nn := nn shr 1;
      guess := n^.sigh shr 2;
      shifted := 1;
      end
    else
      begin
      guess := n^.sigh shr 1;
      shifted := 0;
      end;

    Inc(guess, $40000000);
    guess := guess * $aaaaaaaa;
    guess := guess shl 1;
    guess32 := guess shr 32;
    if ((guess32 and $8000)) = 0 then
      guess32 := $8000;
    halfn := nn shr 1;

    guess32 := halfn div guess32 + (guess32 shr 1);
    guess32 := halfn div guess32 + (guess32 shr 1);
    guess32 := halfn div guess32 + (guess32 shr 1);


(*
 * Now that an estimate accurate to about 30 bits has been obtained,
 * we improve it to 60 bits or so.
 *
 * The strategy from now on is to compute new estimates from
 *      guess ::= guess + (n - guess^2) / (2 * guess)
 *)

    work := guess32;
    work := nn - work * guess32;
    work := work shl 28;       (* 29 - 1 *)
    work := work div guess32;
    work := work shl 3;        (* 29 + 3 := 32 *)
    Inc(work, u64(guess32) shl 32);

    if (work = 0) then  (* This happens in one or two special cases *)
      work := ($ffffffffffffffff);

    guess := work;

    (* guess is now accurate to about 60 bits *)


    if (work > 0) then
      begin
{$ifdef PARANOID}
      if ( (n^.sigh !:= $ffffffff) @ and (n^.sigl !:= $ffffffff) )
	begin
	  EXCEPTION(EX_INTERNAL|$213);
	end;
{$endif}
      (* We know the answer here. *)
      Result := FPU_round(n, $7fffffff, 0, control_w, sign);
      Exit;
      end;

    (* Refine the guess to significantly more than 64 bits. *)

    (* First, square the current guess. *)

    guess32 := guess shr 32;
    work32  := guess;

    (* lower 32 times lower 32 *)
    lowr := work32;
    lowr := lowr * work32;

    (* lower 32 times upper 32 *)
    mid := guess32;
    mid := mid * work32;

    (* upper 32 times upper 32 *)
    upr := guess32;
    upr := upr * guess32;

    (* upper 32 bits of the middle product times 2 *)
    Inc(Upr, mid shr (32 - 1));

    (* lower 32 bits of the middle product times 2 *)
    work32 := mid shl 1;

    (* upper 32 bits of the lower product *)
    mid32 := lowr shr 32;
    mid32 := mid32 + work32;
    if (mid32 < work32) then
      Inc(upr);

    (* We now have the first 96 bits (truncated) of the square of the guess *)

    diff := upr - nn;
    diff32 := mid32 - ne;
    if (diff32 > mid32) then
      Dec(diff);

    if (s64(diff) < 0) then
      begin
      (* The difference is negative, negate it. *)
      diff32 := -s32(diff32);
      diff := not diff;
      if (diff32 = 0) then
        Inc(diff);
{$ifdef PARANOID}
      if ( (diff shr 32) !:= 0 )
	begin
	  EXCEPTION(EX_INTERNAL|$207);
	end;
{$endif}

      diff := diff shl 32;
      diff := diff or diff32;
      work32 := diff div guess32;
      work := work32;
      work := work shl 32;

      diff := diff mod guess32;
      diff := diff shl 32;
      work32 := diff div guess32;

      work := work or work32;

      work := work shr 1;
      work32 := work shr 32;

      Inc(guess, work32);
      guess32 := work;        (* The next 32 bits *)
      (* The guess should now be good to about 90 bits *)
      end
    else
      begin
      (* The difference is positive. *)
      diff := diff shl 32;
      diff := diff or diff32;

      work32 := diff div guess32;
      work := work32;
      work := work shl 32;

      diff := diff mod guess32;
      diff := diff shl 32;
      work32 := diff div guess32;

      work := work or work32;

      work := work shr 1;
      work32 := work shr 32;

      guess32 := work;        (* The last 32 bits (of 96) *)
      guess32 := -guess32;
      if (guess32) <> 0 then
        Dec(guess);
      Dec(guess, work32);       (* The first 64 bits *)
      (* The guess should now be good to about 90 bits *)
      end;


    n^.exp := 0;

    if (guess32 >= u32(-ERR_MARGIN)) then
      begin
      (* Nearly exact, we round the 64 bit result upward. *)
      Inc(guess);
      end
    else
      if ((guess32 > ERR_MARGIN) and ((guess32 < $8000 - ERR_MARGIN) or
        (guess32 > $8000 + ERR_MARGIN))) then
        begin
        (* We have enough accuracy to decide rounding *)
        n^.sigh := guess; // significand(Pfpu_reg(n)) := guess;
        Result  := FPU_round(n, guess32, 0, control_w, sign);
        Exit;
        end;

    if ((guess32 <= ERR_MARGIN) or (guess32 >= u32(-ERR_MARGIN))) then
      begin
      (*
       * This is an easy case because x^1/2 is monotonic.
       * We need just find the square of our estimate, compare it
       * with the argument, and deduce whether our estimate is
       * above, below, or exact. We use the fact that the estimate
       * is known to be accurate to about 90 bits.
       *)


      (* We compute the lower 64 bits of the 128 bit product *)
      work32 := guess;
      lowr := work32;
      lowr := lowr * work32;

      uwork  := guess shr 32;
      work32 := guess;
      uwork  := uwork * work32;
      uwork  := uwork shl 33;   (* 33 := 32+1 (for two times the product) *)

      Inc(lowr, uwork);

      (* We need only look at bits 65..96 of the square of guess. *)
      if (shifted) <> 0 then
        work32 := lowr shr 31
      else
        work32 := lowr shr 32;

{$ifdef PARANOID}
      if ( ((s32)work32 > 3*ERR_MARGIN) or ((s32)work32 < -3*ERR_MARGIN) )
	begin
	  EXCEPTION(EX_INTERNAL|$214);
	end;
{$endif}

      n^.sigh := guess;
      if (s32(work32) > 0) then
        begin
        (* guess is too large *)
        Dec(n^.sigh);
        Result := FPU_round(n, $ffffff00, 0, control_w, sign);
        Exit;
        end
      else
        if (s32(work32) < 0) then
          begin
          (* guess is a little too small *)
          Result := FPU_round(n, $000000ff, 0, control_w, sign);
          Exit;
          end

        else
          if (u32(lowr) <> 0) then
            begin

            (* guess is too large *)
            Dec(n^.sigh);
            Result := FPU_round(n, $ffffff00, 0, control_w, sign);
            Exit;
            end;

      (* Our guess is precise. *)
      Result := FPU_round(n, 0, 0, control_w, sign);
      Exit;
      end;

  (* Very similar to the case above, but the last bit is near 0.5.
     We handle this just like the case above but we shift everything
     by one bit. *)


    uwork := guess;
    uwork := uwork shl 1;
    uwork := uwork or 1;      (* add the half bit *)

    (* We compute the lower 64 bits of the 128 bit product *)
    work32 := uwork;
    lowr := work32;
    lowr := lowr * work32;

    work32 := uwork shr 32;
    uwork  := uwork and $ffffffff;
    uwork  := uwork * work32;
    uwork  := uwork shl 33;   (* 33 := 32+1 (for two times the product) *)

    Inc(lowr, uwork);

  (* We now have the 64 bits. The lowest 32 bits of lowr
                    are not all zero (the lsb is 1). *)

    (* We need only look at bits 65..96 of the square of guess. *)
    if (shifted) <> 0 then
      work32 := lowr shr 31
    else
      work32 := lowr shr 32;

{$ifdef PARANOID}
  if ( ((s32)work32 > 4*3*ERR_MARGIN) or ((s32)work32 < -4*3*ERR_MARGIN) )
    begin
      EXCEPTION(EX_INTERNAL|$215);
    end;
{$endif}

    n^.sigh := guess;
    if (s32(work32) < 0) then
      begin
      (* guess plus half bit is a little too small *)
      Result := FPU_round(n, $800000ff, 0, control_w, sign);
      Exit;
      end
    else (* Note that the lower 64 bits of the product are not all zero *)
      begin
      (* guess plus half bit is too large *)
      Result := FPU_round(n, $7fffff00, 0, control_w, sign);
      Exit;
      end;

  (*
    Note that the result of a square root cannot have precisely a half bit
    of a least significant place (it is left as an exercise for the reader
    to prove this! (hint: 65 bit*65 bit :=> n bits)).
  *)

  end;

function FPU_shrx (arg1: pointer; arg2: u32): intg;
  var
    x: u32;
  begin
    if (arg2 >= 64) then
      begin
      if (arg2 >= 96) then
        begin
        pu64(arg1)^ := 0;
        Result := 0;
        Exit;
        end;
      Dec(arg2, 64);
      x := pu64(arg1)^ shr 32;
      pu64(arg1)^ := 0;

      if (arg2) <> 0 then
        begin
        Result := x shr arg2;
        exit;
        end
      else
        begin
        Result := x;
        Exit;
        end;
      end;

    if (arg2 < 32) then
      begin
      if (arg2 = 0) then
        begin
        Result := 0;
        Exit;
        end;

      x := pu64(arg1)^ shl (32 - arg2);
      end
    else
      if (arg2 > 32) then
        begin
        x := pu64(arg1)^ shr (arg2 - 32);
        end
      else
        begin
        (* arg2 = 32 *)
        x := pu64(arg1)^;
        end;

    pu64(arg1)^ := pu64(arg1)^ shr arg2;

    Result := x;
  end;

function FPU_shrxs (arg1: pointer; arg2: u32): intg;
  var
    x, bits: u32;
    lost: u64;
    v: int64;
  begin

    if (arg2 >= 64) then
      begin
      if (arg2 >= 96) then
        begin
        bits := u32(pu64(arg1)^ <> 0);
        pu64(arg1)^ := 0;
        Result := ifthen(bits <> 0, 1, 0);
        Exit;
        end;
      Dec(arg2, 64);
      lost := pu64(arg1)^ shl (32 - arg2);
      x := pu64(arg1)^ shr 32;
      pu64(arg1)^ := 0;

      if (arg2) <> 0 then
        x := x shr arg2;

      if (lost) <> 0 then
        x := x or 1;

      Result := x;
      end;

    if (arg2 < 32) then
      begin
      if (arg2 = 0) then
        (* No bits are lost *)
        begin
        Result := 0;
        Exit;
        end;

      (* No bits are lost *)
      x := pu64(arg1)^ shl (32 - arg2);
      end
    else
      if (arg2 > 32) then
        begin
        bits := pu64(arg1)^;
        bits := bits shl (64 - arg2);
        x := pu64(arg1)^ shr (arg2 - 32);
        if (bits) <> 0 then
          x := x or 1;
        end
      else
        begin
        (* arg2 = 32 *)
        (* No bits are lost *)
        x := pu64(arg1)^;
        end;

    v := pu64(arg1)^;
    v := v shr arg2;
    pu64(arg1)^ := (v * -1) + 2;

    if (x and $7fffffff) <> 0 then
      x := x or 1;

    Result := x;
  end;

function FPU_div_small (x: pu64; y: u32): intg;
  var
    retval: u32;
  begin

    retval := x^ mod y;

    x^ := x^ div y;

    Result := retval;
  end;

function exponent (x: Pfpu_reg): longword;
  begin
    Result := (x^.exp and $7fff) - EXTENDED_Ebias;
  end;

function top: u32;
  begin
    Result := I387.soft.ftop;
  end;

function no_ip_update: u_char;
  begin
    Result := I387.soft.no_update;
  end;

function register_base: pu_char;
  begin
    Result := pu_char(parray64(@I387.soft.st_space));
  end;

function fpu_register (index: integer): PFPU_REG;
  begin
    Result := PFPU_REG(longint(register_base) + (sizeof(FPU_REG) * index));
  end;

procedure setcc (cc: u32);
  begin
    I387.soft.swd := I387.soft.swd and not (SW_C0 or SW_C1 or SW_C2 or SW_C3);
    I387.soft.swd := I387.soft.swd or cc and (SW_C0 or SW_C1 or SW_C2 or SW_C3);
  end;

procedure poppop;
  begin
    FPU_pop;
    FPU_pop;
  end;

procedure FPU_illegal;
  begin
    math_abort(nil, 0);
  end;

procedure pop_0 ();
  begin
    FPU_settag0(TAG_Empty);
    Inc(I387.soft.ftop);
  end;

function status_word: intg;
  begin
    Result := ((I387.soft.swd and not SW_Top and $ffff) or
      ((top shl SW_Top_Shift) and SW_Top));
  end;

procedure push;
  begin
    Inc(I387.soft.ftop);
  end;

procedure add_two_Xsig (dest: PXsig; const x2: PXsig; exp: ps32);
  var
    ovfl: intg;
  begin
    ovfl := 0;

    dest^.lsw := dest^.lsw + x2^.lsw;
    if (dest^.lsw < x2^.lsw) then
      begin
      Inc(dest^.midw);
      if (dest^.midw = 0) then
        begin
        Inc(dest^.msw);
        if (dest^.msw = 0) then
          ovfl := 1;
        end;
      end;
    dest^.midw := dest^.midw + x2^.midw;
    if (dest^.midw < x2^.midw) then
      begin
      Inc(dest^.msw);
      if (dest^.msw = 0) then
        ovfl := 1;
      end;
    dest^.msw := dest^.msw + x2^.msw;
    if (dest^.msw < x2^.msw) then
      ovfl := 1;
    if (ovfl) <> 0 then
      begin
      Inc(exp^);
      dest^.lsw := dest^.lsw shr 1;
      if (dest^.midw and 1) <> 0 then
        dest^.lsw := dest^.lsw or $8000;
      dest^.midw := dest^.midw shr 1;
      if (dest^.msw and 1) <> 0 then
        dest^.midw := dest^.midw or $8000;
      dest^.msw := dest^.msw shr 1;
      dest^.msw := dest^.msw or $8000;
      end;
  end;

procedure negate_Xsig (x: PXsig);
  begin
    x^.lsw  := not x^.lsw;
    x^.midw := not x^.midw;
    x^.msw  := not x^.msw;
    Inc(x^.lsw);
    if (x^.lsw = 0) then
      begin
      Inc(x^.midw);
      if (x^.midw = 0) then
        Inc(x^.msw);
      end;
  end;

function XSIG_LL (x: Xsig): pu64;
  begin
    Result := pu64(x.msw);
  end;

function LL_MSW (x: Pointer): pu64;
  begin
    Result := pu64(x);
  end;

procedure add_Xsig_Xsig (dest: PXsig; const x2: PXsig);
  begin
    dest^.lsw := dest^.lsw + x2^.lsw;
    if (dest^.lsw < x2^.lsw) then
      begin
      Inc(dest^.midw);
      if (dest^.midw = 0) then
        Inc(dest^.msw);
      end;
    dest^.midw := dest^.midw + x2^.midw;
    if (dest^.midw < x2^.midw) then
      begin
      Inc(dest^.msw);
      end;
    dest^.msw := dest^.msw + x2^.msw;
  end;

function mul_32_32 (const arg1: u32; const arg2: u32): u32;
  begin
    Result := (u64(arg1) * arg2) shr 32;
  end;

function access_limit: pu_char;
  begin
    Result := pu_char(@I387.soft.alimit);
  end;

end.
