jclmath.pas

East make Tray Icon in delphi

  if (N and 1) = 0 then   //Zero or even
  begin
    Result := False;
    Exit;
  end;
  if PrimeSet = nil then // initialize look-up table
     InitPrimeSet;
  if N <= PrimeCacheLimit then // do look-up
    Result := PrimeSet.GetBit(N div 2)
  else
  begin // calculate
    R := N;
    Max := Round(Sqrt (R));
    if Max > PrimeCacheLimit then
    begin
      raise EJclMathError.CreateRes(@RsUnexpectedValue);
      Exit;
    end;
    I := 1;
    repeat
      Inc(I,2);
      if PrimeSet.GetBit(I div 2) then
        if N mod I = 0 then
        begin
          Result := False;
          Exit;
        end;
    until I >= Max;
    Result := True;
  end;
end;

{ Rabin-Miller Strong Primality Test }

function IsPrimeRM(N: Cardinal): Boolean;
asm
        TEST  EAX,1            // Odd(N) ??
        JNZ   @@1
        CMP   EAX,2            // N == 2 ??
        SETE  AL
        RET
@@1:    CMP   EAX,73
        JBE   @@C
        PUSH  ESI
        PUSH  EDI
        PUSH  EBX
        PUSH  EBP
        PUSH  EAX              // save N as Param for @@5
        LEA   EBP,[EAX - 1]    // M == N -1, Exponent
        MOV   ECX,32           // calc remaining Bits of M and shift M'
        MOV   ESI,EBP
@@2:    DEC   ECX
        SHL   ESI,1
        JNC   @@2
        PUSH  ECX              // save Bits as Param for @@5
        PUSH  ESI              // save M' as Param for @@5
        CMP   EAX,08A8D7Fh     // N >= 9080191 ??
        JAE   @@3
// now if (N < 9080191) and SPP(31, N) and SPP(73, N) then N is prime
        MOV   EAX,31
        CALL  @@5
        JC    @@4
        MOV   EAX,73
        PUSH  OFFSET @@4
        JMP   @@5
// now if (N < 4759123141) and SPP(2, N) and SPP(7, N) and SPP(61, N) then N is prime
@@3:    MOV   EAX,2
        CALL  @@5
        JC    @@4
        MOV   EAX,7
        CALL  @@5
        JC    @@4
        MOV   EAX,61
        CALL  @@5
@@4:    SETNC AL
        ADD   ESP,4 * 3
        POP   EBP
        POP   EBX
        POP   EDI
        POP   ESI
        RET
// do a Strong Pseudo Prime Test
@@5:    MOV   EBX,[ESP + 12]   // N on stack
        MOV   ECX,[ESP +  8]   // remaining Bits
        MOV   ESI,[ESP +  4]   // M'
        MOV   EDI,EAX          // T = b, temp. Base
@@6:    DEC   ECX
        MUL   EAX
        DIV   EBX
        MOV   EAX,EDX
        SHL   ESI,1
        JNC   @@7
        MUL   EDI
        DIV   EBX
        AND   ESI,ESI
        MOV   EAX,EDX
@@7:    JNZ   @@6
        CMP   EAX,1            // b^((N -1)(2^s)) mod N ==  1 mod N ??
        JE    @@A
@@8:    CMP   EAX,EBP          // b^((N -1)(2^s)) mod N == -1 mod N ??
        JE    @@A
        DEC   ECX              // second part to 2^s
        JNG   @@9
        MUL   EAX
        DIV   EBX
        CMP   EDX,1
        MOV   EAX,EDX
        JNE   @@8
@@9:    STC
@@A:    RET
@@B:    DB    3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73
@@C:    MOV   EDX,OFFSET @@B
        MOV   ECX,19
@@D:    CMP   AL,[EDX + ECX]
        JE    @@E
        DEC   ECX
        JNL   @@D
@@E:    SETE  AL
end;

function PrimeFactors(N: Cardinal): TDynCardinalArray;
var
  I, L, Max: Cardinal;
  R: Extended;
begin
  SetLength(Result, 0);
  if N <= 1 then
    Exit
  else
  begin
    if PrimeSet = nil then
      InitPrimeSet;
    L := 0;
    R := N;
    R := Sqrt(R);
    Max := Round(R);         // only one factor can be > Sqrt (N)
    if N mod 2 = 0 then      // test even at first
    begin                    // 2 is a prime factor
      Inc(L);
      SetLength(Result, L);
      Result[L - 1] := 2;
      repeat
        N := N div 2;
        if N = 1 then        // no more factors
          Exit;
      until N mod 2 <> 0;
    end;
    I := 3;                  // test all odd factors
    repeat
      if (N mod I = 0) and IsPrime(I)  then
      begin                  // I is a prime factor
        Inc(L);
        SetLength(Result, L);
        Result[L - 1] := I;
        repeat
          N := N div I;
          if N = 1 then      // no more factors
            Exit;
        until N mod I <> 0;
      end;
      Inc(I, 2);
    until I > Max;
    Inc(L);                  // final factor (> Sqrt(N))
    SetLength(Result, L);
    Result[L - 1] := N;
  end;
