fmath.pas

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

{ ----------------------------------------------------------------------
  Elementary functions
  ---------------------------------------------------------------------- }

  function Expo(X : Float) : Float; overload;
  begin
    MathError := FN_OK;
    if X < MINLOG then
      Expo := DefaultVal(FN_UNDERFLOW, 0.0)
    else if X > MAXLOG then
      Expo := DefaultVal(FN_OVERFLOW, MAXNUM)
    else
      Expo := Exp(X);
  end;

  function Exp2(X : Float) : Float;
  var
    XLn2 : Float;
  begin
    MathError := FN_OK;
    XLn2 := X * LN2;
    if XLn2 < MINLOG then
      Exp2 := DefaultVal(FN_UNDERFLOW, 0.0)
    else if XLn2 > MAXLOG then
      Exp2 := DefaultVal(FN_OVERFLOW, MAXNUM)
    else
      Exp2 := Exp(XLn2);
  end;

  function Exp10(X : Float) : Float;
  var
    XLn10 : Float;
  begin
    MathError := FN_OK;
    XLn10 := X * LN10;
    if XLn10 < MINLOG then
      Exp10 := DefaultVal(FN_UNDERFLOW, 0.0)
    else if XLn10 > MAXLOG then
      Exp10 := DefaultVal(FN_OVERFLOW, MAXNUM)
    else
      Exp10 := Exp(XLn10);
  end;

  function Log(X : Float) : Float; overload;
  begin
    MathError := FN_OK;
    if X < 0.0 then
      Log := DefaultVal(FN_DOMAIN, - MAXNUM)
    else if X = 0.0 then
      Log := DefaultVal(FN_SING, - MAXNUM)
    else
      Log := Ln(X);
  end;

  function Log10(X : Float) : Float;
  begin
    MathError := FN_OK;
    if X < 0.0 then
      Log10 := DefaultVal(FN_DOMAIN, - MAXNUM)
    else if X = 0.0 then
      Log10 := DefaultVal(FN_SING, - MAXNUM)
    else
      Log10 := Ln(X) * INVLN10;
  end;

  function Log2(X : Float) : Float;
  begin
    MathError := FN_OK;
    if X < 0.0 then
      Log2 := DefaultVal(FN_DOMAIN, - MAXNUM)
    else if X = 0.0 then
      Log2 := DefaultVal(FN_SING, - MAXNUM)
    else
      Log2 := Ln(X) * INVLN2;
  end;

  function LogA(X, A : Float) : Float;
  var
    Y : Float;
  begin
    Y := Log(X);
    if MathError = FN_OK then
      if A = 1.0 then
        Y := DefaultVal(FN_SING, Sgn(Y) * MAXNUM)
      else
        Y := Y / Log(A);
    LogA := Y;
  end;

