jclmath.pas

East make Tray Icon in delphi

          FSQRT
          FYL2X
          FWAIT
  end;

begin
  DomainCheck(Abs(X) >= 1.0);
  Result := FArcTanH(X);
end;

function CosH(X: Float): Float;
{$IFDEF PUREPASCAL}
begin
  Result := 0.5 * (Exp(X) + Exp(-X));
end;
{$ELSE}
const
  RoundDown: Word = $177F;
  OneHalf: Float = 0.5;
var
  ControlWW: Word;
asm
        {$IFDEF PIC}
        CALL    GetGOT
        {$ENDIF PIC}
        FLD     X    { TODO : Legal values for X? }
        FLDL2E
        FMULP   ST(1), ST
        FSTCW   ControlWW
        {$IFDEF PIC}
        FLDCW   [EAX].RoundDown
        {$ELSE}
        FLDCW   RoundDown
        {$ENDIF PIC}
        FLD     ST(0)
        FRNDINT
        FLDCW   ControlWW
        FXCH
        FSUB    ST, ST(1)
        F2XM1
        FLD1
        FADDP   ST(1), ST
        FSCALE
        FST     ST(1)
        FLD1
        FDIVRP  ST(1), ST
        FADDP   ST(1), ST
        {$IFDEF PIC}
        FLD     [EAX].OneHalf
        {$ELSE}
        FLD     OneHalf
        {$ENDIF PIC}
        FMULP   ST(1), ST
        FWAIT
end;
{$ENDIF PUREPASCAL}

function CotH(X: Float): Float;
begin
  Result := 1 / TanH(X);
end;

function CscH(X: Float): Float;
begin
  Result := Exp(X) - Exp(-X);
  DomainCheck(Result = 0.0);
  Result := 2.0 / Result;
end;

function SecH(X: Float): Float;
begin
  Result := Exp(X) + Exp(-X);
  DomainCheck(Result = 0.0);
  Result := 2.0 / Result;
end;

function SinH(X: Float): Float; assembler;
const
  RoundDown: Word = $177F;
  OneHalf: Float = 0.5;
var
  ControlWW: Word;
asm
        {$IFDEF PIC}
        CALL    GetGOT
        {$ENDIF PIC}
        FLD     X  { TODO : Legal values for X? }
        FLDL2E
        FMULP   ST(1), ST
        FSTCW   ControlWW
        {$IFDEF PIC}
        FLDCW   [EAX].RoundDown
        {$ELSE}
        FLDCW   RoundDown
        {$ENDIF PIC}
        FLD     ST(0)
        FRNDINT
        FLDCW   ControlWW
        FXCH
        FSUB    ST, ST(1)
        F2XM1
        FLD1
        FADDP   ST(1), ST
        FSCALE
        FST     ST(1)
        FLD1
        FDIVRP  ST(1), ST
        FSUBP   ST(1), ST
        {$IFDEF PIC}
        FLD     [EAX].OneHalf
        {$ELSE}
        FLD     OneHalf
        {$ENDIF PIC}
        FMULP   ST(1), ST
        FWAIT
end;

function TanH(X: Float): Float;
begin
  if X > MaxTanH then
    Result := 1.0
  else
  begin
    if X < -MaxTanH then
      Result := -1.0
    else
    begin
      Result := Exp(X);
      Result := Result * Result;
      Result := (Result - 1.0) / (Result + 1.0);
    end;
  end;
end;

//=== Coordinate conversion ==================================================

function DegMinSecToFloat(const Degs, Mins, Secs: Float): Float; // obsolete
begin
  Result := Degs + (Mins / 60.0) + (Secs / 3600.0);
end;

procedure FloatToDegMinSec(const X: Float; var Degs, Mins, Secs: Float); // obsolete
var
  Y: Float;
begin
  Degs := System.Int(X);
  Y := Frac(X) * 60;
  Mins := System.Int(Y);
  Secs := Frac(Y) * 60;
end;

//=== Exponential ============================================================

function Exp(const X: Float): Float;
begin
  {$IFDEF MATH_EXT_EXTREMEVALUES}
  if IsSpecialValue(X) then
  begin
    if IsNaN(X) or (X = Infinity) then
      Result := X
    else
      Result := 0;
    Exit;
  end;
  {$ENDIF MATH_EXT_EXTREMEVALUES}

  Result := System.Exp(X);