end;

function IsPrimeFactor(const F, N: Cardinal): Boolean;
begin
  Result := (N mod F = 0) and IsPrime(F);
end;

function IsRelativePrime(const X, Y: Cardinal): Boolean;
begin
  Result := GCD(X, Y) = 1;
end;

procedure SetPrimalityTest(const Method: TPrimalityTestMethod);
begin
  case Method of
    ptTrialDivision:
      IsPrime := IsPrimeTD;
    ptRabinMiller:
      IsPrime := IsPrimeRM;
  end;
end;

//=== Floating point value classification ====================================

const
  fpEmpty = TFloatingPointClass(Ord(High(TFloatingPointClass))+1);

  FPClasses: array [0..6] of TFloatingPointClass =
   (
    fpInvalid,
    fpNaN,
    fpNormal,
    fpInfinite,
    fpZero,
    fpEmpty,    // should not happen
    fpDenormal
   );

function _FPClass: TFloatingPointClass;
// In: ST(0) Value to examine
//     ECX   address of GOT (PIC only)
asm
        FXAM
        XOR     EDX, EDX
        FNSTSW  AX
        FFREE   ST(0)
        FINCSTP
        BT      EAX, 14 // C3
        RCL     EDX, 1
        BT      EAX, 10 // C2
        RCL     EDX, 1
        BT      EAX, 8  // C0
        RCL     EDX, 1
        {$IFDEF PIC}
        MOVZX   EAX, TFloatingPointClass([ECX].FPClasses[EDX])
        {$ELSE}
        MOVZX   EAX, TFloatingPointClass(FPClasses[EDX])
        {$ENDIF PIC}
end;

function FloatingPointClass(const Value: Single): TFloatingPointClass; overload;
asm
        {$IFDEF PIC}
        CALL    GetGOT
        MOV     ECX, EAX
        {$ENDIF PIC}
        FLD     Value
        CALL    _FPClass
end;

function FloatingPointClass(const Value: Double): TFloatingPointClass; overload;
asm
        {$IFDEF PIC}
        CALL    GetGOT
        MOV     ECX, EAX
        {$ENDIF PIC}
        FLD     Value
        CALL    _FPClass
end;

function FloatingPointClass(const Value: Extended): TFloatingPointClass; overload;
asm
        {$IFDEF PIC}
        CALL    GetGOT
        MOV     ECX, EAX
        {$ENDIF PIC}
        FLD     Value
        CALL    _FPClass
end;

//=== NaN and Infinity support ===============================================

function IsInfinite(const Value: Single): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpInfinite;
end;

function IsInfinite(const Value: Double): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpInfinite;
end;

function IsInfinite(const Value: Extended): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpInfinite;
end;

const
  sSignBit = 31;
  dSignBit = 63;
  xSignBit = 79;

type
  TSingleBits = set of 0..sSignBit;
  TDoubleBits = set of 0..dSignBit;
  TExtendedBits = set of 0..xSignBit;

  sFractionBits = 0..22; // Single type fraction bits
  dFractionBits = 0..51; // Double type fraction bits
  xFractionBits = 0..62; // Extended type fraction bits

  sExponentBits = 23..sSignBit-1;
  dExponentBits = 52..dSignBit-1;
  xExponentBits = 64..xSignBit-1;

  QWord = Int64;

  PExtendedRec = ^TExtendedRec;
  TExtendedRec = packed record
    Significand: QWord;
    Exponent: Word;
  end;

const
  ZeroTag = $3FFFFF;
  InvalidTag = TNaNTag($80000000);
  NaNTagMask = $3FFFFF;

  sNaNQuietFlag = High(sFractionBits);
  dNaNQuietFlag = High(dFractionBits);
  xNaNQuietFlag = High(xFractionBits);

  dNaNTagShift = High(dFractionBits) - High(sFractionBits);
  xNaNTagShift = High(xFractionBits) - High(sFractionBits);

  sNaNBits = $7F800000;
  dNaNBits = $7FF0000000000000;

  sQuietNaNBits = sNaNBits or (1 shl sNaNQuietFlag);
  dQuietNaNBits = dNaNBits or (Int64(1) shl dNaNQuietFlag);

function IsNaN(const Value: Single): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpNaN;
end;

function IsNaN(const Value: Double): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpNaN;
end;

function IsNaN(const Value: Extended): Boolean; overload;
begin
  Result := FloatingPointClass(Value) = fpNaN;
end;

procedure CheckNaN(const Value: Single); overload;
var
  SaveExMask: T8087Exceptions;
begin
  SaveExMask := Mask8087Exceptions([emInvalidOp]);
  try
    if FloatingPointClass(Value) <> fpNaN then
      raise EJclMathError.CreateRes(@RsNoNaN);
  finally
    SetMasked8087Exceptions(SaveExMask);
  end;
