{$I stdefine.inc}
{
    References:
      http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html
      Search for IEEE-754...
}

{AH   2-02-05 stuff to handle Comp/tQuad
     28.02.08 renamed from int64 to mathx
     17.05.09 Removed the stringsx uses clause by adding HexStr[]
}

// the const seems not to work with tquad (it's always copied to stack)
// if it should work someday the code needs to be changed

Unit
  mathx; {Also see Long32 in the VP RTL}

interface

uses math {MaxComp};

type
  CReal = Extended;
  PCReal = ^CReal;

  TQuad      = Comp;
  TQuadRec   = packed record
		Lo, Hi : longint;
	       end;
  Integer64  = TQuad;

  tCondition = (o,no,b,ae,e,ne,be,a,s,ns,p,np,l,ge,le,g);

{* storing of Extended and retrieving of tQuad rounds/truncs automatically
   fld tw ptr x -> fistp qw ptr y
 * tQuad := Extended; is much shorter than inline or function call, but
 * Extended := TQuad(Extended); is not allowed
}
function trunc64( const Value : Extended ) : tQuad;
function round64( const Value : Extended ) : tQuad;
function or64(const Op1, Op2: TQuad) : tQuad;
function and64(const Op1, Op2: TQuad) : tQuad;
function xor64(const Op1, Op2: TQuad) : tQuad;
function shl64(const Op1, Op2: TQuad) : tQuad;
function shr64(const Op1, Op2: TQuad) : tQuad;
function not64(const Op1: TQuad) : tQuad;
function neg64(const Op1: TQuad) : tQuad;
function add64(const Op1, Op2: TQuad) : tQuad;
function sub64(const Op1, Op2: TQuad) : tQuad;
function mul64(const Op1, Op2: TQuad) : tQuad;
function div32(Dividend, Divisor: LongInt): LongInt;
function div64(const Op1, Op2: TQuad) : tQuad;
function idiv32(Dividend, Divisor: LongInt): LongInt;
function idiv64(const Op1, Op2: TQuad) : tQuad;
function mod64(const Op1, Op2: TQuad) : tQuad;
function divmod32(Operand, Divisor: LongInt; var Modulus: LongInt): LongInt;

function cmp64(Op1, Op2: TQuad; Condition: tCondition) : Boolean;

function q2hex(const Value : tQuad) : shortstring;
function q2bin(const Value : tQuad) : shortstring;
function q2oct(const Value : tQuad) : shortstring;

procedure Int64Init;
procedure Int64Done;

{From math.pas, but converted to 64-Bit}
function Ceil64(X: Extended): Extended;
function Floor64(X: Extended): Extended;

implementation

{ For conversion to hex strings }
const
  HexStr : string[16] = '0123456789ABCDEF';

{ 32-Bit operands code}

{ The standard 32-bit division in VP (BP?) is an integer division, see idiv32() }
function div32(Dividend, Divisor: LongInt): LongInt; assembler; {$FRAME-} {$USES ecx,edx}
asm
    mov   eax, [Dividend]
    xor   edx, edx
    div   dword [Divisor]
end;

function idiv32(Dividend, Divisor: LongInt): LongInt; assembler; {$FRAME-} {$USES ecx,edx}
asm
    mov   eax, [Dividend]
    cdq
    idiv  dword [Divisor]
end;

{ From VP \vp\source\rtl\long32.pas:
  Unsigned long division and modulus.  Returns Modulus in the Var parameter}
function divmod32(Operand, Divisor: LongInt; var Modulus: LongInt): LongInt; assembler; {$FRAME-} {$USES ecx,edx}
asm
      mov   eax,   [Operand]
      mov   ecx,   [Divisor]
      xor   edx,   edx
      div   ecx
      mov   ecx,   [Modulus]
      mov   [ecx], edx
end;

{ 64-Bit operands code}

