kolmath.pas | searchcode

/kolmath.pas

http://github.com/rofl0r/KOL
Pascal | 1631 lines | 1082 code | 186 blank | 363 comment | 102 complexity | dd033e8c48025758ec25a54e46fbd3fb MD5 | raw file

{=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-



        KKKKK    KKKKK    OOOOOOOOO    LLLLL

        KKKKK    KKKKK  OOOOOOOOOOOOO  LLLLL

        KKKKK    KKKKK  OOOOO   OOOOO  LLLLL

        KKKKK  KKKKK    OOOOO   OOOOO  LLLLL

        KKKKKKKKKK      OOOOO   OOOOO  LLLLL

        KKKKK  KKKKK    OOOOO   OOOOO  LLLLL

        KKKKK    KKKKK  OOOOO   OOOOO  LLLLL

        KKKKK    KKKKK  OOOOOOOOOOOOO  LLLLLLLLLLLLL

        KKKKK    KKKKK    OOOOOOOOO    LLLLLLLLLLLLL



  Key Objects Library (C) 2000 by Kladov Vladimir.



  mailto: bonanzas@xcl.cjb.net

  Home: http://kol.nm.ru

        http://xcl.cjb.net

        http://xcl.nm.ru



 =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-}

{

  This code is grabbed from standard math.pas unit,

  provided by Borland Delphi. This unit is for working with

  engineering (mathematical) functions. The main difference

  is that err unit specially designed to handle exceptions

  for KOL is used instead of SysUtils. This allows to make

  size of the executable smaller for about 5K. though this

  value is insignificant for project made with VCL, it can

  be more than 15% of executable file size made with KOL.  

}



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

{                                                       }

{       Borland Delphi Runtime Library                  }

{       Math Unit                                       }

{                                                       }

{       Copyright (C) 1996,99 Inprise Corporation       }

{                                                       }

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



unit kolmath;



{ This unit contains high-performance arithmetic, trigonometric, logorithmic,

  statistical and financial calculation routines which supplement the math

  routines that are part of the Delphi language or System unit. }



{$N+,S-}



{$I KOLDEF.INC}



interface



uses err, kol;



const   { Ranges of the IEEE floating point types, including denormals }

  MinSingle        =     1.5e-45;

  MaxSingle        =     3.4e+38;

  MinDouble        =     5.0e-324;

  MaxDouble        =     1.7e+308;

  MinExtended      =     3.4e-4932;

  MaxExtended      =     1.1e+4932;

  MinComp          =     -9.223372036854775807e+18;

  MaxComp          =     9.223372036854775807e+18;



{-----------------------------------------------------------------------

References:



1) P.J. Plauger, "The Standard C Library", Prentice-Hall, 1992, Ch. 7.

2) W.J. Cody, Jr., and W. Waite, "Software Manual For the Elementary

   Functions", Prentice-Hall, 1980.

3) Namir Shammas, "C/C++ Mathematical Algorithms for Scientists and Engineers",

   McGraw-Hill, 1995, Ch 8.

4) H.T. Lau, "A Numerical Library in C for Scientists and Engineers",

   CRC Press, 1994, Ch. 6.

5) "Pentium(tm) Processor User's Manual, Volume 3: Architecture

   and Programming Manual", Intel, 1994



All angle parameters and results of trig functions are in radians.



Most of the following trig and log routines map directly to Intel 80387 FPU

floating point machine instructions.  Input domains, output ranges, and

error handling are determined largely by the FPU hardware.

Routines coded in assembler favor the Pentium FPU pipeline architecture.

-----------------------------------------------------------------------}



function EAbs( D: Double ): Double;

function EMax( const Values: array of Double ): Double;

function EMin( const Values: array of Double ): Double;

function iMax( const Values: array of Integer ): Integer;

function iMin( const Values: array of Integer ): Integer;



{ Trigonometric functions }

function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }

function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }



{ ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.

  IN: |Y| < 2^64, |X| < 2^64, X <> 0   OUT: [-PI..PI] radians }

function ArcTan2(Y, X: Extended): Extended;



{ SinCos is 2x faster than calling Sin and Cos separately for the same angle }

procedure SinCos(Theta: Extended; var Sin, Cos: Extended) register;

function Tan(X: Extended): Extended;

function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }

function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }



{ Angle unit conversion routines }

function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}

function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }

function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }

function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }

function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }

function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }



{ Hyperbolic functions and inverses }

function Cosh(X: Extended): Extended;

function Sinh(X: Extended): Extended;

function Tanh(X: Extended): Extended;

function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }

function ArcSinh(X: Extended): Extended;

function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }



{ Logorithmic functions }

function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }

function Log10(X: Extended): Extended;                     { Log base 10 of X}

function Log2(X: Extended): Extended;                      { Log base 2 of X }

function LogN(Base, X: Extended): Extended;                { Log base N of X }



{ Exponential functions }



{ IntPower: Raise base to an integral power.  Fast. }

//function IntPower(Base: Extended; Exponent: Integer): Extended register;

// -- already defined in kol.pas



{ Power: Raise base to any power.

  For fractional exponents, or |exponents| > MaxInt, base must be > 0. }

function Power(Base, Exponent: Extended): Extended;





{ Miscellaneous Routines }



{ Frexp:  Separates the mantissa and exponent of X. }

procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer) register;



{ Ldexp: returns X*2**P }

function Ldexp(X: Extended; P: Integer): Extended register;



{ Ceil: Smallest integer >= X, |X| < MaxInt }

function Ceil(X: Extended):Integer;



{ Floor: Largest integer <= X,  |X| < MaxInt }

function Floor(X: Extended): Integer;



{ Poly: Evaluates a uniform polynomial of one variable at value X.

    The coefficients are ordered in increasing powers of X:

    Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }

function Poly(X: Extended; const Coefficients: array of Double): Extended;



{-----------------------------------------------------------------------

Statistical functions.



Common commercial spreadsheet macro names for these statistical and

financial functions are given in the comments preceding each function.

-----------------------------------------------------------------------}



{ Mean:  Arithmetic average of values.  (AVG):  SUM / N }