end;

function Power(const Base, Exponent: Float): Float;
var
  IsAnInteger, IsOdd: Boolean;
begin
  if (Exponent = 0.0) or (Base = 1.0) then
    Result := 1
  else
  if Base = 0.0 then
  begin
    if Exponent > 0.0 then
      Result := 0.0
    else
      {$IFDEF MATH_EXT_EXTREMEVALUES}
      Result := Infinity;
      {$ELSE}
      raise EJclMathError.CreateRes(@RsPowerInfinite);
      {$ENDIF MATH_EXT_EXTREMEVALUES}
  end
  else
  if Base > 0.0 then
    Result := Exp(Exponent * System.Ln(Base))
  else
  begin
    IsAnInteger := (Frac(Exponent) = 0.0);
    if IsAnInteger then
    begin
      Result := Exp(Exponent * System.Ln(Abs(Base)));
      IsOdd := Abs(Round(ModFloat(Exponent, 2))) = 1;
      if IsOdd then
        Result := -Result;
    end
    else
      raise EJclMathError.CreateRes(@RsPowerComplex);
  end;
end;

function PowerInt(const X: Float; N: Integer): Float;
var
  M: Integer;
  T: Float;
  Xc: Float;
begin
  if X = 0.0 then
  begin
    if N = 0 then
      Result := 1.0
    else
    if N > 0 then
      Result := 0.0
    else
      Result := MaxFloatingPoint;
    Exit;
  end;

  if N = 0 then
  begin
    Result := 1.0;
    Exit;
  end;

  // Legendre's algorithm for minimizing the number of multiplications
  T := 1.0;
  M := Abs(N);
  Xc := X;
  repeat
    if Odd(M) then
    begin
      Dec(M);
      T := T * Xc;
    end
    else
    begin
      M := M div 2;
      Xc := Sqr(Xc);
    end;
  until M = 0;

  if N > 0 then
    Result := T
  else
    Result := 1.0 / T;
end;

function TenToY(const Y: Float): Float;
begin
  if Y = 0.0 then
    Result := 1.0
  else
    Result := Exp(Y * Ln10);
end;

function TruncPower(const Base, Exponent: Float): Float;
begin
  if Base > 0 then
    Result := Power(Base, Exponent)
  else
    Result := 0;
end;

function TwoToY(const Y: Float): Float;
begin
  if Y = 0.0 then
    Result := 1.0
  else
    Result := Exp(Y * Ln2);
end;

//=== Floating point support routines ========================================

function IsFloatZero(const X: Float): Boolean;
begin
  Result := Abs(X) < PrecisionTolerance;
end;

function FloatsEqual(const X, Y: Float): Boolean;
begin
  try
    if Y = 0 then
      // catch exact equality
      Result := (X = Y) or (Abs(1 - Y/X ) <= PrecisionTolerance)
    else
      // catch exact equality
      Result := (X = Y) or (Abs(1 - X/Y ) <= PrecisionTolerance);
  except
    Result := False;  // catch real rare overflow e.g.  1.0e3000/1.0e-3000
  end
end;

function MaxFloat(const X, Y: Float): Float;
begin
  if X < Y then
    Result := Y
  else
    Result := X;
end;

function MinFloat(const X, Y: Float): Float;
begin
  if X > Y then
    Result := Y
  else
    Result := X;
end;

function ModFloat(const X, Y: Float): Float;
var
  Z: Float;
begin
  Result := X / Y;
  Z := System.Int(Result);
  if Frac(Result) < 0.0 then
    Z := Z - 1.0;
  Result := X - Y * Z;
end;

function RemainderFloat(const X, Y: Float): Float;
begin
  Result := X - System.Int(X / Y) * Y;
end;

procedure SwapFloats(var X, Y: Float);
var
  T: Float;
begin
  T := X;
  X := Y;
  Y := T;
end;

procedure CalcMachineEpsSingle;
var
  One: Single;
  T: Single;
begin
  One := 1.0;
  EpsSingle := One;
  repeat
    EpsSingle := 0.5 * EpsSingle;
    T := One + EpsSingle;
  until One = T;
  EpsSingle := 2.0 * EpsSingle;
  ThreeEpsSingle := 3.0 * EpsSingle;