// from system.pas
// 80x87 Status Word
const
  mSW_IE	= $0001;	// Invalid Operation exception
  wSW_IE	= 1;
  mSW_DE	= $0002;	// Denormalized Operand exception
  wSW_DE	= 1;
  mSW_ZE	= $0004;	// Zero-Divide exception
  wSW_ZE	= 1;
  mSW_OE	= $0008;	// Overflow exception
  wSW_OE	= 1;
  mSW_UE	= $0010;	// Underflow exception
  wSW_UE	= 1;
  mSW_PE	= $0020;	// Precision exception
  wSW_PE	= 1;
  mSW_SF	= $0040;	// Stack flag (387+ only)
  wSW_SF	= 1;
  mSW_ES	= $0080;	// Error summary
  wSW_ES	= 1;
  mSW_C0	= $0100;	// Condition bit 0
  wSW_C0	= 1;
  sSW_C0	= 8;
  mSW_C1	= $0200;	// Condition bit 1
  wSW_C1	= 1;
  mSW_C2	= $0400;	// Condition bit 2
  wSW_C2	= 1;
  mSW_ST	= $3800;	// Stack top
  wSW_ST	= 3;
  sSW_ST	= 11;
  mSW_C3	= $4000;	// Condition bit 3
  wSW_C3	= 1;
  mSW_B 	= $8000;	// Busy bit
  wSW_B 	= 1;

// 80x87 Control Word

  mCW_IM	= $0001;	// Invalid Operation mask
  wCW_IM	= 1;    	// Bit = 1 if Exception is masked
  mCW_DM	= $0002;	// Denormalized Operand mask
  wCW_DM	= 1;
  mCW_ZM	= $0004;	// Zero-Divide mask
  wCW_ZM	= 1;
  mCW_OM	= $0008;	// Overflow mask
  wCW_OM	= 1;
  mCW_UM	= $0010;	// Underflow mask
  wCW_UM	= 1;
  mCW_PM	= $0020;	// Precision mask
  wCW_PM	= 1;
  mCW_PC	= $0300;	// Precision control
  wCW_PC	= 2;
  sCW_PC	= 8;
  mCW_RC	= $0C00;	// Rounding control
  wCW_RC	= 2;
  sCW_RC	= 10;
  mCW_IC	= $1000;	// Infinity control
  sCW_IC	= 12;
  wCW_IC	= 1;

  IC_Projective = 0;    	// Projective closure (387 doesn't support it)
  IC_Affine     = 1;    	// Affine mode

  RC_Nearest    = 0;    	// Rounding to nearest (the default)
  RC_Down       = 1;    	// Rounding down (towards "-" infinity)
  RC_Up 	= 2;    	// Rounding up (towards "+" infinity)
  RC_To_Zero    = 3;    	// Rounding toward zero.

  PC_Single     = 0;    	// Round to single precision
  PC_Reserved   = 1;    	// Reserved ( should not be specified)
  PC_Double     = 2;    	// Round to double precision
  PC_Extended   = 3;    	// Round to extended precision (the default)

  TAG_Valid     = 0;    	// Tag values
  TAG_Zero      = 1;
  TAG_Spec      = 2;
  TAG_Empty     = 3;

  {Default control word for 64-Bit operations}
  CWDefault: Word = (IC_Affine   shl sCW_IC) or    // $1000 Affine mode
		    (PC_Extended shl sCW_PC) or    // $0300 Round to extended
				     mCW_PM  or    // $0020 Precision Masked
				     mCW_UM  or    // $0010 Underflow Masked
				     mCW_OM  or    // $0008 Overflow Masked
				     mCW_ZM  or    // $0004 Zero-Divide Masked
				     mCW_DM  or    // $0002 Denormalized Operand Masked
				     mCW_IM ;      // $0001 Invalid Operation Masked

  {Control word for truncation}
  CWChop: Word = (IC_Affine   shl sCW_IC) or    // Affine mode
		 (RC_To_Zero  shl sCW_RC) or    // Round towards 0
		 (PC_Extended shl sCW_PC) or    // Round to extended
				  mCW_PM  or    // Masked
				  mCW_UM  or    // Masked
				  mCW_OM  or    // Masked
				  mCW_ZM  or    // Masked
				  mCW_DM  or    // Masked
				  mCW_IM ;      // Masked

