fmath.pas

delphi矩阵运算、回归分析等数学计算

      begin
        ExpX := Exp(X);
        Sinh := 0.5 * (ExpX - 1.0 / ExpX);
      end;
  end;

  function Cosh(X : Float) : Float;
  var
    ExpX : Float;
  begin
    MathError := FN_OK;
    if (X < MINLOG) or (X > MAXLOG) then
      Cosh := DefaultVal(FN_OVERFLOW, MAXNUM)
    else
      begin
        ExpX := Exp(X);
        Cosh := 0.5 * (ExpX + 1.0 / ExpX);
      end;
  end;

  procedure SinhCosh(X : Float; var SinhX, CoshX : Float);
  var
    ExpX, ExpMinusX : Float;
  begin
    MathError := FN_OK;
    if (X < MINLOG) or (X > MAXLOG) then
      begin
        CoshX := DefaultVal(FN_OVERFLOW, MAXNUM);
        SinhX := Sgn(X) * CoshX;
      end
    else
      begin
        ExpX := Exp(X);
        ExpMinusX := 1.0 / ExpX;
        SinhX := 0.5 * (ExpX - ExpMinusX);
        CoshX := 0.5 * (ExpX + ExpMinusX);
      end;
  end;

  function Tanh(X : Float) : Float;
  var
    SinhX, CoshX : Float;
  begin
    SinhCosh(X, SinhX, CoshX);
    Tanh := SinhX / CoshX;
  end;

  function ArcSinh(X : Float) : Float;
  begin
    MathError := FN_OK;
    ArcSinh := Ln(X + Sqrt(Sqr(X) + 1.0));
  end;

  function ArcCosh(X : Float) : Float;
  begin
    MathError := FN_OK;
    if X < 1.0 then
      ArcCosh := DefaultVal(FN_DOMAIN, 0.0)
    else
      ArcCosh := Ln(X + Sqrt(Sqr(X) - 1.0));
  end;

  function ArcTanh(X : Float) : Float;
  begin
    MathError := FN_OK;
    if (X < - 1.0) or (X > 1.0) then
      ArcTanh := DefaultVal(FN_DOMAIN, Sgn(X) * MAXNUM)
    else if (X = - 1.0) or (X = 1.0) then
      ArcTanh := DefaultVal(FN_SING, Sgn(X) * MAXNUM)
    else
      ArcTanh := 0.5 * Ln((1.0 + X) / (1.0 - X));
  end;