function Mean(const Data: array of Double): Extended;



{ Sum: Sum of values.  (SUM) }

function Sum(const Data: array of Double): Extended register;

function SumInt(const Data: array of Integer): Integer register;

function SumOfSquares(const Data: array of Double): Extended;

procedure SumsAndSquares(const Data: array of Double;

  var Sum, SumOfSquares: Extended) register;



{ MinValue: Returns the smallest signed value in the data array (MIN) }

function MinValue(const Data: array of Double): Double;

function MinIntValue(const Data: array of Integer): Integer;



function Min(A,B: Integer): Integer;

{$IFDEF _D4orHigher}

overload;

function Min(A,B: I64): I64; overload;

function Min(A,B: Single): Single; overload;

function Min(A,B: Double): Double; overload;

function Min(A,B: Extended): Extended; overload;

{$ENDIF}



{ MaxValue: Returns the largest signed value in the data array (MAX) }

function MaxValue(const Data: array of Double): Double;

function MaxIntValue(const Data: array of Integer): Integer;



function Max(A,B: Integer): Integer;

{$IFDEF _D4orHigher}

overload;

function Max(A,B: I64): I64; overload;

function Max(A,B: Single): Single; overload;

function Max(A,B: Double): Double; overload;

function Max(A,B: Extended): Extended; overload;

{$ENDIF}



{ Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }

function StdDev(const Data: array of Double): Extended;



{ MeanAndStdDev calculates Mean and StdDev in one call. }

procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);



{ Population Standard Deviation (STDP): Sqrt(PopnVariance).

  Used in some business and financial calculations. }

function PopnStdDev(const Data: array of Double): Extended;



{ Variance (VARS): TotalVariance / (N-1). aka Sample Variance }

function Variance(const Data: array of Double): Extended;



{ Population Variance (VAR or VARP): TotalVariance/ N }

function PopnVariance(const Data: array of Double): Extended;



{ Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }

function TotalVariance(const Data: array of Double): Extended;



{ Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }

function Norm(const Data: array of Double): Extended;



{ MomentSkewKurtosis: Calculates the core factors of statistical analysis:

  the first four moments plus the coefficients of skewness and kurtosis.

  M1 is the Mean.  M2 is the Variance.

  Skew reflects symmetry of distribution: M3 / (M2**(3/2))

  Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }

procedure MomentSkewKurtosis(const Data: array of Double;

  var M1, M2, M3, M4, Skew, Kurtosis: Extended);



{ RandG produces random numbers with Gaussian distribution about the mean.

  Useful for simulating data with sampling errors. }

function RandG(Mean, StdDev: Extended): Extended;



{-----------------------------------------------------------------------

Financial functions.  Standard set from Quattro Pro.



Parameter conventions:



From the point of view of A, amounts received by A are positive and

amounts disbursed by A are negative (e.g. a borrower's loan repayments

are regarded by the borrower as negative).



Interest rates are per payment period.  11% annual percentage rate on a

loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667



-----------------------------------------------------------------------}



type

  TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);



{ Double Declining Balance (DDB) }

function DoubleDecliningBalance(Cost, Salvage: Extended;

  Life, Period: Integer): Extended;



{ Future Value (FVAL) }

function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:

  Extended; PaymentTime: TPaymentTime): Extended;



{ Interest Payment (IPAYMT)  }

function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,

  FutureValue: Extended; PaymentTime: TPaymentTime): Extended;



{ Interest Rate (IRATE) }

function InterestRate(NPeriods: Integer;

  Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;



{ Internal Rate of Return. (IRR) Needs array of cash flows. }

function InternalRateOfReturn(Guess: Extended;

  const CashFlows: array of Double): Extended;



{ Number of Periods (NPER) }

function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;

  PaymentTime: TPaymentTime): Extended;



{ Net Present Value. (NPV) Needs array of cash flows. }

function NetPresentValue(Rate: Extended; const CashFlows: array of Double;

  PaymentTime: TPaymentTime): Extended;



{ Payment (PAYMT) }

function Payment(Rate: Extended; NPeriods: Integer;

  PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;



{ Period Payment (PPAYMT) }

function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;

  PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;



{ Present Value (PVAL) }

function PresentValue(Rate: Extended; NPeriods: Integer;

  Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;



{ Straight Line depreciation (SLN) }

function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;



{ Sum-of-Years-Digits depreciation (SYD) }

function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;



{type

  EInvalidArgument = class(EMathError) end;}



implementation



{$IFNDEF _D2orD3}

uses SysConst;

{$ENDIF}



function EAbs( D: Double ): Double;

begin

  Result := D;

  if Result < 0.0 then

    Result := -Result;

end;



function EMax( const Values: array of Double ): Double;

var I: Integer;

begin

  Result := Values[ 0 ];

  for I := 1 to High( Values ) do

    if Result < Values[ I ] then Result := Values[ I ];

end;



function EMin( const Values: array of Double ): Double;

var I: Integer;

begin

  Result := Values[ 0 ];

  for I := 1 to High( Values ) do

    if Result > Values[ I ] then Result := Values[ I ];

end;



function iMax( const Values: array of Integer ): Integer;

var I: Integer;

begin

  Result := Values[ 0 ];

  for I := 1 to High( Values ) do

    if Result < Values[ I ] then Result := Values[ I ];

end;



function iMin( const Values: array of Integer ): Integer;

var I: Integer;

begin

  Result := Values[ 0 ];

  for I := 1 to High( Values ) do

    if Result > Values[ I ] then Result := Values[ I ];

end;



function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;

  var CompoundRN: Extended): Extended; Forward;

function Compound(R: Extended; N: Integer): Extended; Forward;

function RelSmall(X, Y: Extended): Boolean; Forward;



type

  TPoly = record

    Neg, Pos, DNeg, DPos: Extended

  end;



const

  MaxIterations = 15;



procedure ArgError(const Msg: string);

begin

  raise Exception.Create(e_Math_InvalidArgument, Msg);

end;



function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180 }

begin

  Result := Degrees * (PI / 180);

end;



function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }

begin

  Result := Radians * (180 / PI);

end;



function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }

begin

  Result := Grads * (PI / 200);

end;



function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI}

begin

  Result := Radians * (200 / PI);

end;



function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }

begin

  Result := Cycles * (2 * PI);

end;



function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }

begin

  Result := Radians / (2 * PI);

end;



function LnXP1(X: Extended): Extended;

{ Return ln(1 + X).  Accurate for X near 0. }

asm

        FLDLN2

        MOV     AX,WORD PTR X+8               { exponent }

        FLD     X

        CMP     AX,$3FFD                      { .4225 }

        JB      @@1

        FLD1

        FADD

        FYL2X

        JMP     @@2

@@1:

        FYL2XP1

@@2:

        FWAIT

end;



{ Invariant: Y >= 0 & Result*X**Y = X**I.  Init Y = I and Result = 1. }

{function IntPower(X: Extended; I: Integer): Extended;

var

  Y: Integer;

begin

  Y := Abs(I);

  Result := 1.0;

  while Y > 0 do begin

    while not Odd(Y) do

    begin

      Y := Y shr 1;

      X := X * X

    end;

    Dec(Y);

    Result := Result * X

  end;

  if I < 0 then Result := 1.0 / Result

end;

}

(* -- already defined in kol.pas

function IntPower(Base: Extended; Exponent: Integer): Extended;

asm

        mov     ecx, eax

        cdq

        fld1                      { Result := 1 }

        xor     eax, edx

        sub     eax, edx          { eax := Abs(Exponent) }

        jz      @@3

        fld     Base

        jmp     @@2

@@1:    fmul    ST, ST            { X := Base * Base }

@@2:    shr     eax,1

        jnc     @@1

        fmul    ST(1),ST          { Result := Result * X }

        jnz     @@1

        fstp    st                { pop X from FPU stack }

        cmp     ecx, 0

        jge     @@3

        fld1

        fdivrp                    { Result := 1 / Result }

@@3:

        fwait

end;

*)



function Compound(R: Extended; N: Integer): Extended;

{ Return (1 + R)**N. }

begin

  Result := IntPower(1.0 + R, N)

end;



function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;

  var CompoundRN: Extended): Extended;

{ Set CompoundRN to Compound(R, N),

	return (1+Rate*PaymentTime)*(Compound(R,N)-1)/R;

}

begin

  if R = 0.0 then

  begin

    CompoundRN := 1.0;

    Result := N;

  end

  else

  begin

    { 6.1E-5 approx= 2**-14 }

    if EAbs(R) < 6.1E-5 then

    begin

      CompoundRN := Exp(N * LnXP1(R));

      Result := N*(1+(N-1)*R/2);

    end

    else

    begin

      CompoundRN := Compound(R, N);

      Result := (CompoundRN-1) / R

    end;

    if PaymentTime = ptStartOfPeriod then

      Result := Result * (1 + R);

  end;

end; {Annuity2}





procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);

{ Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.

  Accumulate positive and negative terms separately. }

var

  I: Integer;

  Neg, Pos, DNeg, DPos: Extended;

begin

  Neg := 0.0;

  Pos := 0.0;

  DNeg := 0.0;

  DPos := 0.0;

	for I := High(A) downto  Low(A) do

  begin

    DNeg := X * DNeg + Neg;

    Neg := Neg * X;

    DPos := X * DPos + Pos;

    Pos := Pos * X;

    if A[I] >= 0.0 then

      Pos := Pos + A[I]

    else

      Neg := Neg + A[I]

  end;

  Poly.Neg := Neg;

  Poly.Pos := Pos;

  Poly.DNeg := DNeg * X;

  Poly.DPos := DPos * X;

end; {PolyX}





function RelSmall(X, Y: Extended): Boolean;

{ Returns True if X is small relative to Y }

const

  C1: Double = 1E-15;

  C2: Double = 1E-12;

begin

  Result := EAbs(X) < (C1 + C2 * EAbs(Y))

end;



{ Math functions. }



function ArcCos(X: Extended): Extended;

begin

  Result := ArcTan2(Sqrt(1 - X*X), X);

end;



function ArcSin(X: Extended): Extended;

begin

  Result := ArcTan2(X, Sqrt(1 - X*X))

end;



function ArcTan2(Y, X: Extended): Extended;

asm

        FLD     Y

        FLD     X

        FPATAN

        FWAIT

end;



function Tan(X: Extended): Extended;

{  Tan := Sin(X) / Cos(X) }

asm

        FLD    X

        FPTAN

        FSTP   ST(0)      { FPTAN pushes 1.0 after result }

        FWAIT

end;



function CoTan(X: Extended): Extended;

{ CoTan := Cos(X) / Sin(X) = 1 / Tan(X) }

asm

        FLD   X

        FPTAN

        FDIVRP

        FWAIT

end;



function Hypot(X, Y: Extended): Extended;

{ formula: Sqrt(X*X + Y*Y)

  implemented as:  |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision

var

  Temp: Extended;

begin

  X := Abs(X);

  Y := Abs(Y);

  if X > Y then

  begin

    Temp := X;

    X := Y;

    Y := Temp;

  end;

  if X = 0 then

    Result := Y

  else         // Y > X, X <> 0, so Y > 0

    Result := Y * Sqrt(1 + Sqr(X/Y));

end;

}

asm

        FLD     Y

        FABS

        FLD     X

        FABS

        FCOM

        FNSTSW  AX

        TEST    AH,$45

        JNZ      @@1        // if ST > ST(1) then swap

        FXCH    ST(1)      // put larger number in ST(1)

@@1:    FLDZ

        FCOMP

        FNSTSW  AX

        TEST    AH,$40     // if ST = 0, return ST(1)

        JZ      @@2

        FSTP    ST         // eat ST(0)

        JMP     @@3

@@2:    FDIV    ST,ST(1)   // ST := ST / ST(1)

        FMUL    ST,ST      // ST := ST * ST

        FLD1

        FADD               // ST := ST + 1

        FSQRT              // ST := Sqrt(ST)

        FMUL               // ST(1) := ST * ST(1); Pop ST