function trunc64( const Value : Extended ) : tQuad; assembler; {$FRAME-}{$USES none}
var CtrlWord : word;
//const FixNum : Extended = 0.0000000000000000001;
//      Zero : Longint = 0;
asm
  fstcw   word ptr [CtrlWord]        // Save control word
  fldcw   word ptr [CWChop]          // Set Rounding towards zero
  fld     [Value]       	     // ST(0)=Value

//  fcom   [Zero]
//  fstsw  ax
//  fld     [FixNum]
//  sahf
//  jae @@0
//  fchs
//@@0:
//  faddp   st(1),st(0)

  frndint
  fldcw   word ptr [CtrlWord]        // Restore previous control word
  // result is popped after return by VP
end;

function round64( const Value : Extended ) : tQuad; assembler; {$FRAME-}{$USES none}
//tquad := round(Value); // crash
var x: TQuad;
asm
  fld   [Value] 		     // ST(0)=Value
////  frndint   			 // assume default: rounding to nearest
  fistp qword ptr x
  fild qword ptr x
  // result is popped after return by VP
end;

function or64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  mov   eax, TQuadRec[Op2].Lo
  or    TQuadRec[Op1].Lo, eax
  mov   eax, TQuadRec[Op2].Hi
  or    TQuadRec[Op1].Hi, eax
  fild  [Op1]
end;

function and64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  mov   eax, TQuadRec[Op2].Lo
  and    TQuadRec[Op1].Lo, eax
  mov   eax, TQuadRec[Op2].Hi
  and    TQuadRec[Op1].Hi, eax
  fild  [Op1]
end;

function xor64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  mov   eax, TQuadRec[Op2].Lo
  xor   TQuadRec[Op1].Lo, eax
  mov   eax, TQuadRec[Op2].Hi
  xor   TQuadRec[Op1].Hi, eax
  fild  [Op1]
end;

function shl64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES ecx}
asm
  mov   ecx, TQuadRec[Op2].Lo // get the shift count
  cmp   cl, 63
  jg    @zero

  mov   eax, TQuadRec[Op1].Lo // get target operand

  cmp   cl, 31
  jg    @zerolo

  shld  TQuadRec[Op1].Hi, eax, cl
  shl   TQuadRec[Op1].Lo, cl

  jmp   @finish
@zerolo:
  sub   cl, 32
  shl   eax, cl
  mov   TQuadRec[Op1].Lo, 0
  mov   TQuadRec[Op1].Hi, eax
  jmp   @finish
@zero:
  fldz
  jmp   @exit
@finish:
  fild  [Op1]
@exit:
end;

function shr64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  mov   ecx, TQuadRec[Op2].Lo // get the shift count
  cmp   cl, 63
  jg    @zero

  mov   eax, TQuadRec[Op1].Hi // get target operand

  cmp   cl, 31
  jg    @zerohi

  shrd  TQuadRec[Op1].Lo, eax, cl
  shr   TQuadRec[Op1].Hi, cl

  jmp   @finish
@zerohi:
  sub   cl, 32
  shr   eax, cl
  mov   TQuadRec[Op1].Hi, 0
  mov   TQuadRec[Op1].Lo, eax
  jmp   @finish
@zero:
  fldz
  jmp   @exit
@finish:
  fild  [Op1]
@exit:
end;