{ ----------------------------------------------------------------------
  Complex functions
  ---------------------------------------------------------------------- }

  function Cmplx(X, Y : Float) : Complex;
  var
    Z : Complex;
  begin
    Z.X := X;
    Z.Y := Y;
    Cmplx := Z;
  end;

  function Polar(R, Theta : Float) : Complex;
  begin
    Polar := Cmplx(R * Cos(Theta), R * Sin(Theta));
  end;

  function CAbs(Z : Complex) : Float;
  begin
    CAbs := Pythag(Z.X, Z.Y);
  end;

  function CArg(Z : Complex) : Float;
  begin
    CArg := ArcTan2(Z.Y, Z.X);
  end;

  function CNeg(A : Complex) : Complex;
  begin
    CNeg := Cmplx(-A.X, -A.Y);
  end;

  function CConj(A : Complex) : Complex;
  begin
    CConj := Cmplx(A.X, -A.Y);
  end;

  function CAdd(A, B : Complex) : Complex;
  begin
    CAdd := Cmplx(A.X + B.X, A.Y + B.Y);
  end;

  function CSub(A, B : Complex) : Complex;
  begin
    CSub := Cmplx(A.X - B.X, A.Y - B.Y);
  end;

  function CMult(A, B : Complex) : Complex;
  begin
    CMult := Cmplx(A.X * B.X - A.Y * B.Y, A.X * B.Y + A.Y * B.X);
  end;

  function CDiv(A, B : Complex) : Complex;
  var
    Temp : Float;
  begin
    if (B.X = 0.0) and (B.Y = 0.0) then
      begin
        MathError := FN_OVERFLOW;
        CDiv := C_infinity;
        Exit;
      end;
    Temp := Sqr(B.X) + Sqr(B.Y);
    CDiv := Cmplx((A.X * B.X + A.Y * B.Y) / Temp,
                  (A.Y * B.X - A.X * B.Y) / Temp);
  end;

  function CRoot(Z : Complex; K, N : Integer) : Complex;
  { CRoot can calculate all 'N' roots of 'A' by varying 'K' from 0..N-1 }
  var
    R, Theta : Float;
  begin
    if (N <= 0) or (K < 0) or (K >= N) then
      begin
        MathError := FN_DOMAIN;
        CRoot := C_zero;
        Exit;
      end;
    R := CAbs(Z);
    Theta := CArg(Z);
    if R = 0.0 then
      Croot := C_zero
    else
      Croot := Polar(Power(R, 1.0 / N), FixAngle((Theta + K * TWOPI) / N));
  end;

  function CSqrt(Z : Complex) : Complex;
  var
    R, Theta : Float;
  begin
    R := CAbs(Z);
    Theta := CArg(Z);
    if R = 0.0 then
      CSqrt := C_zero
    else
      CSqrt := Polar(Sqrt(R), FixAngle(0.5 * Theta));
  end;

  function CCos(Z : Complex) : Complex;
  var
    SinX, CosX, SinhY, CoshY : Float;
  begin
    SinX := Sin(Z.X);
    CosX := Cos(Z.X);
    SinhCosh(Z.Y, SinhY, CoshY);  { Called here to set MathError }
    CCos := Cmplx(CosX * CoshY, - SinX * SinhY);
  end;

  function CSin(Z : Complex) : Complex;
  var
    SinX, CosX, SinhY, CoshY : Float;
  begin
    SinX := Sin(Z.X);
    CosX := Cos(Z.X);
    SinhCosh(Z.Y, SinhY, CoshY);  { Called here to set MathError }
    CSin := Cmplx(SinX * CoshY, CosX * SinhY);
  end;

  function CArcTan(Z : Complex) : Complex;
  var
    XX, Yp1, Ym1 : Float;
  begin
    if (Z.X = 0.0) and (Abs(Z.Y) = 1.0) then  { Z = +/- i }
      begin
        MathError := FN_SING;
        CArcTan := Cmplx(0.0, Sgn(Z.Y) * MAXNUM);
        Exit;
      end;
    XX := Sqr(Z.X);
    Yp1 := Z.Y + 1.0;
    Ym1 := Z.Y - 1.0;
    CArcTan := Cmplx(0.5 * (ArcTan2(Z.X, - Ym1) - ArcTan2(- Z.X, Yp1)),
                     0.25 * Log((XX + Sqr(Yp1)) / (XX + Sqr(Ym1))));
  end;

  function Expo(Z : Complex) : Complex; overload;
  var
    ExpX : Float;
  begin
    ExpX := Expo(Z.X);
    if MathError = FN_OK then
      Expo := Cmplx(ExpX * Cos(Z.Y), ExpX * Sin(Z.Y))
    else
      Expo := Cmplx(ExpX, 0.0);
  end;

  function Log(Z : Complex) : Complex; overload;
  var
    R, Theta, LnR : Float;
  begin
    R := CAbs(Z);
    Theta := CArg(Z);
    LnR := Log(R);
    if MathError = FN_OK then
      Log := Cmplx(LnR, Theta)
    else
      Log := Cmplx(- MAXNUM, 0.0);
  end;

  function Power(Z : Complex; N : Integer) : Complex; overload;
  var
    R, Theta : Float;
  begin
    R := CAbs(Z);
    Theta := CArg(Z);
    if R = 0.0 then
      if N = 0 then
        Power := C_one
      else if N > 0 then
        Power := C_zero
      else
        begin
          MathError := FN_SING;
          Power := C_infinity;
        end
    else
      Power := Polar(IntPower(R, N), FixAngle(N * Theta));
  end;

  function Power(Z : Complex; X : Float) : Complex; overload;
  var
    R, Theta : Float;
  begin
    R := CAbs(Z);
    Theta := CArg(Z);
    if R = 0.0 then
      if X = 0.0 then
        Power := C_one
      else if X > 0.0 then
        Power := C_zero
      else
        begin
          MathError := FN_SING;
          Power := C_infinity;
        end
    else
      Power := Polar(Power(R, X), FixAngle(X * Theta));
  end;

  function Power(A, B : Complex) : Complex; overload;
  begin
    if (A.X = 0.0) and (A.Y = 0.0) then
      if (B.X = 0.0) and (B.Y = 0.0) then
        Power := C_one                         { lim a^a = 1 as a -> 0 }
      else
        Power := C_zero                        { 0^b = 0, b > 0 }
    else
      Power := Expo(CMult(B, Log(A)));         { a^b = exp(b * ln(a)) }
  end;

  function Tan(Z : Complex) : Complex; overload;
  var
    X2, Y2, SinX2, CosX2, SinhY2, CoshY2, Temp : Float;
  begin
    X2 := 2.0 * Z.X;
    Y2 := 2.0 * Z.Y;
    SinX2 := Sin(X2);
    CosX2 := Cos(X2);
    SinhCosh(Y2, SinhY2, CoshY2);
    if MathError = FN_OK then
      Temp := CosX2 + CoshY2
    else
      Temp := CoshY2;
    if Temp <> 0.0 then
      Tan := Cmplx(SinX2 / Temp, SinhY2 / Temp)
    else
      begin                  { Z = Pi/2 + k*Pi }
        MathError := FN_SING;
        Tan := C_infinity;
      end;
  end;

  function ArcSin(Z : Complex) : Complex; overload;
  var
    Rp, Rm, S, T, X2, XX, YY : Float;
    B                        : Complex;
  begin
    B := Cmplx(Z.Y, - Z.X);  { Y - i*X }
    X2 := 2.0 * Z.X;
    XX := Sqr(Z.X);
    YY := Sqr(Z.Y);
    S := XX + YY + 1.0;
    Rp := 0.5 * Sqrt(S + X2);
    Rm := 0.5 * Sqrt(S - X2);
    T := Rp + Rm;
    ArcSin := Cmplx(ArcSin(Rp - Rm), Sgn(B) * Log(T + Sqrt(Sqr(T) - 1.0)));
  end;

  function ArcCos(Z : Complex) : Complex; overload;
  begin
    ArcCos := CSub(C_pi_div_2, ArcSin(Z));  { Pi/2 - ArcSin(Z) }
  end;

  function Sinh(Z : Complex) : Complex; overload;
  var
    SinhX, CoshX : Float;
  begin
    SinhCosh(Z.X, SinhX, CoshX);
    Sinh := Cmplx(SinhX * Cos(Z.Y), CoshX * Sin(Z.Y));
  end;

  function Cosh(Z : Complex) : Complex; overload;
  var
    SinhX, CoshX : Float;
  begin
    SinhCosh(Z.X, SinhX, CoshX);
    Cosh := Cmplx(CoshX * Cos(Z.Y), SinhX * Sin(Z.Y))
  end;

  function Tanh(Z : Complex) : Complex; overload;
  var
    X2, Y2, SinY2, CosY2, SinhX2, CoshX2, Temp : Float;
  begin
    X2 := 2.0 * Z.X;
    Y2 := 2.0 * Z.Y;
    SinY2 := Sin(Y2);
    CosY2 := Cos(Y2);
    SinhCosh(X2, SinhX2, CoshX2);
    if MathError = FN_OK then
      Temp := CoshX2 + CosY2
    else
      Temp := CoshX2;
    if Temp <> 0.0 then
      Tanh := Cmplx(SinhX2 / Temp, SinY2 / Temp)
    else
      begin                  { Z = i * (Pi/2 + k*Pi) }
        MathError := FN_SING;
        Tanh := Cmplx(0.0, MAXNUM);
      end;
  end;

  function ArcSinh(Z : Complex) : Complex; overload;
  { ArcSinh(Z) = -i*ArcSin(i*Z) }
  begin
    ArcSinh := Cneg(CMult(C_i, ArcSin(CMult(C_i, Z))));
  end;

  function ArcCosh(Z : Complex) : Complex; overload;
  { ArcCosh(Z) = CSgn(Y + i(1-X))*i*ArcCos(Z) where Z = X+iY }
  var
    Temp : Complex;
  begin
    Temp := CMult(C_i, ArcCos(Z));
    if Sgn(Cmplx(Z.Y, 1.0 - Z.X)) = -1 then
      Temp := CNeg(Temp);
    ArcCosh := Temp;
  end;

  function ArcTanh(Z : Complex) : Complex; overload;
  { ArcTanh(Z) = -i*ArcTan(i*Z) }
  begin
    if (Abs(Z.X) = 1.0) and (Z.Y = 0.0) then  { Z = +/- 1 }
      begin
        MathError := FN_SING;
        ArcTanh := Cmplx(Sgn(Z.X) * MAXNUM, 0.0);
      end
    else
      ArcTanh := CNeg(CMult(C_i, CArcTan(CMult(C_i, Z))));
  end;

begin
  MathError := FN_OK;
  SgnZeroIsOne := True;

  { Initialize constants for Lambert's function }
  NBITS := Round(- Log2(MACHEP));
  with LC do
    begin
      EM := - Exp(- 1.0);
      EM9 := - Exp(- 9.0);
      C13 := 1.0 / 3.0;
      C23 := 2.0 * C13;
      EM2 := 2.0 / EM;
      D12 := - EM2;
      X0 := Power(MACHEP, 1.0 / 6.0) * 0.5;
      X1 := (1.0 - 17.0 * Power(MACHEP, 2.0 / 7.0)) * EM;
      AN3 := 8.0 / 3.0;
      AN4 := 135.0 / 83.0;
      AN5 := 166.0 / 39.0;
      AN6 := 3167.0 / 3549.0;
      S21 := 2.0 * SQRT2 - 3.0;
      S22 := 4.0 - 3.0 * SQRT2;
      S23 := SQRT2 - 2.0;
    end;
end.