{ ----------------------------------------------------------------------
  Power functions
  ---------------------------------------------------------------------- }

  function PowerTests(X, Y : Float; var Res : Float) : Boolean;
  { Tests the cases X=0, Y=0 and Y=1. Returns X^Y in Res }
  begin
    if X = 0.0 then
      begin
        PowerTests := True;
        if Y = 0.0 then       { 0^0 = lim  X^X = 1 }
          Res := 1.0          {       X->0         }
        else if Y > 0.0 then
          Res := 0.0          { 0^Y = 0 }
        else
          Res := DefaultVal(FN_SING, MAXNUM);
      end
    else if Y = 0.0 then
      begin
        Res := 1.0;           { X^0 = 1 }
        PowerTests := True;
      end
    else if Y = 1.0 then
      begin
        Res := X;             { X^1 = X }
        PowerTests := True;
      end
    else
      PowerTests := False;
  end;

  function IntPower(X : Float; N : Integer) : Float;
  { Computes X^N by repeated multiplications }
  const
    InverseMaxNum = 1.0 / MAXNUM;
  var
    T      : Float;
    M      : Integer;
    Invert : Boolean;
  begin
    if PowerTests(X, N, T) then
      begin
        IntPower := T;
        Exit;
      end;

    Invert := (N < 0);    { Test if inverting is needed }
    if 1.0 < Abs(X) then  { Test for 0 ..|x| .. 1 }
      begin
        X := 1.0 / X;
        Invert := not Invert;
      end;

    { Legendre's algorithm for
      minimizing the number of multiplications }
    T := 1.0; M := Abs(N);
    while 0 < M do
      begin
        if Odd(M) then T := T * X;
        X := Sqr(X);
        M := M div 2;
      end;

    if Invert then
      if Abs(T) < InverseMaxNum then  { Only here overflow }
        T := DefaultVal(FN_OVERFLOW, Sgn(T) * MAXNUM)
      else
        T := 1.0 / T;

    IntPower := T;
  end;

  function Power(X : Float; N : Integer) : Float; overload;
  begin
    Power := IntPower(X, N);
  end;

  function Power(X, Y : Float) : Float; overload;
  { Computes X^Y = Exp(Y * Ln(X)), for X > 0
    Resorts to IntPower if Y is integer }
  var
    Res  : Float;
    YLnX : Float;
  begin
    if PowerTests(X, Y, Res) then
      Power := Res
    else if (Abs(Y) < MaxInt) and (Trunc(Y) = Y) then  { Integer exponent }
      Power := IntPower(X, Trunc(Y))
    else if X <= 0.0 then
      Power := DefaultVal(FN_DOMAIN, 0.0)
    else
      begin
        YLnX := Y * Ln(X);
        if YLnX < MINLOG then
          Power := DefaultVal(FN_UNDERFLOW, 0.0)
        else if YLnX > MAXLOG then
          Power := DefaultVal(FN_OVERFLOW, MAXNUM)
        else
          Power := Exp(YLnX);
      end;
  end;

  function Pythag(X, Y : Float) : Float;
  { Computes Sqrt(X^2 + Y^2) without destructive underflow or overflow }
  var
    AbsX, AbsY : Float;
  begin
    MathError := FN_OK;
    AbsX := Abs(X);
    AbsY := Abs(Y);
    if AbsX > AbsY then
      Pythag := AbsX * Sqrt(1.0 + Sqr(AbsY / AbsX))
    else if AbsY = 0.0 then
      Pythag := 0.0
    else
      Pythag := AbsY * Sqrt(1.0 + Sqr(AbsX / AbsY));
  end;

{ ----------------------------------------------------------------------
  Lambert's function
  ---------------------------------------------------------------------- }

var
  NBITS : Integer;

var
  LC : record
         EM, EM9, C13, C23, EM2, D12, X0, X1,
         AN3, AN4, AN5, AN6, S21, S22, S23 : Float;
       end;

  function LambertW(X : Float; UBranch, Offset : Boolean) : Float;
  var
    I, NITER : Integer;
    AN2, DELX, ETA, RETA, T, TEMP, TEMP2, TS, WAPR, XX, ZL, ZN : Float;

  begin
    if Offset then
      begin
        DELX := X;
        if DELX < 0.0 then
          begin
            LambertW := DefaultVal(FN_DOMAIN, 0.0);
            Exit;
          end;
        XX := X + LC.EM;
      end
    else
      begin
        if X < LC.EM then
          begin
            LambertW := DefaultVal(FN_DOMAIN, 0.0);
            Exit;
          end;

        if X = LC.EM then
          begin
            LambertW := - 1.0;
            Exit;
          end;

        XX := X;
        DELX := XX - LC.EM;
      end;

    if UBranch then
      begin
        if Abs(XX) <= LC.X0 then
          begin
            LambertW := XX / (1.0 + XX / (1.0 + XX / (2.0 + XX / (0.6 + 0.34 * XX))));
            Exit;
          end;

        if XX <= LC.X1 then
          begin
            RETA := Sqrt(LC.D12 * DELX);
            LambertW := RETA / (1.0 + RETA / (3.0 + RETA / (RETA / (LC.AN4 +
                        RETA / (RETA * LC.AN6 + LC.AN5)) + LC.AN3))) - 1.0;
            Exit;
          end;

        if XX <= 20.0 then
          begin
            RETA := SQRT2 * Sqrt(1.0 - XX / LC.EM);
            AN2 := 4.612634277343749 * Sqrt(Sqrt(RETA + 1.09556884765625));
            WAPR := RETA / (1.0 + RETA / (3.0 + (LC.S21 * AN2 + LC.S22) * RETA / (LC.S23 * (AN2 + RETA)))) - 1.0;
          end
        else
          begin
            ZL := Ln(XX);
            WAPR := Ln(XX / Ln(Power(XX / ZL, Exp(- 1.124491989777808 / (0.4225028202459761 + ZL)))));
          end
      end
    else
      begin
        if XX >= 0.0 then
          begin
            LambertW := DefaultVal(FN_DOMAIN, 0.0);
            Exit;
          end;

        if XX <= LC.X1 then
          begin
            RETA := Sqrt(LC.D12 * DELX);
            LambertW := RETA / (RETA / (3.0 + RETA / (RETA / (LC.AN4 +
                        RETA / (RETA * LC.AN6 - LC.AN5)) - LC.AN3)) - 1.0) - 1.0;
            Exit;
          end;

        if XX <= LC.EM9 then
          begin
            ZL := Ln(- XX);
            T := - 1.0 - ZL;
            TS := Sqrt(T);
            WAPR := ZL - (2.0 * TS) / (SQRT2 + (LC.C13 - T / (2.7E2 + TS * 127.0471381349219)) * TS);
          end
        else
          begin
            ZL := Ln(- XX);
            ETA := 2.0 - LC.EM2 * XX;
            WAPR := Ln(XX / Ln(- XX / ((1.0 - 0.5043921323068457 * (ZL + 1.0)) * (Sqrt(ETA) + ETA / 3.0) + 1.0)));
          end
      end;

    if NBITS < 52 then NITER := 1 else NITER := 2;

    for I := 1 to NITER do
      begin
        ZN := Ln(XX / WAPR) - WAPR;
        TEMP := 1.0 + WAPR;
        TEMP2 := TEMP + LC.C23 * ZN;
        TEMP2 := 2.0 * TEMP * TEMP2;
        WAPR := WAPR * (1.0 + (ZN / TEMP) * (TEMP2 - ZN) / (TEMP2 - 2.0 * ZN));
      end;

    LambertW := WAPR;
  end;

