fmath.pas

Delphi 的数学控件

{ **********************************************************************
  *                            Unit FMATH.PAS                          *
  *                             Version 3.0d                           *
  *                      (c) J. Debord, August 2003                    *
  **********************************************************************
    This unit implements some mathematical functions in Delphi Pascal
  **********************************************************************
  Notes:

  1) The default real type is DOUBLE (8-byte real).
     Other types may be selected by defining the symbols:

       SINGLEREAL   (Single precision, 4 bytes)
       EXTENDEDREAL (Extended precision, 10 bytes)

  2) Error handling: The global variable MathError returns the error code
     from the last function evaluation. It must be checked immediately
     after a function call:

       Y := f(X);        (* f is one of the functions of the library *)
       if MathError = FN_OK then ...

     The possible error codes are:

     ---------------------------------------------
     Error code   Value  Significance
     ---------------------------------------------
     FN_OK          0    No error
     FN_DOMAIN     -1    Argument domain error
     FN_SING       -2    Positive singularity
     FN_OVERFLOW   -3    Overflow range error
     FN_UNDERFLOW  -4    Underflow range error
     FN_TLOSS      -5    Total loss of precision
     FN_PLOSS      -6    Partial loss of precision
     ---------------------------------------------

     When an error occurs, a default value is attributed to the function.

     The standard functions Exp and Ln have been redefined according to
     the above conventions as Expo and Log.

  **********************************************************************
  References:

  Complex functions: modified Pascal code from E. F. Glynn
                     (http://www.efg2.com/Lab/)

  Lambert's function: translated Fortran code from D. A. Barry et al.
                      (http://www.netlib.org/toms/743)

  Thanks to Volker Walter <[e-mail]>
  for improving the Power functions.
  ********************************************************************** }

unit fmath;

interface

{ ----------------------------------------------------------------------
  Floating point type (Default = Double)
  ---------------------------------------------------------------------- }

{$IFDEF SINGLEREAL}
  type Float = Single;
{$ELSE}
{$IFDEF EXTENDEDREAL}
  type Float = Extended;
{$ELSE}
  {$DEFINE DOUBLEREAL}
  type Float = Double;
{$ENDIF}
{$ENDIF}

{ ----------------------------------------------------------------------
  Mathematical constants
  ---------------------------------------------------------------------- }

const
  PI         = 3.14159265358979323846;  { Pi }
  LN2        = 0.69314718055994530942;  { Ln(2) }
  LN10       = 2.30258509299404568402;  { Ln(10) }
  LNPI       = 1.14472988584940017414;  { Ln(Pi) }
  INVLN2     = 1.44269504088896340736;  { 1/Ln(2) }
  INVLN10    = 0.43429448190325182765;  { 1/Ln(10) }
  TWOPI      = 6.28318530717958647693;  { 2*Pi }
  PIDIV2     = 1.57079632679489661923;  { Pi/2 }
  SQRTPI     = 1.77245385090551602730;  { Sqrt(Pi) }
  SQRT2PI    = 2.50662827463100050242;  { Sqrt(2*Pi) }
  INVSQRT2PI = 0.39894228040143267794;  { 1/Sqrt(2*Pi) }
  LNSQRT2PI  = 0.91893853320467274178;  { Ln(Sqrt(2*Pi)) }
  LN2PIDIV2  = 0.91893853320467274178;  { Ln(2*Pi)/2 }
  SQRT2      = 1.41421356237309504880;  { Sqrt(2) }
  SQRT2DIV2  = 0.70710678118654752440;  { Sqrt(2)/2 }
  GOLD       = 1.61803398874989484821;  { Golden Mean = (1 + Sqrt(5))/2 }
  CGOLD      = 0.38196601125010515179;  { 2 - GOLD }

{ ----------------------------------------------------------------------
  Machine-dependent constants
  ---------------------------------------------------------------------- }

{$IFDEF SINGLEREAL}
const
  MACHEP = 1.192093E-7;               { Floating point precision: 2^(-23) }
  MAXNUM = 3.402823E+38;              { Max. floating point number: 2^128 }
  MINNUM = 1.175495E-38;              { Min. floating point number: 2^(-126) }
  MAXLOG = 88.72283;                  { Max. argument for Exp = Ln(MAXNUM) }
  MINLOG = -87.33655;                 { Min. argument for Exp = Ln(MINNUM) }
  MAXFAC = 33;                        { Max. argument for Factorial }
  MAXGAM = 34.648;                    { Max. argument for Gamma }
  MAXLGM = 1.0383E+36;                { Max. argument for LnGamma }
{$ELSE}
{$IFDEF DOUBLEREAL}
const
  MACHEP = 2.220446049250313E-16;     { 2^(-52) }
  MAXNUM = 1.797693134862315E+308;    { 2^1024 }
  MINNUM = 2.225073858507202E-308;    { 2^(-1022) }
  MAXLOG = 709.7827128933840;
  MINLOG = -708.3964185322641;
  MAXFAC = 170;
  MAXGAM = 171.624376956302;
  MAXLGM = 2.556348E+305;
{$ELSE}
{$IFDEF EXTENDEDREAL}
const
  MACHEP = 1.08420217248550444E-19;   { 2^(-63) }
  MAXNUM = 1.18973149535723103E+4932; { 2^16384 }
  MINNUM = 3.36210314311209558E-4932; { 2^(-16382) }
  MAXLOG = 11356.5234062941439;
  MINLOG = - 11355.137111933024;
  MAXFAC = 1754;
  MAXGAM = 1755.455;
  MAXLGM = 1.04848146839019521E+4928;
{$ENDIF}
{$ENDIF}
{$ENDIF}

{ ----------------------------------------------------------------------
  Error codes for mathematical functions
  ---------------------------------------------------------------------- }

const
  FN_OK        =   0;  { No error }
  FN_DOMAIN    = - 1;  { Argument domain error }
  FN_SING      = - 2;  { Singularity }
  FN_OVERFLOW  = - 3;  { Overflow range error }
  FN_UNDERFLOW = - 4;  { Underflow range error }
  FN_TLOSS     = - 5;  { Total loss of precision }
  FN_PLOSS     = - 6;  { Partial loss of precision }

{ ----------------------------------------------------------------------
  Global variables
  ---------------------------------------------------------------------- }

var
  MathError    : Integer;  { Error code from the latest function evaluation }
  SgnZeroIsOne : Boolean;  { Indicates if Sgn(0) returns 1 (default) or 0 }

{ ----------------------------------------------------------------------
  Functional type
  ---------------------------------------------------------------------- }

type
  TFunc = function(X : Float) : Float;

{ ----------------------------------------------------------------------
  Complex numbers
  ---------------------------------------------------------------------- }

type
  Complex = record
    X, Y : Float;
  end;

{ ----------------------------------------------------------------------
  Complex constants
  ---------------------------------------------------------------------- }

const
  C_infinity : Complex = (X : MAXNUM; Y : 0.0);
  C_zero     : Complex = (X : 0.0;    Y : 0.0);
  C_one      : Complex = (X : 1.0;    Y : 0.0);
  C_i        : Complex = (X : 0.0;    Y : 1.0);
  C_pi       : Complex = (X : PI;     Y : 0.0);
  C_pi_div_2 : Complex = (X : PIDIV2; Y : 0.0);

{ ----------------------------------------------------------------------
  Min, max, sign and exchange
  ---------------------------------------------------------------------- }

function Min(X, Y : Integer) : Integer; overload;
function Max(X, Y : Integer) : Integer; overload;

function Min(X, Y : Float) : Float; overload;
function Max(X, Y : Float) : Float; overload;

function Sgn(X : Float)     : Integer; overload;  { Sign }
function Sgn(Z : Complex)   : Integer; overload;  { Complex sign }
function DSgn(A, B : Float) : Float;              { Sgn(B) * |A| }

procedure Swap(var X, Y : Float);   overload;   { Exchange 2 reals }
procedure Swap(var X, Y : Integer); overload;   { Exchange 2 integers }
procedure Swap(var W, Z : Complex); overload;   { Exchange 2 complexes }

{ ----------------------------------------------------------------------
  Complex numbers
  ---------------------------------------------------------------------- }

function Cmplx(X, Y : Float)     : Complex;  { X + i.Y }
function Polar(R, Theta : Float) : Complex;  { R.exp(i.Theta) }

function CAbs(Z : Complex) : Float;          { |Z| }
function CArg(Z : Complex) : Float;          { Arg(Z) }

function CNeg(A : Complex)     : Complex;    { -A }
function CConj(A : Complex)    : Complex;    { A* }
function CAdd(A, B : Complex)  : Complex;    { A + B }
function CSub(A, B : Complex)  : Complex;    { A - B }
function CMult(A, B : Complex) : Complex;    { A * B }
function CDiv(A, B : Complex)  : Complex;    { A / B }

function CSqrt(Z : Complex)                 : Complex;  { Sqrt(Z) }
function CRoot(Z : Complex; K, N : Integer) : Complex;  { Z^(1/N), K=0..N-1 }

function CSin(Z : Complex)    : Complex;     { Sin(Z) }
function CCos(Z : Complex)    : Complex;     { Cos(Z) }
function CArcTan(Z : Complex) : Complex;     { ArcTan(Z) }

{ ----------------------------------------------------------------------
  Logarithms, exponentials, and powers (real or complex argument)
  ---------------------------------------------------------------------- }

function Expo(X : Float)   : Float;   overload;  { Exponential }
function Expo(Z : Complex) : Complex; overload;

function Log(X : Float)   : Float;   overload;   { Natural log }
function Log(Z : Complex) : Complex; overload;

function Power(X : Float; N : Integer)   : Float;   overload; { X^N }
function Power(X, Y : Float)             : Float;   overload; { X^Y, X >= 0 }
function Power(Z : Complex; N : Integer) : Complex; overload; { Z^N }
function Power(Z : Complex; X : Float)   : Complex; overload; { Z^X }
function Power(A, B : Complex)           : Complex; overload; { A^B }

{ ----------------------------------------------------------------------
  Other logarithms and exponentials (real argument only)
  ---------------------------------------------------------------------- }

function Exp2(X : Float)    : Float;  { 2^X }
function Exp10(X : Float)   : Float;  { 10^X }
function Log2(X : Float)    : Float;  { Log, base 2 }
function Log10(X : Float)   : Float;  { Decimal log }
function LogA(X, A : Float) : Float;  { Log, base A }

{ ----------------------------------------------------------------------
  Lambert's function
  ---------------------------------------------------------------------- }
function LambertW(X : Float; UBranch, Offset : Boolean) : Float;
{ ----------------------------------------------------------------------
  Lambert's W function: Y = W(X) ==> X = Y * Exp(Y)    X >= -1/e
  ----------------------------------------------------------------------
  X       = Argument
  UBranch = TRUE  for computing the upper branch (X >= -1/e, W(X) >= -1)
            FALSE for computing the lower branch (-1/e <= X < 0, W(X) <= -1)
  Offset  = TRUE  for computing W(X - 1/e), X >= 0
            FALSE for computing W(X)
  ---------------------------------------------------------------------- }

{ ----------------------------------------------------------------------
  Trigonometric functions (real or complex argument)
  ---------------------------------------------------------------------- }

function Tan(X : Float)   : Float;   overload;     { Tangent }
function Tan(Z : Complex) : Complex; overload;

function ArcSin(X : Float)   : Float;   overload;  { Arc sinus }
function ArcSin(Z : Complex) : Complex; overload;

function ArcCos(X : Float)   : Float;   overload;  { Arc cosinus }
function ArcCos(Z : Complex) : Complex; overload;

{ ----------------------------------------------------------------------
  Other trigonometric functions (real argument only)
  ---------------------------------------------------------------------- }

function Pythag(X, Y : Float)    : Float;  { Sqrt(X^2 + Y^2) }
function ArcTan2(Y, X : Float)   : Float;  { Angle (Ox, OM) with M(X,Y) }
function FixAngle(Theta : Float) : Float;  { Set Theta in -Pi..Pi }

{ ----------------------------------------------------------------------
  Hyperbolic functions (real or complex argument)
  ---------------------------------------------------------------------- }

function Sinh(X : Float)   : Float;   overload;     { Hyperbolic sine }
function Sinh(Z : Complex) : Complex; overload;

function Cosh(X : Float)   : Float;   overload;     { Hyperbolic cosine }
function Cosh(Z : Complex) : Complex; overload;

function Tanh(X : Float)   : Float;   overload;     { Hyperbolic tangent }
function Tanh(Z : Complex) : Complex; overload;

function ArcSinh(X : Float)   : Float;   overload;  { Inverse hyperbolic sine }
function ArcSinh(Z : Complex) : Complex; overload;

function ArcCosh(X : Float)   : Float;   overload;  { Inverse hyperbolic cosine }
function ArcCosh(Z : Complex) : Complex; overload;

function ArcTanh(X : Float)   : Float;   overload;  { Inverse hyperbolic tangent }
function ArcTanh(Z : Complex) : Complex; overload;

{ ----------------------------------------------------------------------
  Other hyperbolic function (real argument only)
  ---------------------------------------------------------------------- }

procedure SinhCosh(X : Float; var SinhX, CoshX : Float);  { Sinh & Cosh }

{ ********************************************************************** }

implementation

{ ----------------------------------------------------------------------
  Error handling function
  ---------------------------------------------------------------------- }

  function DefaultVal(ErrCode : Integer; Default : Float) : Float;
  { Sets the global variable MathError and returns
    the default value according to the error code }
  begin
    MathError := ErrCode;
    DefaultVal := Default;
  end;

{ ----------------------------------------------------------------------
  Min, max, sign and exchange
  ---------------------------------------------------------------------- }

  function Min(X, Y : Integer) : Integer; overload;
  begin
    if X < Y then Min := X else Min := Y;
  end;

  function Max(X, Y : Integer) : Integer; overload;
  begin
    if X > Y then Max := X else Max := Y;
  end;

  function Min(X, Y : Float) : Float; overload;
    begin
    if X < Y then Min := X else Min := Y;
  end;

  function Max(X, Y : Float) : Float; overload;
  begin
    if X > Y then Max := X else Max := Y;
  end;

  function Sgn0(X : Float) : Integer;
  begin
    if X > 0.0 then
      Sgn0 := 1
    else if X = 0.0 then
      Sgn0 := 0
    else
      Sgn0 := - 1;
  end;

  function Sgn(X : Float) : Integer;
  begin
    if X > 0.0 then
      Sgn := 1
    else if X < 0.0 then
      Sgn := - 1
    else
      Sgn := Ord(SgnZeroIsOne);
  end;

  function Sgn(Z : Complex) : Integer; overload;
  begin
    if Z.X > 0.0 then
      Sgn := 1
    else if Z.X < 0.0 then
      Sgn := - 1
    else if Z.Y > 0.0 then
      Sgn := 1
    else if Z.Y < 0.0 then
      Sgn := - 1
    else
      Sgn := Ord(SgnZeroIsOne);
  end;

  function DSgn(A, B : Float) : Float;
  begin
    if B < 0.0 then DSgn := - Abs(A) else DSgn := Abs(A);
  end;

  procedure Swap(var X, Y : Integer); overload;
  var
    Temp : Integer;
  begin
    Temp := X;
    X := Y;
    Y := Temp;
  end;

  procedure Swap(var X, Y : Float); overload;
  var
    Temp : Float;
  begin
    Temp := X;
    X := Y;
    Y := Temp;
  end;

  procedure Swap(var W, Z : Complex); overload;
  var
    Temp : Complex;
  begin
    Temp := W;
    W := Z;
    Z := Temp;
  end;