end;

procedure CheckNaN(const Value: Double); overload;
var
  SaveExMask: T8087Exceptions;
begin
  SaveExMask := Mask8087Exceptions([emInvalidOp]);
  try
    if FloatingPointClass(Value) <> fpNaN then
      raise EJclMathError.CreateRes(@RsNoNaN);
  finally
    SetMasked8087Exceptions(SaveExMask);
  end;
end;

procedure CheckNaN(const Value: Extended); overload;
var
  SaveExMask: T8087Exceptions;
begin
  SaveExMask := Mask8087Exceptions([emInvalidOp]);
  try
    if FloatingPointClass(Value) <> fpNaN then
      raise EJclMathError.CreateRes(@RsNoNaN);
  finally
    SetMasked8087Exceptions(SaveExMask);
  end;
end;

function GetNaNTag(const NaN: Single): TNaNTag;
var
  Temp: Integer;
begin
  CheckNaN(NaN);
  Temp := PLongint(@NaN)^ and NaNTagMask;
  if sSignBit in TSingleBits(NaN) then
    Result := -Temp
  else
  if Temp = ZeroTag then
    Result := 0
  else
    Result := Temp;
end;

function GetNaNTag(const NaN: Double): TNaNTag;
var
  Temp: Integer;
begin
  CheckNaN(NaN);
  Temp := (PInt64(@NaN)^ shr dNaNTagShift) and NaNTagMask;
  {$IFDEF FPC}
  if Int64(NaN) < 0 then
  {$ELSE}
  if dSignBit in TDoubleBits(NaN) then
  {$ENDIF FPC}
    Result := -Temp
  else
  if Temp = ZeroTag then
    Result := 0
  else
    Result := Temp;
end;

function GetNaNTag(const NaN: Extended): TNaNTag;
var
  Temp: Integer;
begin
  CheckNaN(NaN);
  Temp := (PExtendedRec(@NaN)^.Significand shr xNaNTagShift) and NaNTagMask;

  {$IFDEF FPC}
  if (TExtendedRec(NaN).Exponent and $8000) <> 0 then
  {$ELSE}
  if xSignBit in TExtendedBits(NaN) then
  {$ENDIF FPC}
    Result := -Temp
  else
  if Temp = ZeroTag then
    Result := 0
  else
    Result := Temp;
end;

{$IFDEF MSWINDOWS}

type
  TRealType = (rtUndef, rtSingle, rtDouble, rtExtended);

  { ExceptionInformation record for FPU exceptions under WinNT,
    where documented? }
  PFPUExceptionInfo = ^TFPUExceptionInfo;
  TFPUExceptionInfo = packed record
    Unknown: array [0..7] of Longint;
    ControlWord: Word;
    Dummy1: Word;
    StatusWord: Word;
    Dummy2: Word;
    TagWord: Word;
    Dummy3: Word;
    InstructionPtr: Pointer;
    UnknownW: Word;
    OpCode: Word;  // Note: 5 most significant bits of first opcode byte
                   // (always 11011b) not stored in FPU opcode register
    OperandPtr: Pointer;
    UnknownL: Longint;
  end;

  TExceptObjProc = function(P: PExceptionRecord): Exception;

var
  PrevExceptObjProc: TExceptObjProc;
  ExceptObjProcInitialized: Boolean = False;

function GetExceptionObject(P: PExceptionRecord): Exception;
var
  Tag: TNaNTag;
  FPUExceptInfo: PFPUExceptionInfo;
  OPtr: Pointer;
  OType: TRealType;

  function GetOperandType(OpCode: Word): TRealType;
  var
    NNN: 0..7;
  begin
    Result := rtUndef;
    NNN := (Lo(OpCode) shr 3) and 7;   // NNN field of ModR/M byte
    if Lo(OpCode) <= $BF then
    case Hi(OpCode) of   // 3 least significant bits of first opcode byte
      0:
        Result := rtSingle;
      1:
        if NNN < 4 then
          Result := rtSingle;
      // Extended signaling NaNs don't cause exceptions on FLD/FST(P) ?!
      3:
        if NNN = 5 then
          Result := rtExtended;
      4:
        Result := rtDouble;
      5:
        if NNN = 0 then
          Result := rtDouble;
    end;
  end;

begin
  Tag := InvalidTag; // shut up compiler warning
  OType := rtUndef;
  if P^.ExceptionCode = STATUS_FLOAT_INVALID_OPERATION then
  begin
    FPUExceptInfo := @P^.ExceptionInformation;
    OPtr := FPUExceptInfo^.OperandPtr;
    OType := GetOperandType(FPUExceptInfo^.OpCode);
    case OType of
      rtSingle:
        Tag := GetNaNTag(PSingle(OPtr)^);
      rtDouble:
        Tag := GetNaNTag(PDouble(OPtr)^);