end;

procedure CalcMachineEpsDouble;
var
  One: Double;
  T: Double;
begin
  One := 1.0;
  EpsDouble := One;
  repeat
    EpsDouble := 0.5 * EpsDouble;
    T := One + EpsDouble;
  until One = T;
  EpsDouble := 2.0 * EpsDouble;
  ThreeEpsDouble := 3.0 * EpsDouble;
end;

procedure CalcMachineEpsExtended;
var
  One: Extended;
  T: Extended;
begin
  One := 1.0;
  EpsExtended := One;
  repeat
    EpsExtended := 0.5 * EpsExtended;
    T := One + EpsExtended;
  until One = T;
  EpsExtended := 2.0 * EpsExtended;
  ThreeEpsExtended := 3.0 * EpsExtended;
end;

procedure CalcMachineEps;
begin
  {$IFDEF MATH_EXTENDED_PRECISION}
  CalcMachineEpsExtended;
  Epsilon := EpsExtended;
  ThreeEpsilon := ThreeEpsExtended;
  {$ENDIF MATH_EXTENDED_PRECISION}
  {$IFDEF MATH_DOUBLE_PRECISION}
  CalcMachineEpsDouble;
  Epsilon := EpsDouble;
  ThreeEpsilon := ThreeEpsDouble;
  {$ENDIF MATH_DOUBLE_PRECISION}
  {$IFDEF MATH_SINGLE_PRECISION}
  CalcMachineEpsSingle;
  Epsilon := EpsSingle;
  ThreeEpsilon := ThreeEpsSingle;
  {$ENDIF MATH_SINGLE_PRECISION}
end;

procedure SetPrecisionToleranceToEpsilon;
begin
  CalcMachineEps;
  PrecisionTolerance := Epsilon;
end;

function SetPrecisionTolerance(NewTolerance: Float): Float;
begin
  Result := PrecisionTolerance;
  PrecisionTolerance := NewTolerance;
end;

//=== Miscellaneous ==========================================================

function Ceiling(const X: Float): Integer;
begin
  Result := Integer(Trunc(X));
  if Frac(X) > 0 then
    Inc(Result);
end;

function CommercialRound(const X: Float): Int64;
begin
  Result := Trunc(X);
  if Frac(Abs(X)) >= 0.5 then
    Result := Result + Sgn(X);
end;

const
  PreCompFactsCount = 33; // all factorials that fit in a Single

  {$IFDEF MATH_SINGLE_PRECISION}
  PreCompFacts: array [0..PreCompFactsCount] of Float =
   (
    1.0,
    1.0,
    2.0,
    6.0,
    24.0,
    120.0,
    720.0,
    5040.0,
    40320.0,
    362880.0,
    3628800.0,
    39916800.0,
    479001600.0,
    6227020800.0,
    87178289152.0,
    1307674279936.0,
    20922788478976.0,
    355687414628352.0,
    6.4023735304192E15,
    1.21645096004223E17,
    2.43290202316367E18,
    5.10909408371697E19,
    1.12400072480601E21,
    2.58520174445945E22,
    6.20448454699065E23,
    1.55112110792462E25,
    4.03291499589617E26,
    1.08888704151327E28,
    3.04888371623715E29,
    8.8417630793192E30,
    2.65252889961724E32,
    8.22283968552752E33,
    2.63130869936881E35,
    8.68331850984666E36
   );
  {$ENDIF MATH_SINGLE_PRECISION}
  {$IFDEF MATH_DOUBLE_PRECISION}
  PreCompFacts: array [0..PreCompFactsCount] of Float =
   (
    1.0,
    1.0,
    2.0,
    6.0,
    24.0,
    120.0,
    720.0,
    5040.0,
    40320.0,
    362880.0,
    3628800.0,
    39916800.0,
    479001600.0,
    6227020800.0,
    87178291200.0,
    1307674368000.0,
    20922789888000.0,
    355687428096000.0,
    6.402373705728E15,
    1.21645100408832E17,
    2.43290200817664E18,
    5.10909421717094E19,
    1.12400072777761E21,
    2.5852016738885E22,
    6.20448401733239E23,