@@3:    FWAIT

end;





procedure SinCos(Theta: Extended; var Sin, Cos: Extended);

asm

        FLD     Theta

        FSINCOS

        FSTP    tbyte ptr [edx]    // Cos

        FSTP    tbyte ptr [eax]    // Sin

        FWAIT

end;



{ Extract exponent and mantissa from X }

procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);

{ Mantissa ptr in EAX, Exponent ptr in EDX }

asm

        FLD     X

        PUSH    EAX

        MOV     dword ptr [edx], 0    { if X = 0, return 0 }



        FTST

        FSTSW   AX

        FWAIT

        SAHF

        JZ      @@Done



        FXTRACT                 // ST(1) = exponent, (pushed) ST = fraction

        FXCH



// The FXTRACT instruction normalizes the fraction 1 bit higher than

// wanted for the definition of frexp() so we need to tweak the result

// by scaling the fraction down and incrementing the exponent.



        FISTP   dword ptr [edx]

        FLD1

        FCHS

        FXCH

        FSCALE                  // scale fraction

        INC     dword ptr [edx] // exponent biased to match

        FSTP ST(1)              // discard -1, leave fraction as TOS



@@Done:

        POP     EAX

        FSTP    tbyte ptr [eax]

        FWAIT

end;



function Ldexp(X: Extended; P: Integer): Extended;

  { Result := X * (2^P) }

asm

        PUSH    EAX

        FILD    dword ptr [ESP]

        FLD     X

        FSCALE

        POP     EAX

        FSTP    ST(1)

        FWAIT

end;



function Ceil(X: Extended): Integer;

begin

  Result := Integer(Trunc(X));

  if Frac(X) > 0 then

    Inc(Result);

end;



function Floor(X: Extended): Integer;

begin

  Result := Integer(Trunc(X));

  if Frac(X) < 0 then

    Dec(Result);

end;



{ Conversion of bases:  Log.b(X) = Log.a(X) / Log.a(b)  }



function Log10(X: Extended): Extended;

  { Log.10(X) := Log.2(X) * Log.10(2) }

asm

        FLDLG2     { Log base ten of 2 }

        FLD     X

        FYL2X

        FWAIT

end;



function Log2(X: Extended): Extended;

asm

        FLD1

        FLD     X

        FYL2X

        FWAIT

end;



function LogN(Base, X: Extended): Extended;

{ Log.N(X) := Log.2(X) / Log.2(N) }

asm

        FLD1

        FLD     X

        FYL2X

        FLD1

        FLD     Base

        FYL2X

        FDIV

        FWAIT

end;



function Poly(X: Extended; const Coefficients: array of Double): Extended;

{ Horner's method }

var

  I: Integer;

begin

  Result := Coefficients[High(Coefficients)];

  for I := High(Coefficients)-1 downto Low(Coefficients) do

    Result := Result * X + Coefficients[I];

end;



function Power(Base, Exponent: Extended): Extended;

begin

  if Exponent = 0.0 then

    Result := 1.0               { n**0 = 1 }

  else if (Base = 0.0) and (Exponent > 0.0) then

    Result := 0.0               { 0**n = 0, n > 0 }

  else if (Frac(Exponent) = 0.0) and (EAbs(Exponent) <= MaxInt) then

    Result := IntPower(Base, Integer(Trunc(Exponent)))

  else

    Result := Exp(Exponent * Ln(Base))

end;



{ Hyperbolic functions }



function CoshSinh(X: Extended; Factor: Double): Extended;

begin

  Result := Exp(X) / 2;

  Result := Result + Factor / Result;

end;



function Cosh(X: Extended): Extended;

begin

  Result := CoshSinh(X, 0.25)

end;



function Sinh(X: Extended): Extended;

begin

  Result := CoshSinh(X, -0.25)

end;



const

  MaxTanhDomain = 5678.22249441322; // Ln(MaxExtended)/2



function Tanh(X: Extended): Extended;

begin

  if X > MaxTanhDomain then

    Result := 1.0

  else if X < -MaxTanhDomain then

    Result := -1.0

  else

  begin

    Result := Exp(X);

    Result := Result * Result;

    Result := (Result - 1.0) / (Result + 1.0)

  end;

end;



function ArcCosh(X: Extended): Extended;

begin

  if X <= 1.0 then

    Result := 0.0

  else if X > 1.0e10 then

    Result := Ln(2) + Ln(X)

  else

    Result := Ln(X + Sqrt((X - 1.0) * (X + 1.0)));

end;



function ArcSinh(X: Extended): Extended;

var

  Neg: Boolean;

begin

  if X = 0 then

    Result := 0

  else

  begin

    Neg := (X < 0);

    X := EAbs(X);

    if X > 1.0e10 then

      Result := Ln(2) + Ln(X)

    else

    begin

      Result := X*X;

      Result := LnXP1(X + Result / (1 + Sqrt(1 + Result)));

    end;

    if Neg then Result := -Result;

  end;

end;



function ArcTanh(X: Extended): Extended;

var

  Neg: Boolean;

begin

  if X = 0 then

    Result := 0

  else

  begin

    Neg := (X < 0);

    X := EAbs(X);

    if X >= 1 then

      Result := MaxExtended

    else

      Result := 0.5 * LnXP1((2.0 * X) / (1.0 - X));

    if Neg then Result := -Result;

  end;

end;



{ Statistical functions }



function Mean(const Data: array of Double): Extended;

begin

  Result := SUM(Data) / (High(Data) - Low(Data) + 1)

end;



function MinValue(const Data: array of Double): Double;

var

  I: Integer;

begin

  Result := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do

    if Result > Data[I] then

      Result := Data[I];

end;



function MinIntValue(const Data: array of Integer): Integer;

var

  I: Integer;

begin

  Result := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do

    if Result > Data[I] then

      Result := Data[I];

end;



function Min(A,B: Integer): Integer;

begin

  if A < B then

    Result := A

  else

    Result := B;

end;



{$IFDEF _D4orHigher}

function Min(A,B: I64): I64;

begin

  if Cmp64( A, B ) < 0 then

    Result := A

  else

    Result := B;

end;



function Min(A,B: Single): Single;

begin

  if A < B then

    Result := A

  else

    Result := B;

end;



function Min(A,B: Double): Double;

begin

  if A < B then

    Result := A

  else

    Result := B;

end;



function Min(A,B: Extended): Extended;

begin

  if A < B then

    Result := A

  else

    Result := B;

end;

{$ENDIF}



function MaxValue(const Data: array of Double): Double;

var

  I: Integer;

begin

  Result := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do

    if Result < Data[I] then

      Result := Data[I];

end;



function MaxIntValue(const Data: array of Integer): Integer;

var

  I: Integer;

begin

  Result := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do

    if Result < Data[I] then

      Result := Data[I];

end;



function Max(A,B: Integer): Integer;

begin

  if A > B then

    Result := A

  else

    Result := B;

end;



{$IFDEF _D4orHigher}

function Max(A,B: I64): I64;

begin

  if Cmp64( A, B ) > 0 then

    Result := A

  else

    Result := B;

end;



function Max(A,B: Single): Single;

begin

  if A > B then

    Result := A

  else

    Result := B;

end;



function Max(A,B: Double): Double;

begin

  if A > B then

    Result := A

  else

    Result := B;

end;



function Max(A,B: Extended): Extended;

begin

  if A > B then

    Result := A

  else

    Result := B;

end;

{$ENDIF}



procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);

