unit reg_u_div;

interface

function FPU_u_div(a:PFPU_REG; b:PFPU_REG; dest:PFPU_REG;control_w:u16; sign:u8):intg;

implementation

function FPU_u_div(a:PFPU_REG; b:PFPU_REG; dest:PFPU_REG;control_w:u16; sign:u8):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 $80000000)=0 ) then
    begin
      fpu_EXCEPTION(EX_INTERNAL or $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 := Word(a^.sigh >= divr32);
      rat1 := Trunc(work64 / 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 $80000000;
        	  rat2 :=	  rat2 shr 1;
        	  if ( rat1  and 1 ) <>0 then
        	    rat2 :=	    rat2 or $80000000;
        	  rat1 :=	  rat1 shr 1;
        	  rat1 :=	  rat1 or $80000000;
        	  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(b);

  ovfl:=Word(accum64 >= divr64);
  if ( (ovfl)<>0 ) then
    accum64 :=  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;

  accum64 :=  accum64 - prod64;
  prod64 := rat1 * (u64)b^.sigl;
  accum3 := prod64;
  if ( accum3 )<>0 then
    begin
      accum3 := -accum3;
      dec(accum64);
    end;
  prodh := prod64 shr 32;
  accum64 :=  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);
      accum64 :=      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
      accum64 :=      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;

  accum64 :=  accum64 - prod64;
  prod64 := rat2 * (u64)b^.sigl;
  accum3 := prod64;
  if ( accum3 )
    begin
      accum3 := -accum3;
      dec(accum64);
    end;
  prodh := prod64 shr 32;
  accum64 :=  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);
      accum64 :=      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
      accum64 :=      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;
	end;
{$endif}
    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 := $80000000;
      else
      	rem := 0;
      rat2 :=      rat2 shr 1;
      if ( rat1  and 1 )<>0 then
      	rat2 :=	rat2 or $80000000;
      rat1 :=      rat1 shr 1;
      rat1 :=      rat1 or $80000000;

      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 $80000000 )<>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
		        accum64 :=		  accum64 - divr64;
      		  if ( accum64 = 0 ) then
      		    rem := $80000000
		            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;

end.
