unit div_xsig_unit;

interface

uses fpu_types;

{$include fpu_const.pas}

procedure div_Xsig (const aa: PXsig; const b: PXsig; dest: PXsig);

implementation

uses mul_Xsig;

procedure div_Xsig (const aa: PXsig; const b: PXsig; dest: PXsig);
  var
    a, xpr, res: Xsig;
    prodh, prodl, den, wd: u32;
    num, prod: u64;
  begin
    a := aa^;

{$ifdef PARANOID}
  if ( (b^.msw  and $8000) = 0 )
    begin
      EXCEPTION(EX_INTERNAL|$240);
      exit;
    end;
{$endif}

    (* Shift a right *)
    a.lsw := a.lsw shr 1;
    if (a.midw and 1) <> 0 then
      a.lsw := a.lsw or $8000;
    a.midw := a.midw shr 1;
    if (a.msw and 1) <> 0 then
      a.midw := a.midw or $8000;
    a.msw := a.msw shr 1;

    num := a.msw;
    num := num shl 32;
    num := num or a.midw;

    den := b^.msw + 1;
    if (den) <> 0 then
      begin
      res.msw := num div den;
      end
    else
      res.msw := a.msw;

    xpr := b^;
    mul32_Xsig(@xpr, res.msw);
    a.msw := a.msw - xpr.msw;
    wd := a.midw;
    a.midw := a.midw - xpr.midw;
    if (a.midw > wd) then
      Dec(a.msw);
    wd := a.lsw;
    a.lsw := a.lsw - xpr.lsw;
    if (a.lsw > wd) then
      begin
      Dec(a.midw);
      if (a.midw = $ffffffff) then
        Dec(a.msw);
      end;

{$ifdef PARANOID}
      if ( a.msw > 1 )
	begin
	  EXCEPTION(EX_INTERNAL|$241);
	end;
{$endif}

    while ((a.msw <> 0) or ((a.midw > b^.msw))) do
      begin
      wd := a.midw;
      a.midw := a.midw - b^.msw;
      if (a.midw > wd) then
        Dec(a.msw);
      wd := a.lsw;
      a.lsw := a.lsw - b^.midw;
      if (a.lsw > wd) then
        begin
        Dec(a.midw);
        if (a.midw = $ffffffff) then
          Dec(a.msw);
        end;
      Inc(res.msw);
      end;

    (* Whew! result.msw is now done. *)

    num := a.midw;
    num := num shl 32;
    num := num or a.lsw;

    if (den) <> 0 then
      begin
      res.midw := num div den;
      end
    else
      res.midw := a.midw;

    prod := res.midw;
    prod := prod * b^.msw;
    a.midw := a.midw - prod shr 32;
    prodl := prod;
    wd := a.lsw;
    a.lsw := a.lsw - prodl;
    if (a.lsw > wd) then
      Dec(a.midw);

    prod := res.midw;
    prod := prod * b^.midw;
    prodh := prod shr 32;
    wd := a.lsw;
    a.lsw := a.lsw - prodh;
    if (a.lsw > wd) then
      Dec(a.midw);

{$ifdef PARANOID}
      if ( a.midw > 1 ) then
	begin
          fpu_EXCEPTION(EX_INTERNAL or $242);
	end;
{$endif}

    while ((a.midw <> 0) or (a.lsw > b^.msw)) do
      begin
      wd := a.lsw;
      a.lsw := a.lsw - b^.msw;
      if (a.lsw > wd) then
        Dec(a.midw);
      Inc(res.midw);
      end;


    (* Now result.msw is done, the lsw is next... *)

    num := a.lsw;
    num := num shl 32;

    if (den) <> 0 then
      begin
      res.lsw := num div den;
      end
    else
      res.lsw := a.lsw;

    prod  := res.lsw;
    prod  := prod * b^.msw;
    a.lsw := a.lsw - prod shr 32;

{$ifdef PARANOID}
  if ( a.lsw > 2 ) then
    begin
      EXCEPTION(EX_INTERNAL|$243);
    end;
{$endif}

    res.lsw := res.lsw - a.lsw;

    (* Hey! we're done. *)

    dest^ := res;

  end;

end.