var

  S: Extended;

  N,I: Integer;

begin

  N := High(Data)- Low(Data) + 1;

  if N = 1 then

  begin

    Mean := Data[0];

    StdDev := Data[0];

    Exit;

  end;

  Mean := Sum(Data) / N;

  S := 0;               // sum differences from the mean, for greater accuracy

  for I := Low(Data) to High(Data) do

    S := S + Sqr(Mean - Data[I]);

  StdDev := Sqrt(S / (N - 1));

end;



procedure MomentSkewKurtosis(const Data: array of Double;

  var M1, M2, M3, M4, Skew, Kurtosis: Extended);

var

  Sum, SumSquares, SumCubes, SumQuads, OverN, Accum, M1Sqr, S2N, S3N: Extended;

  I: Integer;

begin

  OverN := 1 / (High(Data) - Low(Data) + 1);

  Sum := 0;

  SumSquares := 0;

  SumCubes := 0;

  SumQuads := 0;

  for I := Low(Data) to High(Data) do

  begin

    Sum := Sum + Data[I];

    Accum := Sqr(Data[I]);

    SumSquares := SumSquares + Accum;

    Accum := Accum*Data[I];

    SumCubes := SumCubes + Accum;

    SumQuads := SumQuads + Accum*Data[I];

  end;

  M1 := Sum * OverN;

  M1Sqr := Sqr(M1);

  S2N := SumSquares * OverN;

  S3N := SumCubes * OverN;

  M2 := S2N - M1Sqr;

  M3 := S3N - (M1 * 3 * S2N) + 2*M1Sqr*M1;

  M4 := (SumQuads * OverN) - (M1 * 4 * S3N) + (M1Sqr*6*S2N - 3*Sqr(M1Sqr));

  Skew := M3 * Power(M2, -3/2);   // = M3 / Power(M2, 3/2)

  Kurtosis := M4 / Sqr(M2);

end;



function Norm(const Data: array of Double): Extended;

begin

  Result := Sqrt(SumOfSquares(Data));

end;



function PopnStdDev(const Data: array of Double): Extended;

begin

  Result := Sqrt(PopnVariance(Data))

end;



function PopnVariance(const Data: array of Double): Extended;

begin

  Result := TotalVariance(Data) / (High(Data) - Low(Data) + 1)

end;



function RandG(Mean, StdDev: Extended): Extended;

{ Marsaglia-Bray algorithm }

var

  U1, S2: Extended;

begin

  repeat

    U1 := 2*Random - 1;

    S2 := Sqr(U1) + Sqr(2*Random-1);

  until S2 < 1;

  Result := Sqrt(-2*Ln(S2)/S2) * U1 * StdDev + Mean;

end;



function StdDev(const Data: array of Double): Extended;

begin

  Result := Sqrt(Variance(Data))

end;



procedure RaiseOverflowError; forward;



function SumInt(const Data: array of Integer): Integer;

{var

  I: Integer;

begin

  Result := 0;

  for I := Low(Data) to High(Data) do

    Result := Result + Data[I]

end; }

asm  // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1

     // loop unrolled 4 times, 5 clocks per loop, 1.2 clocks per datum

      PUSH EBX

      MOV  ECX, EAX         // ecx = ptr to data

      MOV  EBX, EDX

      XOR  EAX, EAX

      AND  EDX, not 3

      AND  EBX, 3

      SHL  EDX, 2

      JMP  @Vector.Pointer[EBX*4]

@Vector:

      DD @@1

      DD @@2

      DD @@3

      DD @@4

@@4:

      ADD  EAX, [ECX+12+EDX]

      JO   RaiseOverflowError

@@3:

      ADD  EAX, [ECX+8+EDX]

      JO   RaiseOverflowError

@@2:

      ADD  EAX, [ECX+4+EDX]

      JO   RaiseOverflowError

@@1:

      ADD  EAX, [ECX+EDX]

      JO   RaiseOverflowError

      SUB  EDX,16

      JNS  @@4

      POP  EBX

end;



procedure RaiseOverflowError;

begin

  raise Exception.Create(e_IntOverflow, SIntOverflow);

end;



function SUM(const Data: array of Double): Extended;

{var

  I: Integer;

begin

  Result := 0.0;

  for I := Low(Data) to High(Data) do

    Result := Result + Data[I]

end; }

asm  // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1

     // Uses 4 accumulators to minimize read-after-write delays and loop overhead

     // 5 clocks per loop, 4 items per loop = 1.2 clocks per item

       FLDZ

       MOV      ECX, EDX

       FLD      ST(0)

       AND      EDX, not 3

       FLD      ST(0)

       AND      ECX, 3

       FLD      ST(0)

       SHL      EDX, 3      // count * sizeof(Double) = count * 8

       JMP      @Vector.Pointer[ECX*4]

