//ÛßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÛ
//Û                                                       Û
//Û      Virtual Pascal Runtime Library.  Version 2.1.    Û
//Û      Mathematical function library                    Û
//Û      ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÛ
//Û      Copyright (C) 1996 NITEK Corporation             Û
//Û      Included in Virtual Pascal with permission       Û
//Û                                                       Û
//ßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßß

unit Math;

{ This unit contains high-performance arithmetic, trigonometric, logorithmic,
  statistical and financial calculation routines which supplement the math
  routines that are part of the Virtual Pascal/2 language or System unit. }

{&X+}

interface

uses SysUtils;

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;

  {-----------------------------------------------------------------------
  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.
  -----------------------------------------------------------------------}

  { 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);
  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;

  { 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);

  { Ldexp: returns X*2**P }
  function Ldexp(X: Extended; P: Integer): Extended;

  { Ceil: Smallest integer >= X, |X| < MaxInt }
  function Ceil(X: Extended):LongInt;

  { Floor: Largest integer <= X,  |X| < MaxInt }
  function Floor(X: Extended): LongInt;

  { 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;
  function SumOfSquares(const Data: array of Double): Extended;
  procedure SumsAndSquares(const Data: array of Double;
    var Sum, SumOfSquares: Extended);

  { MinValue: Returns the smallest signed value in the data array (MIN) }
  function MinValue(const Data: array of Double): Double;

  { MaxValue: Returns the largest signed value in the data array (MAX) }
  function MaxValue(const Data: array of Double): Double;

  { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  function StdDev(const Data: array of Double): Extended;

  { MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
    calculating them separately.  Less accurate when the mean is very large
    (> 10e7) or the variance is very small. }
  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

  type
    TPoly = record
      Neg, Pos, DNeg, DPos: Extended
    end;

  const
    MaxIterations = 15;

  procedure ArgError;
  begin
    raise EInvalidArgument.Create('');
  end;

  function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }
  {&Frame-} {&Uses none}
  asm
    fld     X         // c
    fld     st        // c c
    fmul    st,st     // c*c c
    fld1              // 1 c*c c
    fsubrp  st(1),st  // (1-c*c) c
    fsqrt             // s c
    fxch    st(1)     // c s
    fpatan
  end;

  function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }
  {&Frame-} {&Uses none}
  asm
    fld     X         // s
    fld     st        // s s
    fmul    st,st     // s*s s
    fld1              // 1 s*s s
    fsubrp  st(1),st  // (1-s*s) s
    fsqrt             // c s
    fpatan
  end;

  { 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;
  {&Frame-} {&Uses none}
  asm
    fld     Y
    fld     X
    fpatan
  end;

  { SinCos is 2x faster than calling Sin and Cos separately for the same angle }
  procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
  {&Frame-} {&Uses none}
  asm
    fld     Theta
    fsincos
    mov     eax, Cos
    fstp    tbyte ptr [eax]
    mov     eax, Sin
    fstp    tbyte ptr [eax]
    fwait
  end;

  function Tan(X: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld     X
    fsincos
    fdivp   st(1),st
  end;

  function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }
  {&Frame-} {&Uses none}
  asm
    fld     X
    fsincos
    fdivrp  st(1),st
  end;

  function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }
  {&Frame-} {&Uses none}
  asm
    fld     X
    fmul    st,st
    fld     Y
    fmul    st,st
    faddp   st(1),st
    fsqrt
  end;

  { Angle unit conversion routines }
  function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}
  begin
    DegToRad := Degrees * PI / 180;
  end;

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

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

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

  function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  const
    TwoPI = 2*PI;
  begin
    CycleToRad := Cycles * TwoPI
  end;

  function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  const
    TwoPI = 2*PI;
  begin
    RadToCycle := Radians / TwoPI;
  end;


  { Hyperbolic functions and inverses }

  function Cosh(X: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld      X
    fldl2e
    fmulp    st(1),st
    fabs
    fld      st
    frndint
    fsub     st(1),st
    fxch     st(1)
    ftst
    fstsw    ax
    sahf
    jb       @L0
    f2xm1
    jmp      @L1
  @L0:
    fchs
    f2xm1
    fld1
    fadd     st,st(1)
    fdivp    st(1),st
    fchs
  @L1:
    fld1
    faddp    st(1),st
    fscale
    fstp     st(1)
    fld1
    fdiv     st,st(1)
    faddp    st(1),st
    fld1
    fchs
    fxch     st(1)
    fscale
    fstp     st(1)
  end;

  function Sinh(X: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld      X
    ftst
    fstsw    ax
    mov      dx,ax
    fldl2e
    fmulp    st(1),st
    fabs
    fld      st
    frndint
    fsub     st(1),st
    fchs
    fld1
    fscale
    fxch     st(1)
    fchs
    fxch     st(2)
    ftst
    fstsw    ax
    sahf
    jb       @L0
    f2xm1
    jmp      @L1
  @L0:
    fchs
    f2xm1
    fld1
    fadd     st,st(1)
    fdivp    st(1),st
    fchs
  @L1:
    fld1
    fadd     st,st(1)
    fxch     st(2)
    fld1
    fsubp    st(1),st
    fsubp    st(1),st
    fdivr    st(1),st
    fxch     st(1)
    fxch     st(2)
    fxch     st(1)
    fscale
    fstp     st(1)
    faddp    st(1),st
    fld1
    fchs
    fxch     st(1)
    fscale
    fstp     st(1)
    mov      ax,dx
    sahf
    jae      @L2
    fchs
  @L2:
  end;

  function Tanh(X: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld      X
    ftst
    fstsw    ax
    mov      dx,ax
    fldl2e
    fmulp    st(1),st
    fadd     st,st
    fabs
    fld      st
    frndint
    fsub     st(1),st
    fchs
    fld1
    fscale
    fstp     st(1)
    fld1
    fsubp    st(1),st
    fxch     st(1)
    ftst
    fstsw    ax
    sahf
    jb       @L0
    f2xm1
    jmp      @L1
  @L0:
    fchs
    f2xm1
    fld1
    fadd     st,st(1)
    fdivp    st(1),st
    fchs
  @L1:
    fld      st
    fsub     st,st(2)
    fxch     st(1)
    faddp    st(2),st
    fld1
    fadd     st(2),st
    faddp    st(2),st
    fdivrp   st(1),st
    mov      ax,dx
    sahf
    jae      @L2
    fchs
  @L2:
  end;

  function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }
  {&Frame-} {&Uses none}
  asm
    fld     x
    fld     st        // x x
    fld1              // 1 x x
    fsub    st(2),st  // 1 x (x-1)
    faddp   st(1),st  // (x+1) (x-1)
    fld     st(1)     // (x-1) (x+1) (x-1)
    fmulp   st(1),st  // (x-1)*(x+1) (x-1)
    fsqrt             // sqrt((x+1)*(x-1)) (x-1)
    faddp   st(1),st  // z, now ArcCosh(X) = ln(1+z)
    fldln2            // loge(2) z
    fxch    st(1)     // z loge(2)
    fyl2xp1           // fyl2xp1(fldln2,z)
  end;

  function ArcSinh(X: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld     x
    ftst
    fstsw   ax
    sahf
    jne     @L0
    fstp    st
    fldz
    jmp     @L2
  @L0:
    fabs               // ax
    fld1               // 1 ax
    fdiv    st,st(1)   // 1/ax ax
    fld     st         // 1/ax 1/ax ax
    fmul    st,st      // 1/ax^2 1/ax ax
    fld1               // 1 1/aX^2 1/ax ax
    faddp   st(1),st   // 1+1/ax^2 1/ax ax
    fsqrt              // sqrt(1+1/ax^2) 1/ax ax
    faddp   st(1),st   // 1/ax+sqrt(1+1/ax^2) ax
    fdivr   st,st(1)
    faddp   st(1),st   // z
    fldln2
    jae     @L1
    fchs
  @L1:
    fxch    st(1)
    fyl2xp1
  @L2:
  end;

  function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }
  {&Frame-} {&Uses none}
  asm
    fld     x
    ftst
    fstsw   ax
    sahf
    fabs              // ax
    fld1              // 1 ax
    fsub    st,st(1)  // 1-ax ax
    fdivp   st(1),st  // ax/(1-ax)
    fadd    st,st     // z
    fldln2
    jae     @L0
    fchs
  @L0:
    fxch    st(1)
    fyl2xp1
  end;


  { Logorithmic functions }

  function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }
  {&Frame-} {&Uses none}
  asm
    fld1
    fld     X
    fyl2xp1
    fldln2
    fmulp   st(1),st
  end;

  function Log10(X: Extended): Extended;                     { Log base 10 of X}
  {&Frame-} {&Uses none}
  asm
    fldlg2
    fld      X
    fyl2x
  end;

  function Log2(X: Extended): Extended;                      { Log base 2 of X }
  {&Frame-} {&Uses none}
  asm
    fld1
    fld     X
    fyl2x
  end;

  function LogN(Base, X: Extended): Extended;                { Log base N of X }
  {&Frame-} {&Uses none}
  asm
    fld1
    fld     X
    fyl2x
    fld1
    fld     Base
    fyl2x
    fdivp   st(1),st
  end;


  { Exponential functions }

  { IntPower: Raise base to an integral power.  Fast. }
  function IntPower(Base: Extended; Exponent: Integer): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Exponent
    cdq
    xor     eax,edx
    sub     eax,edx          // eax := Abs(Exponent)
    mov     ecx,eax
    fld     Base
    fld1
    jecxz   @L1
  @L0:
    fmul    st,st(1)
    loop    @L0
  @L1:
    cmp     dword ptr Exponent,0
    jge     @L2
    fld1
    fdivrp                   // Result := 1 / Result
  @L2:
    fstp    st(1)
  end;

  { Power: Raise base to any power.
    For fractional exponents, or exponents > MaxInt, base must be > 0. }
  function Power(Base, Exponent: Extended): Extended;
  {&Frame-} {&Uses none}
  asm
    fld      Exponent
    fld      Base
    fabs
    fyl2x
    fld      st
    frndint
    fsub     st(1),st
    fxch     st(1)
    ftst
    fstsw    ax
    sahf
    jb       @L0
    f2xm1
    jmp      @L1
  @L0:
    fchs
    f2xm1
    fld1
    fadd     st,st(1)
    fdivp    st(1),st
    fchs
  @L1:
    fld1
    faddp    st(1),st
    fscale
    fstp     st(1)
  end;


  { Miscellaneous Routines }

  { Frexp:  Separates the mantissa and exponent of X. }
  procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
  {&Frame-} {&Uses none}
  asm
    fld     X
    fxtract
    mov     eax,Mantissa
    fstp    tbyte ptr [eax]
    mov     eax,Exponent
    fistp   word ptr [eax]
    fwait
  end;

  { Ldexp: returns X*2**P }
  function Ldexp(X: Extended; P: Integer): Extended;
  {&Frame-} {&Uses none}
  asm
    fild    P
    fld     X
    fscale
    fstp    st(1)
  end;

  { Ceil: Smallest integer >= X, |X| < MaxInt }
  function Ceil(X: Extended):LongInt;
  {&Frame-} {&Uses none}
  asm
    fld     x
    fstcw   word ptr x      { save rounding control }
    fwait
    mov     ax, word ptr x
    and     ax, 0f3ffh
    or      ax, 00800h
    xchg    ax, word ptr x
    fldcw   word ptr x      { set rounding toward positive infinity }
    frndint                 { round x }
    fistp   dword ptr x+4   { store as LongInt }
    fwait
    xchg    ax, word ptr x
    fldcw   word ptr x      { reset rounding control }
    mov     eax,dword ptr x+4
  end;

  { Floor: Largest integer <= X,  |X| < MaxInt }
  function Floor(X: Extended): LongInt;
  {&Frame-} {&Uses none}
  asm
    fld     x
    fstcw   word ptr x      { save rounding control }
    fwait
    mov     ax, word ptr x
    and     ax, 0f3ffh
    or      ax, 00400h
    xchg    ax, word ptr x
    fldcw   word ptr x      { set rounding toward negative infinity }
    frndint                 { round x }
    fistp   dword ptr x+4   { store as LongInt }
    fwait
    xchg    ax, word ptr x
    fldcw   word ptr x      { reset rounding control }
    mov     eax,dword ptr x+4
  end;

  { 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;
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Coefficients-4  // ecx = High(Coefficients)
    mov     eax,ecx
    shl     eax,3
    fld     X
    add     eax,Coefficients
    fld     qword ptr [eax]
    jecxz   @L1
  @L0:
    sub     eax,8
    fmul    st,st(1)
    fadd    qword ptr [eax]
    loop    @L0
  @L1:
    fstp    st(1)
  end;

  {-----------------------------------------------------------------------
  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;
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Data-4          // ecx = High(Data)
    inc     ecx
    mov     Data-4,ecx
    mov     eax,Data
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4
  end;

  { Sum: Sum of values.  (SUM) }
  function Sum(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Data-4
    inc     ecx
    mov     eax,Data
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
  end;

  function SumOfSquares(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Data-4
    inc     ecx
    mov     eax,Data
    fldz
  @L0:
    fld     qword ptr [eax]
    fmul    st,st
    add     eax,8
    faddp   st(1),st
    loop    @L0
  end;

  procedure SumsAndSquares(const Data: array of Double;
                           var Sum, SumOfSquares: Extended);
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Data-4
    inc     ecx
    mov     eax,Data
    fldz
    fldz
  @L0:
    fld     qword ptr [eax]
    fadd    st(2),st
    fmul    st,st
    add     eax,8
    faddp   st(1),st
    loop    @L0
    mov     eax,SumOfSquares
    fstp    tbyte ptr [eax]
    mov     eax,Sum
    fstp    tbyte ptr [eax]
    fwait
  end;

  { MinValue: Returns the smallest signed value in the data array (MIN) }
  function MinValue(const Data: array of Double): Double;
  {&Frame-} {&Uses ebx}
  asm
    mov     ebx,Data
    mov     ecx,Data-4
    fld     qword ptr [ebx]
    jecxz   @L2
  @L0:
    add     ebx,8
    fcom    qword ptr [ebx]
    fstsw   ax
    sahf
    jbe     @L1
    fstp    st
    fld     qword ptr [ebx]
  @L1:
    loop    @L0
  @L2:
  end;

  { MaxValue: Returns the largest signed value in the data array (MAX) }
  function MaxValue(const Data: array of Double): Double;
  {&Frame-} {&Uses ebx}
  asm
    push    ebx
    mov     ebx,Data
    mov     ecx,Data-4
    fld     qword ptr [ebx]
    jecxz   @L2
  @L0:
    add     ebx,8
    fcom    qword ptr [ebx]
    fstsw   ax
    sahf
    jae     @L1
    fstp    st
    fld     qword ptr [ebx]
  @L1:
    loop    @L0
  @L2:
    pop     ebx
  end;

  { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  function StdDev(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4 // xb = sx/N
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
    dec     dword ptr Data-4
    fidiv   dword ptr Data-4
    fsqrt
  end;

  { MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
    calculating them separately.  Less accurate when the mean is very large
    (> 10e7) or the variance is very small. }
  procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4
    mov     eax,Mean
    fld     st
    fstp    tbyte ptr [eax]
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
    dec     dword ptr Data-4
    fidiv   dword ptr Data-4
    fsqrt
    mov     eax,StdDev
    fstp    tbyte ptr [eax]
    fwait
  end;

  { Population Standard Deviation (STDP): Sqrt(PopnVariance).
    Used in some business and financial calculations. }
  function PopnStdDev(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4 // xb = sx/N
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
    fidiv   dword ptr Data-4
    fsqrt
  end;

  { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
  function Variance(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4 // xb = sx/N
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
    dec     dword ptr Data-4
    fidiv   dword ptr Data-4
  end;

  { Population Variance (VAR or VARP): TotalVariance/ N }
  function PopnVariance(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4 // xb = sx/N
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
    fidiv   dword ptr Data-4
    fwait
  end;

  { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
  function TotalVariance(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    mov     ecx,Data-4
    inc     ecx
    mov     Data-4,ecx
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4 // xb = sx/N
    mov     eax,Data
    mov     ecx,Data-4
    fldz
  @L1:
    fld     qword ptr [eax]
    fsub    st,st(2)
    fmul    st,st
    faddp   st(1),st
    add     eax,8
    loop @L1
    fstp    st(1)
  end;

  { Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
  function Norm(const Data: array of Double): Extended;
  {&Frame-} {&Uses none}
  asm
    mov     eax,Data
    fldz
    mov     ecx,Data-4
    inc     ecx
  @L0:
    fld     qword ptr [eax]
    fmul    st,st
    add     eax,8
    faddp   st(1),st
    loop    @L0
    fsqrt
  end;

  { 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);
  {&Frame-} {&Uses none}
  asm
    mov     ecx,Data-4
    inc     ecx
    // calculate the mean
    mov     Data-4,ecx
    mov     eax,Data
    fldz
  @L0:
    fadd    qword ptr [eax]
    add     eax,8
    loop    @L0
    fidiv   dword ptr Data-4    // avg
    // now calculate first (absolute), second, third and fourth moments
    // of the deviation from the mean
    mov     ecx,Data-4
    mov     eax,Data
    fldz
    fldz
    fldz
    fldz                        // M1 M2 M3 M4 avg
  @L1:
    fld     qword ptr [eax]     // x M1 M2 M3 M4 avg
    fsub    st,st(5)            // dev M1 M2 M3 M4 avg
    fld     st                  // dev dev M1 M2 M3 M4 avg
    fabs                        // adev dev M1 M2 M3 M4 avg
    fadd    st(2),st            // adev dev M1+ M2 M3 M4 avg
    fmul    st,st               // dev^2 dev M1+ M2 M3 M4 avg
    fmul    st(1),st            // dev^2 dev^3 M1+ M2 M3 M4 avg
    fadd    st(3),st            // dev^2 dev^3 M1+ M2+ M3 M4 avg
    fmul    st,st               // dev^4 dev^3 M1+ M2+ M3 M4 avg
    faddp   st(5),st            // dev^3 M1+ M2+ M3 M4+ avg
    faddp   st(3),st            // M1+ M2+ M3+ M4+ avg
    add     eax,8
    loop    @L1
    mov     eax,M1
    fstp    tbyte ptr [eax]     // M2 M3 M4 avg
    mov     eax,M2
    fstp    tbyte ptr [eax]     // M3 M4 avg
    mov     eax,M3
    fstp    tbyte ptr [eax]     // M4 avg
    mov     eax,M4
    fstp    tbyte ptr [eax]     // avg
    fstp    st
    mov     eax,M2
    fld     tbyte ptr [eax]
    dec     dword ptr Data-4
    fidiv   dword ptr Data-4    // var
    ftst                        // var = 0?
    fstsw   ax
    sahf
    je      @Err
    fld     st                  // var var
    fsqrt                       // sdev var
    mov     eax,M3
    fld     tbyte ptr [eax]     // M3 sdev var
    inc     dword ptr Data-4
    fidiv   dword ptr Data-4    // M3/N sdev var
    fdiv    st,st(1)            // M3/(N*sdev) var
    fdiv    st,st(2)            // Skew=M3/(N*sdev*var) var
    mov     eax,Skew
    fstp    tbyte ptr [eax]     // var
    fmul    st,st               // var*var
    mov     eax,M4
    fld     tbyte ptr [eax]     // M4 var*var
    fidiv   dword ptr Data-4    // M4/N var*var
    fdivrp  st(1),st            // M4/(N*var*var)
    mov     dword ptr Data-4,3
    fisub   dword ptr Data-4    // Kurtosis = M4/(N*var*var) - 3
    mov     eax,Kurtosis
    fstp    tbyte ptr [eax]
    jmp     @Done
  @Err:
    // no skew/kurtosis when var=0
    fldz
    mov     eax,Skew
    fstp    tbyte ptr [eax]
    fldz
    mov     eax,Kurtosis
    fstp    tbyte ptr [eax]
  @Done:
  end;

  { RandG produces random numbers with Gaussian distribution about the mean.
    Useful for simulating data with sampling errors. }
  function RandG(Mean, StdDev: Extended): Extended;
  var
    v1,v2,rsq,fac : Extended;
  begin
    repeat
      v1 := 2.0 * random - 1.0;
      v2 := 2.0 * random - 1.0;
      rsq := sqr(v1) + sqr(v2);
    until (rsq < 1.0) and (rsq <> 0.0);
    fac := sqrt(-2.0*ln(rsq)/rsq);
    RandG := Mean + StdDev*v2*fac;
  end;

  {-----------------------------------------------------------------------
  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

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

  { Double Declining Balance (DDB) }
  function DoubleDecliningBalance(Cost, Salvage: Extended;
    Life, Period: Integer): Extended;
  VAR
    i : Integer;
    c,DDB : Extended;
  begin
    DDB := 0.0;
    for i := 1 TO Period do begin
      c := Cost*(1.0 - 2.0/Life);
      if c < Salvage then c := Salvage;
      DDB := Cost - c;
      Cost := c;
    end;
    DoubleDecliningBalance := DDB;
  end;

  { Future Value (FVAL) }
  function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
    Extended; PaymentTime: TPaymentTime): Extended;
  VAR
    r, rN, rN1, s, p : Extended;
  begin
    r   := 1.0 + Rate;
    rN  := power(r,NPeriods);
    rN1 := r*rN;
    s   := -Payment*(rN1 - r)/(r - 1.0);
    if PaymentTime = ptEndOfPeriod then
      s := s/r;
    p   := -PresentValue*rN;
    FutureValue := s + p;
  end;

  { Interest Payment (IPAYMT)  }
  function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
    FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  var
    rf,r,pf,s,PayOf,PayOf0 : Extended;

    function FofPn(n: Integer) : Extended;
    begin
      FofPn := PresentValue*power(r,n);
    end;

    function FofAn(n: Integer) : Extended;
    begin
      result := pf*(power(r,n)-1.0)/rf;
      if PaymentTime = ptStartOfPeriod then
        result := r * result;
    end;

  begin
    rf  := Rate;
    r   := 1.0+rf;
    // calculate monthly payment pf
    s   := power(r,NPeriods);
    pf  := PresentValue*rf*(s/(s-1.0));
    // calculate outstanding principle
    PayOf := FofPn(Period) - FofAn(Period);
    // calculate previous outstanding principle
    PayOf0 := FofPn(Period-1) - FofAn(Period-1);
    InterestPayment := PayOf0 - PayOf;
  end;

  { Interest Rate (IRATE) }
  function InterestRate(NPeriods: Integer;
    Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

    function e(r: Extended) : Extended;
    var rN, s0, s1, s2 : Extended;
    begin
      rN := power(r,NPeriods);
      if PaymentTime = ptStartOfPeriod then begin
        s0 := rN*(PresentValue + Payment);
        s2 := 0.0;
      end
      else begin
        s0 := rN*PresentValue;
        s2 := Payment;
      end;
      if r = 1.0 then
        s1 := Payment*(NPeriods*rN/r - 1)
      else
        s1 := Payment*(rN-r)/(r-1);
      e := FutureValue + s0 + s1 + s2;
    end;

  var
    r0, r1, dr, e0, e1, de, eok : Extended;
    count : LongInt;
  begin
    if abs(FutureValue) > abs(PresentValue) then
      eok := abs(FutureValue)
    else
      eok := abs(PresentValue);
    r1 := 1.001;
    e1 := e(r1);
    count := 0;
    repeat
      e0 := e1;
      r0 := r1;
      de := ( e(r0+0.000001) - e0 ) / 0.000001;
      dr := -e0/de;
      if dr > 0.05 then dr := 0.05;
      if dr < -0.05 then dr := -0.05;
      e1 := e(r0+dr);
      while abs(e1) > abs(e0) do begin
        dr := 0.5*dr;
        e1 := e(r0+dr);
      end;
      r1 := r0 + dr;
      inc(count);
    until (abs(e1) < 1.0e-14*eok) or (count > 50);
    InterestRate := r1 - 1.0;
  end;

  { Discounted cash flow functions. }

  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) * (1 - Compound(R, N)) / R }
  begin
    if R = 0.0 then
    begin
      CompoundRN := 1.0;
      Result := -N
    end
    else
    begin
      { 6.1E-5 approx= 2**-14 }
      if Abs(R) < 6.1E-5 then
      begin
        CompoundRN := Exp(N * LnXP1(R));
        Result := (CompoundRN - 1.0) / -R;
      end
      else
      begin
        CompoundRN := Compound(R, N);
        Result := (1 - CompoundRN) / R
      end;
  //    if PaymentTime = ptEndOfPeriod then Result := Result + R
      if PaymentTime = ptStartOfPeriod then Result := Result * (1 + R);
    end
  end;

  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 := Low(A) to High(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;

  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 := Abs(X) < (C1 + C2 * Abs(Y))
  end;

  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;
      Result := K
    end;

  begin
    Result := 0;
    K := ConditionP(CashFlows);
    if K < 0 then ArgError;
    if K = 0 then
    begin
      if Guess <= -1.0 then ArgError;
      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;
      if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then
      begin
        Result := -1.0;
        Exit;
      end;
      with Poly do
        Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);
      T := T - Y;
      if RelSmall(Y, T) then
      begin
        Result := Exp(-T) - 1.0;
        Exit;
      end
    end;
    ArgError;
  end;

  { Number of Periods (NPER) }
  function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
    PaymentTime: TPaymentTime): Extended;
  { If Rate = 0 then nper := -(pv + fv) / pmt
    else t := pv + pmt * (1 + rate*pmtTime) / rate
         nper := LnXP1(-(pv + vf) / t) / LnXP1(rate)
  }
  var
    T: Extended;
  begin
    if Frac(Rate) = 0.0 then
      Result := -(PresentValue + FutureValue) / Payment
    else
    begin
      T := PresentValue + Payment * (1 + Integer(PaymentTime) * Rate) / Rate;
      Result := LnXP1(-(PresentValue + FutureValue) / T) / LnXP1(Rate)
    end
  end;

  { Net Present Value. (NPV) Needs array of cash flows. }
  function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
    PaymentTime: TPaymentTime): Extended;
  var
    X, Y: Extended;
    I: Integer;
  begin
    if Rate <= -1.0 then ArgError;
    X := 1.0 / (1.0 + Rate);
    Y := CashFlows[High(CashFlows)];
    for I := High(CashFlows) - 1 downto Low(CashFlows) do
      Y := X * Y + CashFlows[I];
    if PaymentTime = ptStartOfPeriod then Y := X * Y;
    Result := Y;
  end;

  { Payment (PAYMT) }
  function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,
    FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):
    Extended;
  var
    Pmt, Rate1, T1, T2, T3: Extended;
    M: Integer;

    function Annuity(R: Extended; N: Integer): Extended;
    begin
      Result := (1 - Compound(R, N)) / R
    end;

  begin
    M := Period;
    Rate1 := 0.0;
    if PaymentTime = ptEndOfPeriod then
    begin
      Inc(M);
      Rate1 := Rate
    end;
    T1 := (1 + Rate1) * Annuity(Rate, NPeriods);
    T2 := (1 + Rate1) * Annuity(Rate, M);
    T3 := (1 + Rate1) * Annuity(Rate, NPeriods - M);
    Pmt := (FutureValue + PresentValue * Compound(Rate, NPeriods)) / T1;
    IntPmt := Rate *
              (FutureValue * Compound(Rate, M) - PresentValue * T2 * T3) / T1;
    Result := Pmt - IntPmt
  end;

  function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:
    Extended; PaymentTime: TPaymentTime): Extended;
  var
    Annuity, CompoundRN: Extended;
  begin
    Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
    if CompoundRN > 1.0E16 then
      Result := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)
    else
      Result := (PresentValue * CompoundRN + FutureValue) / Annuity
  end;

  { Period Payment (PPAYMT) }
  function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
    PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  var
    Junk: Extended;
  begin
    if (Period < 1) or (Period > NPeriods) then ArgError;
    Result := PaymentParts(Period, NPeriods, Rate, PresentValue, FutureValue,
                           PaymentTime, Junk);
  end;

  { Present Value (PVAL) }
  function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:
    Extended; PaymentTime: TPaymentTime): Extended;
  var
    Annuity, CompoundRN: Extended;
  begin
    Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
    if CompoundRN > 1.0E16 then
      Result := -(Payment / Rate * Integer(PaymentTime) * Payment)
    else
      Result := (Payment * Annuity - FutureValue) / CompoundRN
  end;

  { Straight Line depreciation (SLN) }
  function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  { Spreads depreciation linearly over life. }
  begin
    if Life < 1 then ArgError;
    Result := (Cost - Salvage) / Life
  end;

  { Sum-of-Years-Digits depreciation (SYD) }
  function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  { SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }
  var
    X1, X2: Extended;
  begin
    if (Period < 1) or (Life < Period) then ArgError;
    X1 := 2 * (Life - Period + 1);
    X2 := Life * (Life + 1);
    Result := (Cost - Salvage) * X1 / X2
  end;

end.