function not64(const Op1: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  not   TQuadRec[Op1].Lo
  not   TQuadRec[Op1].Hi
  fild  [Op1]
end;

function neg64(const Op1: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  not   TQuadRec[Op1].Lo
  not   TQuadRec[Op1].Hi
  fild  [Op1]
  fld1
  faddp st(1), st(0) // stack cleanup
end;

function add64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm // clean up stack after operation -> addp
  fild  [Op2]
  fild  [Op1]
  faddp st(1),st(0)
end;

function sub64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm // clean up stack after operation -> subp
  fild  [Op1]
  fild  [Op2]
  fsubp st(1),st(0)
end;

function mul64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm // clean up stack after operation -> mulp
  fild  [Op2]
  fild  [Op1]
  fmulp st(1),st(0)
end;

function div64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm // clean up stack after operation -> divp
  fild  [Op1]
  fild  [Op2]
  fdivp st(1),st(0)
end;

function idiv64(const Op1, Op2: TQuad) : tQuad;
begin // integer division: whole numbers only
  Result := Trunc64(div64(Op1, Op2));
end;

function mod64(const Op1, Op2: TQuad) : tQuad; assembler; {$FRAME-}{$USES none}
asm
  fild  [Op2]
  fild  [Op1]
  fprem
  ffree st(1)   // clean up stack
end;

function cmp64(Op1: TQuad; Op2: TQuad; Condition: tCondition) : Boolean; assembler; {$FRAME-}{$USES edx}
asm
  fild  [Op2]
  fild  [Op1]
  fcompp
  fnstsw ax
  mov   edx, [Condition]
  shl   edx, 2  	 {*4}
  add   edx, [Condition] {*5}
  add   edx, Offset @Start
  sahf
  jmp   edx
@Start:
  seto  al     {\ 5    }
  jmp   @Done  {/ bytes}
  setno al
  jmp   @Done
  setb  al
  jmp   @Done
  setae al
  jmp   @Done
  sete  al
  jmp   @Done
  setne al
  jmp   @Done
  setbe al
  jmp   @Done
  seta  al
  jmp   @Done
  sets  al
  jmp   @Done
  setns al
  jmp   @Done
  setp  al
  jmp   @Done
  setnp al
  jmp   @Done
  setl  al
  jmp   @Done
  setge al
  jmp   @Done
  setle al
  jmp   @Done
  setg  al
@Done:
end;

function GetNValue64(L: tQuad; A: Integer) : shortstring;
var i,j: byte;
begin
  Result := '';
  if L = 0 then Result := '0' else
  while (L <> 0) do with tQuadRec(L) do
  begin
    Result := HexStr[(lo and ((1 shl A) - 1)) + 1] + Result;
    L := shr64(L, A);
  end;
end;

function q2hex(const Value: TQuad) : shortstring;
begin
  Result := GetNValue64(Value, 4);
end;

function q2bin(const Value: TQuad) : shortstring;
begin
  Result := GetNValue64(Value, 1);
end;

function q2oct(const Value: TQuad) : shortstring;
begin
  Result := GetNValue64(Value, 3);
end;

var FPU_CW : Word;

procedure Int64Init; assembler; {&USES None} {&FRAME-}
asm
  fstcw   word ptr [FPU_CW]           // Save control word
  fldcw   word ptr [CWDefault]        // Load CW
end;

procedure Int64Done; assembler; {&USES None} {&FRAME-}
asm
  fldcw   word ptr [FPU_CW]           // Load control word
end;

{ Ceil: Smallest integer >= X, |X| < MaxInt }
function Ceil64(X: Extended): Extended;
{&Frame-} {&Uses none}
asm
  fld     x
  fstcw   word ptr x      { save rounding control }
  fwait
  mov     ax, word ptr x
  and     ax, not mCW_RC  {Clear the rounding options}
  or      ax, RC_Up shl sCW_RC
  xchg    ax, word ptr x
  fldcw   word ptr x      { set rounding toward positive infinity }
  frndint       	  { round x }
  xchg    ax, word ptr x
  fldcw   word ptr x      { reset rounding control }
end;

{ Floor: Largest integer <= X, |X| < MaxInt }
function Floor64(X: Extended): Extended;
{&Frame-} {&Uses none}
asm
  fld     x
  fstcw   word ptr x      { save rounding control }
  fwait
  mov     ax, word ptr x
  and     ax, not mCW_RC  {Clear the rounding options}
  or      ax, RC_Down shl sCW_RC
  xchg    ax, word ptr x
  fldcw   word ptr x      { set rounding toward negative infinity }
  frndint       	  { round x }
  xchg    ax, word ptr x
  fldcw   word ptr x      { reset rounding control }
end;

end.