@Vector:

       DD @@1

       DD @@2

       DD @@3

       DD @@4

@@4:   FADD     qword ptr [EAX+EDX+24]    // 1

       FXCH     ST(3)                     // 0

@@3:   FADD     qword ptr [EAX+EDX+16]    // 1

       FXCH     ST(2)                     // 0

@@2:   FADD     qword ptr [EAX+EDX+8]     // 1

       FXCH     ST(1)                     // 0

@@1:   FADD     qword ptr [EAX+EDX]       // 1

       FXCH     ST(2)                     // 0

       SUB      EDX, 32

       JNS      @@4

       FADDP    ST(3),ST                  // ST(3) := ST + ST(3); Pop ST

       FADD                               // ST(1) := ST + ST(1); Pop ST

       FADD                               // ST(1) := ST + ST(1); Pop ST

       FWAIT

end;



function SumOfSquares(const Data: array of Double): Extended;

var

  I: Integer;

begin

  Result := 0.0;

  for I := Low(Data) to High(Data) do

    Result := Result + Sqr(Data[I]);

end;



procedure SumsAndSquares(const Data: array of Double; var Sum, SumOfSquares: Extended);

{var

  I: Integer;

begin

  Sum := 0;

  SumOfSquares := 0;

  for I := Low(Data) to High(Data) do

  begin

    Sum := Sum + Data[I];

    SumOfSquares := SumOfSquares + Data[I]*Data[I];

  end;

end;  }

asm  // IN:  EAX = ptr to Data

     //      EDX = High(Data) = Count - 1

     //      ECX = ptr to Sum

     // Est. 17 clocks per loop, 4 items per loop = 4.5 clocks per data item

       FLDZ                 // init Sum accumulator

       PUSH     ECX

       MOV      ECX, EDX

       FLD      ST(0)       // init Sqr1 accum.

       AND      EDX, not 3

       FLD      ST(0)       // init Sqr2 accum.

       AND      ECX, 3

       FLD      ST(0)       // init/simulate last data item left in ST

       SHL      EDX, 3      // count * sizeof(Double) = count * 8

       JMP      @Vector.Pointer[ECX*4]

@Vector:

       DD @@1

       DD @@2

       DD @@3

       DD @@4

@@4:   FADD                            // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4

       FLD     qword ptr [EAX+EDX+24]  // Load Data1

       FADD    ST(3),ST                // Sum := Sum + Data1

       FMUL    ST,ST                   // Data1 := Sqr(Data1)

@@3:   FLD     qword ptr [EAX+EDX+16]  // Load Data2

       FADD    ST(4),ST                // Sum := Sum + Data2

       FMUL    ST,ST                   // Data2 := Sqr(Data2)

       FXCH                            // Move Sqr(Data1) into ST(0)

       FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data1); Pop Data1

@@2:   FLD     qword ptr [EAX+EDX+8]   // Load Data3

       FADD    ST(4),ST                // Sum := Sum + Data3

       FMUL    ST,ST                   // Data3 := Sqr(Data3)

       FXCH                            // Move Sqr(Data2) into ST(0)

       FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data2); Pop Data2

@@1:   FLD     qword ptr [EAX+EDX]     // Load Data4

       FADD    ST(4),ST                // Sum := Sum + Data4

       FMUL    ST,ST                   // Sqr(Data4)

       FXCH                            // Move Sqr(Data3) into ST(0)

       FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data3); Pop Data3

       SUB     EDX,32

       JNS     @@4

       FADD                         // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4

       POP     ECX

       FADD                         // Sqr1 := Sqr2 + Sqr1; Pop Sqr2

       FXCH                         // Move Sum1 into ST(0)

       MOV     EAX, SumOfSquares

       FSTP    tbyte ptr [ECX]      // Sum := Sum1; Pop Sum1

       FSTP    tbyte ptr [EAX]      // SumOfSquares := Sum1; Pop Sum1

       FWAIT

end;



function TotalVariance(const Data: array of Double): Extended;

var

  Sum, SumSquares: Extended;

begin

  SumsAndSquares(Data, Sum, SumSquares);

  Result := SumSquares - Sqr(Sum)/(High(Data) - Low(Data) + 1);

end;



function Variance(const Data: array of Double): Extended;

begin

  Result := TotalVariance(Data) / (High(Data) - Low(Data))

end;





{ Depreciation functions. }



function DoubleDecliningBalance(Cost, Salvage: Extended; Life, Period: Integer): Extended;

{ dv := cost * (1 - 2/life)**(period - 1)

  DDB = (2/life) * dv

  if DDB > dv - salvage then DDB := dv - salvage

  if DDB < 0 then DDB := 0

}

var

  DepreciatedVal, Factor: Extended;