{ ----------------------------------------------------------------------
  Trigonometric functions
  ---------------------------------------------------------------------- }

  function FixAngle(Theta : Float) : Float;
  begin
    MathError := FN_OK;
    while Theta > PI do
      Theta := Theta - TWOPI;
    while Theta <= - PI do
      Theta := Theta + TWOPI;
    FixAngle := Theta;
  end;

  function Tan(X : Float) : Float; overload;
  var
    SinX, CosX : Float;
  begin
    MathError := FN_OK;
    SinX := Sin(X);
    CosX := Cos(X);
    if CosX = 0.0 then
      Tan := DefaultVal(FN_SING, Sgn(SinX) * MAXNUM)
    else
      Tan := SinX / CosX;
  end;

  function ArcSin(X : Float) : Float;
  var
    A : Float;
  begin
    MathError := FN_OK;
    A := Abs(X);
    if A > 1.0 then
      ArcSin := DefaultVal(FN_DOMAIN, Sgn(X) * PIDIV2)
    else if A = 1.0 then
      ArcSin := Sgn(X) * PIDIV2
    else
      ArcSin := ArcTan(X / Sqrt(1.0 - Sqr(X)));
  end;

  function ArcCos(X : Float) : Float;
  begin
    MathError := FN_OK;
    if X < - 1.0 then
      ArcCos := DefaultVal(FN_DOMAIN, PI)
    else if X > 1.0 then
      ArcCos := DefaultVal(FN_DOMAIN, 0.0)
    else if X = 1.0 then
      ArcCos := 0.0
    else if X = - 1.0 then
      ArcCos := PI
    else
      ArcCos := PIDIV2 - ArcTan(X / Sqrt(1.0 - Sqr(X)));
  end;

  function ArcTan2(Y, X : Float) : Float;
  var
    Theta : Float;
  begin
    MathError := FN_OK;
    if X = 0.0 then
      if Y = 0.0 then
        ArcTan2 := 0.0
      else if Y > 0.0 then
        ArcTan2 := PIDIV2
      else
        ArcTan2 := - PIDIV2
    else
      begin
        { 4th/1st quadrant -PI/2..PI/2 }
        Theta := ArcTan(Y / X);

        { 2nd/3rd quadrants }
        if X < 0.0 then
          if Y >= 0.0 then
            Theta := Theta + PI   { 2nd quadrant:  PI/2..PI }
          else
            Theta := Theta - PI;  { 3rd quadrant: -PI..-PI/2 }
        ArcTan2 := Theta;
      end;
  end;

{ ----------------------------------------------------------------------
  Hyperbolic functions
  ---------------------------------------------------------------------- }

  function Sinh(X : Float) : Float;
  var
    ExpX : Float;
  begin
    MathError := FN_OK;
    if (X < MINLOG) or (X > MAXLOG) then
      Sinh := DefaultVal(FN_OVERFLOW, Sgn(X) * MAXNUM)
    else