begin

  Result := 0;

	if (Period < 1) or (Life < Period) or (Life < 1) or (Cost <= Salvage) then

    Exit;



	{depreciate everything in period 1 if life is only one or two periods}

	if ( Life <= 2 ) then

  begin

		if ( Period = 1 ) then

      DoubleDecliningBalance:=Cost-Salvage

		else

			DoubleDecliningBalance:=0; {all depreciation occurred in first period}

		exit;

  end;

  Factor := 2.0 / Life;



  DepreciatedVal := Cost * IntPower((1.0 - Factor), Period - 1);

	   {DepreciatedVal is Cost-(sum of previous depreciation results)}



  Result := Factor * DepreciatedVal;

	   {Nominal computed depreciation for this period.  The rest of the

			function applies limits to this nominal value. }



	{Only depreciate until total depreciation equals cost-salvage.}

  if Result > DepreciatedVal - Salvage then

    Result := DepreciatedVal - Salvage;



	{No more depreciation after salvage value is reached.  This is mostly a nit.

	 If Result is negative at this point, it's very close to zero.}

  if Result < 0.0 then

    Result := 0.0;

end;



function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;

{ Spreads depreciation linearly over life. }

begin

  if Life < 1 then ArgError('SLNDepreciation');

  Result := (Cost - Salvage) / Life

end;



function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;

{ SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }

{ Note: life*(life+1)/2 = 1+2+3+...+life "sum of years"

				The depreciation factor varies from life/sum_of_years in first period = 1

																			 downto  1/sum_of_years in last period = life.

				Total depreciation over life is cost-salvage.}

var

  X1, X2: Extended;

begin

  Result := 0;

	if (Period < 1) or (Life < Period) or (Cost <= Salvage) then Exit;

  X1 := 2 * (Life - Period + 1);

  X2 := Life * (Life + 1);

  Result := (Cost - Salvage) * X1 / X2

end;



{ Discounted cash flow functions. }



function InternalRateOfReturn(Guess: Extended; const CashFlows: array of Double): Extended;

{

Use Newton's method to solve NPV = 0, where NPV is a polynomial in

x = 1/(1+rate).  Split the coefficients into negative and postive sets:

  neg + pos = 0, so pos = -neg, so  -neg/pos = 1

Then solve:

  log(-neg/pos) = 0



  Let  t = log(1/(1+r) = -LnXP1(r)

  then r = exp(-t) - 1

Iterate on t, then use the last equation to compute r.

}

var

  T, Y: Extended;

  Poly: TPoly;

  K, Count: Integer;



  function ConditionP(const CashFlows: array of Double): Integer;

  { Guarantees existence and uniqueness of root.  The sign of payments

    must change exactly once, the net payout must be always > 0 for

    first portion, then each payment must be >= 0.

    Returns: 0 if condition not satisfied, > 0 if condition satisfied

    and this is the index of the first value considered a payback. }

  var

    X: Double;

    I, K: Integer;

  begin

    K := High(CashFlows);

    while (K >= 0) and (CashFlows[K] >= 0.0) do Dec(K);

    Inc(K);

    if K > 0 then

    begin

      X := 0.0;

      I := 0;

      while I < K do begin

	      X := X + CashFlows[I];

        if X >= 0.0 then

        begin

     	    K := 0;

      	  Break

    	  end;

      	Inc(I)

      end

    end;

    ConditionP := K

  end;



begin

  InternalRateOfReturn := 0;

  K := ConditionP(CashFlows);

  if K < 0 then ArgError('InternalRateOfReturn');

  if K = 0 then

  begin

    if Guess <= -1.0 then ArgError('InternalRateOfReturn');

    T := -LnXP1(Guess)

  end else

    T := 0.0;

  for Count := 1 to MaxIterations do

  begin

    PolyX(CashFlows, Exp(T), Poly);

    if Poly.Pos <= Poly.Neg then ArgError('InternalRateOfReturn');

    if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then

    begin

      InternalRateOfReturn := -1.0;

      Exit;

    end;

    with Poly do

      Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);

    T := T - Y;

    if RelSmall(Y, T) then

    begin

      InternalRateOfReturn := Exp(-T) - 1.0;

      Exit;

    end

  end;

  ArgError('InternalRateOfReturn');

end;



function NetPresentValue(Rate: Extended; const CashFlows: array of Double;

  PaymentTime: TPaymentTime): Extended;

{ Caution: The sign of NPV is reversed from what would be expected for standard

   cash flows!}

var

  rr: Extended;

  I: Integer;

begin

  if Rate <= -1.0 then ArgError('NetPresentValue');

  rr := 1/(1+Rate);

  result := 0;

  for I := High(CashFlows) downto Low(CashFlows) do

    result := rr * result + CashFlows[I];

  if PaymentTime = ptEndOfPeriod then result := rr * result;

end;



{ Annuity functions. }



{---------------

From the point of view of A, amounts received by A are positive and

amounts disbursed by A are negative (e.g. a borrower's loan repayments

are regarded by the borrower as negative).



Given interest rate r, number of periods n:

  compound(r, n) = (1 + r)**n               "Compounding growth factor"

  annuity(r, n) = (compound(r, n)-1) / r   "Annuity growth factor"



Given future value fv, periodic payment pmt, present value pv and type

of payment (start, 1 , or end of period, 0) pmtTime, financial variables satisfy:



  fv = -pmt*(1 + r*pmtTime)*annuity(r, n) - pv*compound(r, n)



For fv, pv, pmt:



  C := compound(r, n)

  A := (1 + r*pmtTime)*annuity(r, n)

  Compute both at once in Annuity2.



  if C > 1E16 then A = C/r, so:

    fv := meaningless

    pv := -pmt*(pmtTime+1/r)

    pmt := -pv*r/(1 + r*pmtTime)

  else

    fv := -pmt(1+r*pmtTime)*A - pv*C

    pv := (-pmt(1+r*pmtTime)*A - fv)/C

    pmt := (-pv*C-fv)/((1+r*pmtTime)*A)

---------------}



function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,

  FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):

  Extended;

var

		Crn:extended; { =Compound(Rate,NPeriods) }

		Crp:extended; { =Compound(Rate,Period-1) }

		Arn:extended; { =AnnuityF(Rate,NPeriods) }



begin

  if Rate <= -1.0 then ArgError('PaymentParts');

  Crp:=Compound(Rate,Period-1);

  Arn:=Annuity2(Rate,NPeriods,PaymentTime,Crn);

  IntPmt:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;

  PaymentParts:=(-FutureValue-PresentValue)*Crp/Arn;

end;



function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:

  Extended; PaymentTime: TPaymentTime): Extended;

var

  Annuity, CompoundRN: Extended;

begin

  if Rate <= -1.0 then ArgError('FutureValue');

  Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);

  if CompoundRN > 1.0E16 then ArgError('FutureValue');

  FutureValue := -Payment * Annuity - PresentValue * CompoundRN

end;



function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,

  FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

var

  Crp:extended; { compound(rate,period-1)}

  Crn:extended; { compound(rate,nperiods)}

  Arn:extended; { annuityf(rate,nperiods)}

begin

  if (Rate <= -1.0)

    or (Period < 1) or (Period > NPeriods) then ArgError('InterestPayment');

	Crp:=Compound(Rate,Period-1);

	Arn:=Annuity2(Rate,Nperiods,PaymentTime,Crn);

	InterestPayment:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;

end;



function InterestRate(NPeriods: Integer;

  Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{

Given:

  First and last payments are non-zero and of opposite signs.

  Number of periods N >= 2.

Convert data into cash flow of first, N-1 payments, last with

first < 0, payment > 0, last > 0.

Compute the IRR of this cash flow:

  0 = first + pmt*x + pmt*x**2 + ... + pmt*x**(N-1) + last*x**N

where x = 1/(1 + rate).

Substitute x = exp(t) and apply Newton's method to

  f(t) = log(pmt*x + ... + last*x**N) / -first

which has a unique root given the above hypotheses.

}

var

  X, Y, Z, First, Pmt, Last, T, ET, EnT, ET1: Extended;

  Count: Integer;

  Reverse: Boolean;



  function LostPrecision(X: Extended): Boolean;

  asm

        XOR     EAX, EAX

        MOV     BX,WORD PTR X+8

        INC     EAX

        AND     EBX, $7FF0

        JZ      @@1

        CMP     EBX, $7FF0

        JE      @@1

        XOR     EAX,EAX

  @@1:

  end;



begin

  Result := 0;

  if NPeriods <= 0 then ArgError('InterestRate');

  Pmt := Payment;

  if PaymentTime = ptEndOfPeriod then

  begin

    X := PresentValue;

    Y := FutureValue + Payment

  end

  else

  begin

    X := PresentValue + Payment;

    Y := FutureValue

  end;

  First := X;

  Last := Y;

  Reverse := False;

  if First * Payment > 0.0 then

  begin

    Reverse := True;

    T := First;

    First := Last;

    Last := T

  end;

  if first > 0.0 then

  begin

    First := -First;

    Pmt := -Pmt;

    Last := -Last

  end;

  if (First = 0.0) or (Last < 0.0) then ArgError('InterestRate');

  T := 0.0;                     { Guess at solution }

  for Count := 1 to MaxIterations do

  begin

    EnT := Exp(NPeriods * T);

    if {LostPrecision(EnT)} ent=(ent+1) then

    begin

      Result := -Pmt / First;

      if Reverse then

        Result := Exp(-LnXP1(Result)) - 1.0;

      Exit;

    end;

    ET := Exp(T);

    ET1 := ET - 1.0;

    if ET1 = 0.0 then

    begin

      X := NPeriods;

      Y := X * (X - 1.0) / 2.0

    end

    else

    begin

      X := ET * (Exp((NPeriods - 1) * T)-1.0) / ET1;

      Y := (NPeriods * EnT - ET - X * ET) / ET1

    end;

    Z := Pmt * X + Last * EnT;

    Y := Ln(Z / -First) / ((Pmt * Y + Last * NPeriods *EnT) / Z);

    T := T - Y;

    if RelSmall(Y, T) then

    begin

      if not Reverse then T := -T;

      InterestRate := Exp(T)-1.0;

      Exit;

    end

  end;

  ArgError('InterestRate');

end;



function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;

  PaymentTime: TPaymentTime): Extended;



{ If Rate = 0 then nper := -(pv + fv) / pmt

  else cf := pv + pmt * (1 + rate*pmtTime) / rate

       nper := LnXP1(-(pv + fv) / cf) / LnXP1(rate) }



var

  PVRPP: Extended; { =PV*Rate+Payment } {"initial cash flow"}

  T:     Extended;



begin



  if Rate <= -1.0 then ArgError('NumberOfPeriods');



{whenever both Payment and PaymentTime are given together, the PaymentTime has the effect

 of modifying the effective Payment by the interest accrued on the Payment}



  if ( PaymentTime=ptStartOfPeriod ) then

    Payment:=Payment*(1+Rate);



{if the payment exactly matches the interest accrued periodically on the

 presentvalue, then an infinite number of payments are going to be

 required to effect a change from presentvalue to futurevalue.  The

 following catches that specific error where payment is exactly equal,

 but opposite in sign to the interest on the present value.  If PVRPP

 ("initial cash flow") is simply close to zero, the computation will

 be numerically unstable, but not as likely to cause an error.}



  PVRPP:=PresentValue*Rate+Payment;

  if PVRPP=0 then ArgError('NumberOfPeriods');



  { 6.1E-5 approx= 2**-14 }

  if ( EAbs(Rate)<6.1E-5 ) then

    Result:=-(PresentValue+FutureValue)/PVRPP

  else

  begin



{starting with the initial cash flow, each compounding period cash flow

 should result in the current value approaching the final value.  The

 following test combines a number of simultaneous conditions to ensure

 reasonableness of the cashflow before computing the NPER.}



    T:= -(PresentValue+FutureValue)*Rate/PVRPP;

    if  T<=-1.0  then ArgError('NumberOfPeriods');

    Result := LnXP1(T) / LnXP1(Rate)

  end;

  NumberOfPeriods:=Result;

end;



function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:

  Extended; PaymentTime: TPaymentTime): Extended;

var

  Annuity, CompoundRN: Extended;

begin

  if Rate <= -1.0 then ArgError('Payment');

  Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);

  if CompoundRN > 1.0E16 then

    Payment := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)

  else

    Payment := (-PresentValue * CompoundRN - FutureValue) / Annuity

end;



function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;

  PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

var

  Junk: Extended;

begin

  if (Rate <= -1.0) or (Period < 1) or (Period > NPeriods) then ArgError('PeriodPayment');

  PeriodPayment := PaymentParts(Period, NPeriods, Rate, PresentValue,

       FutureValue, PaymentTime, Junk);

end;



function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:

  Extended; PaymentTime: TPaymentTime): Extended;

var

  Annuity, CompoundRN: Extended;

begin

  if Rate <= -1.0 then ArgError('PresentValue');

  Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);

  if CompoundRN > 1.0E16 then

    PresentValue := -(Payment / Rate * Integer(PaymentTime) * Payment)

  else

    PresentValue := (-Payment * Annuity - FutureValue) / CompoundRN

end;



end.