home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / vp21beta.zip / ARTLSRC.RAR / MATH.PAS < prev    next >
Pascal/Delphi Source File  |  2000-09-07  |  44KB  |  1,584 lines
  1. //█▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█
  2. //█                                                       █
  3. //█      Virtual Pascal Runtime Library.  Version 2.1.    █
  4. //█      Mathematical function library                    █
  5. //█      ─────────────────────────────────────────────────█
  6. //█      Copyright (C) 1996 NITEK Corporation             █
  7. //█      Included in Virtual Pascal with permission       █
  8. //█                                                       █
  9. //▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
  10.  
  11. unit Math;
  12.  
  13. { This unit contains high-performance arithmetic, trigonometric, logorithmic,
  14.   statistical and financial calculation routines which supplement the math
  15.   routines that are part of the Virtual Pascal/2 language or System unit. }
  16.  
  17. interface
  18.  
  19. uses SysUtils;
  20.  
  21. const   { Ranges of the IEEE floating point types, including denormals }
  22.   MinSingle        =     1.5e-45;
  23.   MaxSingle        =     3.4e+38;
  24.   MinDouble        =     5.0e-324;
  25.   MaxDouble        =     1.7e+308;
  26.   MinExtended      =     3.4e-4932;
  27.   MaxExtended      =     1.1e+4932;
  28.   MinComp          =     -9.223372036854775807e+18;
  29.   MaxComp          =     9.223372036854775807e+18;
  30.  
  31.   {-----------------------------------------------------------------------
  32.   All angle parameters and results of trig functions are in radians.
  33.  
  34.   Most of the following trig and log routines map directly to Intel 80387 FPU
  35.   floating point machine instructions.  Input domains, output ranges, and
  36.   error handling are determined largely by the FPU hardware.
  37.   Routines coded in assembler favor the Pentium FPU pipeline architecture.
  38.   -----------------------------------------------------------------------}
  39.  
  40.   { Trigonometric functions }
  41.   function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }
  42.   function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }
  43.  
  44.   { ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
  45.     IN: |Y| < 2^64, |X| < 2^64, X <> 0   OUT: [-PI..PI] radians }
  46.   function ArcTan2(Y, X: Extended): Extended;
  47.  
  48.   { SinCos is 2x faster than calling Sin and Cos separately for the same angle }
  49.   procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
  50.   function Tan(X: Extended): Extended;
  51.   function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }
  52.   function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }
  53.  
  54.   { Angle unit conversion routines }
  55.   function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}
  56.   function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
  57.   function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
  58.   function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }
  59.   function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  60.   function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  61.  
  62.   { Hyperbolic functions and inverses }
  63.   function Cosh(X: Extended): Extended;
  64.   function Sinh(X: Extended): Extended;
  65.   function Tanh(X: Extended): Extended;
  66.   function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }
  67.   function ArcSinh(X: Extended): Extended;
  68.   function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }
  69.  
  70.   { Logorithmic functions }
  71.   function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }
  72.   function Log10(X: Extended): Extended;                     { Log base 10 of X}
  73.   function Log2(X: Extended): Extended;                      { Log base 2 of X }
  74.   function LogN(Base, X: Extended): Extended;                { Log base N of X }
  75.  
  76.   { Exponential functions }
  77.  
  78.   { IntPower: Raise base to an integral power.  Fast. }
  79.   function IntPower(Base: Extended; Exponent: Integer): Extended;
  80.  
  81.   { Power: Raise base to any power.
  82.     For fractional exponents, or exponents > MaxInt, base must be > 0. }
  83.   function Power(Base, Exponent: Extended): Extended;
  84.  
  85.   { Miscellaneous Routines }
  86.  
  87.   { Frexp:  Separates the mantissa and exponent of X. }
  88.   procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
  89.  
  90.   { Ldexp: returns X*2**P }
  91.   function Ldexp(X: Extended; P: Integer): Extended;
  92.  
  93.   { Ceil: Smallest integer >= X, |X| < MaxInt }
  94.   function Ceil(X: Extended):LongInt;
  95.  
  96.   { Floor: Largest integer <= X,  |X| < MaxInt }
  97.   function Floor(X: Extended): LongInt;
  98.  
  99.   { Poly: Evaluates a uniform polynomial of one variable at value X.
  100.       The coefficients are ordered in increasing powers of X:
  101.       Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
  102.   function Poly(X: Extended; const Coefficients: array of Double): Extended;
  103.  
  104.   {-----------------------------------------------------------------------
  105.   Statistical functions.
  106.  
  107.   Common commercial spreadsheet macro names for these statistical and
  108.   financial functions are given in the comments preceding each function.
  109.   -----------------------------------------------------------------------}
  110.  
  111.   { Mean:  Arithmetic average of values.  (AVG):  SUM / N }
  112.   function Mean(const Data: array of Double): Extended;
  113.  
  114.   { Sum: Sum of values.  (SUM) }
  115.   function Sum(const Data: array of Double): Extended;
  116.   function SumOfSquares(const Data: array of Double): Extended;
  117.   procedure SumsAndSquares(const Data: array of Double;
  118.     var Sum, SumOfSquares: Extended);
  119.  
  120.   { MinValue: Returns the smallest signed value in the data array (MIN) }
  121.   function MinValue(const Data: array of Double): Double;
  122.  
  123.   { MaxValue: Returns the largest signed value in the data array (MAX) }
  124.   function MaxValue(const Data: array of Double): Double;
  125.  
  126.   { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  127.   function StdDev(const Data: array of Double): Extended;
  128.  
  129.   { MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
  130.     calculating them separately.  Less accurate when the mean is very large
  131.     (> 10e7) or the variance is very small. }
  132.   procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  133.  
  134.   { Population Standard Deviation (STDP): Sqrt(PopnVariance).
  135.     Used in some business and financial calculations. }
  136.   function PopnStdDev(const Data: array of Double): Extended;
  137.  
  138.   { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
  139.   function Variance(const Data: array of Double): Extended;
  140.  
  141.   { Population Variance (VAR or VARP): TotalVariance/ N }
  142.   function PopnVariance(const Data: array of Double): Extended;
  143.  
  144.   { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
  145.   function TotalVariance(const Data: array of Double): Extended;
  146.  
  147.   { Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
  148.   function Norm(const Data: array of Double): Extended;
  149.  
  150.   { MomentSkewKurtosis: Calculates the core factors of statistical analysis:
  151.     the first four moments plus the coefficients of skewness and kurtosis.
  152.     M1 is the Mean.  M2 is the Variance.
  153.     Skew reflects symmetry of distribution: M3 / (M2**(3/2))
  154.     Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
  155.   procedure MomentSkewKurtosis(const Data: array of Double;
  156.     var M1, M2, M3, M4, Skew, Kurtosis: Extended);
  157.  
  158.   { RandG produces random numbers with Gaussian distribution about the mean.
  159.     Useful for simulating data with sampling errors. }
  160.   function RandG(Mean, StdDev: Extended): Extended;
  161.  
  162.   {-----------------------------------------------------------------------
  163.   Financial functions.  Standard set from Quattro Pro.
  164.  
  165.   Parameter conventions:
  166.  
  167.   From the point of view of A, amounts received by A are positive and
  168.   amounts disbursed by A are negative (e.g. a borrower's loan repayments
  169.   are regarded by the borrower as negative).
  170.  
  171.   Interest rates are per payment period.  11% annual percentage rate on a
  172.   loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
  173.  
  174.   -----------------------------------------------------------------------}
  175.  
  176.   type
  177.     TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);
  178.  
  179.   { Double Declining Balance (DDB) }
  180.   function DoubleDecliningBalance(Cost, Salvage: Extended;
  181.     Life, Period: Integer): Extended;
  182.  
  183.   { Future Value (FVAL) }
  184.   function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  185.     Extended; PaymentTime: TPaymentTime): Extended;
  186.  
  187.   { Interest Payment (IPAYMT)  }
  188.   function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  189.     FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  190.  
  191.   { Interest Rate (IRATE) }
  192.   function InterestRate(NPeriods: Integer;
  193.     Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  194.  
  195.   { Internal Rate of Return. (IRR) Needs array of cash flows. }
  196.   function InternalRateOfReturn(Guess: Extended;
  197.     const CashFlows: array of Double): Extended;
  198.  
  199.   { Number of Periods (NPER) }
  200.   function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  201.     PaymentTime: TPaymentTime): Extended;
  202.  
  203.   { Net Present Value. (NPV) Needs array of cash flows. }
  204.   function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  205.     PaymentTime: TPaymentTime): Extended;
  206.  
  207.   { Payment (PAYMT) }
  208.   function Payment(Rate: Extended; NPeriods: Integer;
  209.     PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  210.  
  211.   { Period Payment (PPAYMT) }
  212.   function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  213.     PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  214.  
  215.   { Present Value (PVAL) }
  216.   function PresentValue(Rate: Extended; NPeriods: Integer;
  217.     Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  218.  
  219.   { Straight Line depreciation (SLN) }
  220.   function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  221.  
  222.   { Sum-of-Years-Digits depreciation (SYD) }
  223.   function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  224.  
  225.   type
  226.     EInvalidArgument = class(EMathError) end;
  227.  
  228. implementation
  229.  
  230.   type
  231.     TPoly = record
  232.       Neg, Pos, DNeg, DPos: Extended
  233.     end;
  234.  
  235.   const
  236.     MaxIterations = 15;
  237.  
  238.   procedure ArgError;
  239.   begin
  240.     raise EInvalidArgument.Create('');
  241.   end;
  242.  
  243.   function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }
  244.   {&Frame-} {&Uses none}
  245.   asm
  246.     fld     X         // c
  247.     fld     st        // c c
  248.     fmul    st,st     // c*c c
  249.     fld1              // 1 c*c c
  250.     fsubrp  st(1),st  // (1-c*c) c
  251.     fsqrt             // s c
  252.     fxch    st(1)     // c s
  253.     fpatan
  254.   end;
  255.  
  256.   function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }
  257.   {&Frame-} {&Uses none}
  258.   asm
  259.     fld     X         // s
  260.     fld     st        // s s
  261.     fmul    st,st     // s*s s
  262.     fld1              // 1 s*s s
  263.     fsubrp  st(1),st  // (1-s*s) s
  264.     fsqrt             // c s
  265.     fpatan
  266.   end;
  267.  
  268.   { ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
  269.     IN: |Y| < 2^64, |X| < 2^64, X <> 0   OUT: [-PI..PI] radians }
  270.   function ArcTan2(Y, X: Extended): Extended;
  271.   {&Frame-} {&Uses none}
  272.   asm
  273.     fld     Y
  274.     fld     X
  275.     fpatan
  276.   end;
  277.  
  278.   { SinCos is 2x faster than calling Sin and Cos separately for the same angle }
  279.   procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
  280.   {&Frame-} {&Uses none}
  281.   asm
  282.     fld     Theta
  283.     fsincos
  284.     mov     eax, Cos
  285.     fstp    tbyte ptr [eax]
  286.     mov     eax, Sin
  287.     fstp    tbyte ptr [eax]
  288.     fwait
  289.   end;
  290.  
  291.   function Tan(X: Extended): Extended;
  292.   {&Frame-} {&Uses none}
  293.   asm
  294.     fld     X
  295.     fsincos
  296.     fdivp   st(1),st
  297.   end;
  298.  
  299.   function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }
  300.   {&Frame-} {&Uses none}
  301.   asm
  302.     fld     X
  303.     fsincos
  304.     fdivrp  st(1),st
  305.   end;
  306.  
  307.   function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }
  308.   {&Frame-} {&Uses none}
  309.   asm
  310.     fld     X
  311.     fmul    st,st
  312.     fld     Y
  313.     fmul    st,st
  314.     faddp   st(1),st
  315.     fsqrt
  316.   end;
  317.  
  318.   { Angle unit conversion routines }
  319.   function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}
  320.   begin
  321.     DegToRad := Degrees * PI / 180;
  322.   end;
  323.  
  324.   function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
  325.   begin
  326.     RadToDeg := Radians * 180 / PI;
  327.   end;
  328.  
  329.   function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
  330.   begin
  331.     GradToRad := Grads * PI / 200;
  332.   end;
  333.  
  334.   function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }
  335.   begin
  336.     RadToGrad := Radians * 200 / PI;
  337.   end;
  338.  
  339.   function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  340.   const
  341.     TwoPI = 2*PI;
  342.   begin
  343.     CycleToRad := Cycles * TwoPI
  344.   end;
  345.  
  346.   function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  347.   const
  348.     TwoPI = 2*PI;
  349.   begin
  350.     RadToCycle := Radians / TwoPI;
  351.   end;
  352.  
  353.  
  354.   { Hyperbolic functions and inverses }
  355.  
  356.   function Cosh(X: Extended): Extended;
  357.   {&Frame-} {&Uses none}
  358.   asm
  359.     fld      X
  360.     fldl2e
  361.     fmulp    st(1),st
  362.     fabs
  363.     fld      st
  364.     frndint
  365.     fsub     st(1),st
  366.     fxch     st(1)
  367.     ftst
  368.     fstsw    ax
  369.     sahf
  370.     jb       @L0
  371.     f2xm1
  372.     jmp      @L1
  373.   @L0:
  374.     fchs
  375.     f2xm1
  376.     fld1
  377.     fadd     st,st(1)
  378.     fdivp    st(1),st
  379.     fchs
  380.   @L1:
  381.     fld1
  382.     faddp    st(1),st
  383.     fscale
  384.     fstp     st(1)
  385.     fld1
  386.     fdiv     st,st(1)
  387.     faddp    st(1),st
  388.     fld1
  389.     fchs
  390.     fxch     st(1)
  391.     fscale
  392.     fstp     st(1)
  393.   end;
  394.  
  395.   function Sinh(X: Extended): Extended;
  396.   {&Frame-} {&Uses none}
  397.   asm
  398.     fld      X
  399.     ftst
  400.     fstsw    ax
  401.     mov      dx,ax
  402.     fldl2e
  403.     fmulp    st(1),st
  404.     fabs
  405.     fld      st
  406.     frndint
  407.     fsub     st(1),st
  408.     fchs
  409.     fld1
  410.     fscale
  411.     fxch     st(1)
  412.     fchs
  413.     fxch     st(2)
  414.     ftst
  415.     fstsw    ax
  416.     sahf
  417.     jb       @L0
  418.     f2xm1
  419.     jmp      @L1
  420.   @L0:
  421.     fchs
  422.     f2xm1
  423.     fld1
  424.     fadd     st,st(1)
  425.     fdivp    st(1),st
  426.     fchs
  427.   @L1:
  428.     fld1
  429.     fadd     st,st(1)
  430.     fxch     st(2)
  431.     fld1
  432.     fsubp    st(1),st
  433.     fsubp    st(1),st
  434.     fdivr    st(1),st
  435.     fxch     st(1)
  436.     fxch     st(2)
  437.     fxch     st(1)
  438.     fscale
  439.     fstp     st(1)
  440.     faddp    st(1),st
  441.     fld1
  442.     fchs
  443.     fxch     st(1)
  444.     fscale
  445.     fstp     st(1)
  446.     mov      ax,dx
  447.     sahf
  448.     jae      @L2
  449.     fchs
  450.   @L2:
  451.   end;
  452.  
  453.   function Tanh(X: Extended): Extended;
  454.   {&Frame-} {&Uses none}
  455.   asm
  456.     fld      X
  457.     ftst
  458.     fstsw    ax
  459.     mov      dx,ax
  460.     fldl2e
  461.     fmulp    st(1),st
  462.     fadd     st,st
  463.     fabs
  464.     fld      st
  465.     frndint
  466.     fsub     st(1),st
  467.     fchs
  468.     fld1
  469.     fscale
  470.     fstp     st(1)
  471.     fld1
  472.     fsubp    st(1),st
  473.     fxch     st(1)
  474.     ftst
  475.     fstsw    ax
  476.     sahf
  477.     jb       @L0
  478.     f2xm1
  479.     jmp      @L1
  480.   @L0:
  481.     fchs
  482.     f2xm1
  483.     fld1
  484.     fadd     st,st(1)
  485.     fdivp    st(1),st
  486.     fchs
  487.   @L1:
  488.     fld      st
  489.     fsub     st,st(2)
  490.     fxch     st(1)
  491.     faddp    st(2),st
  492.     fld1
  493.     fadd     st(2),st
  494.     faddp    st(2),st
  495.     fdivrp   st(1),st
  496.     mov      ax,dx
  497.     sahf
  498.     jae      @L2
  499.     fchs
  500.   @L2:
  501.   end;
  502.  
  503.   function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }
  504.   {&Frame-} {&Uses none}
  505.   asm
  506.     fld     x
  507.     fld     st        // x x
  508.     fld1              // 1 x x
  509.     fsub    st(2),st  // 1 x (x-1)
  510.     faddp   st(1),st  // (x+1) (x-1)
  511.     fld     st(1)     // (x-1) (x+1) (x-1)
  512.     fmulp   st(1),st  // (x-1)*(x+1) (x-1)
  513.     fsqrt             // sqrt((x+1)*(x-1)) (x-1)
  514.     faddp   st(1),st  // z, now ArcCosh(X) = ln(1+z)
  515.     fldln2            // loge(2) z
  516.     fxch    st(1)     // z loge(2)
  517.     fyl2xp1           // fyl2xp1(fldln2,z)
  518.   end;
  519.  
  520.   function ArcSinh(X: Extended): Extended;
  521.   {&Frame-} {&Uses none}
  522.   asm
  523.     fld     x
  524.     ftst
  525.     fstsw   ax
  526.     sahf
  527.     jne     @L0
  528.     fstp    st
  529.     fldz
  530.     jmp     @L2
  531.   @L0:
  532.     fabs               // ax
  533.     fld1               // 1 ax
  534.     fdiv    st,st(1)   // 1/ax ax
  535.     fld     st         // 1/ax 1/ax ax
  536.     fmul    st,st      // 1/ax^2 1/ax ax
  537.     fld1               // 1 1/aX^2 1/ax ax
  538.     faddp   st(1),st   // 1+1/ax^2 1/ax ax
  539.     fsqrt              // sqrt(1+1/ax^2) 1/ax ax
  540.     faddp   st(1),st   // 1/ax+sqrt(1+1/ax^2) ax
  541.     fdivr   st,st(1)
  542.     faddp   st(1),st   // z
  543.     fldln2
  544.     jae     @L1
  545.     fchs
  546.   @L1:
  547.     fxch    st(1)
  548.     fyl2xp1
  549.   @L2:
  550.   end;
  551.  
  552.   function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }
  553.   {&Frame-} {&Uses none}
  554.   asm
  555.     fld     x
  556.     ftst
  557.     fstsw   ax
  558.     sahf
  559.     fabs              // ax
  560.     fld1              // 1 ax
  561.     fsub    st,st(1)  // 1-ax ax
  562.     fdivp   st(1),st  // ax/(1-ax)
  563.     fadd    st,st     // z
  564.     fldln2
  565.     jae     @L0
  566.     fchs
  567.   @L0:
  568.     fxch    st(1)
  569.     fyl2xp1
  570.   end;
  571.  
  572.  
  573.   { Logorithmic functions }
  574.  
  575.   function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }
  576.   {&Frame-} {&Uses none}
  577.   asm
  578.     fld1
  579.     fld     X
  580.     fyl2xp1
  581.     fldln2
  582.     fmulp   st(1),st
  583.   end;
  584.  
  585.   function Log10(X: Extended): Extended;                     { Log base 10 of X}
  586.   {&Frame-} {&Uses none}
  587.   asm
  588.     fldlg2
  589.     fld      X
  590.     fyl2x
  591.   end;
  592.  
  593.   function Log2(X: Extended): Extended;                      { Log base 2 of X }
  594.   {&Frame-} {&Uses none}
  595.   asm
  596.     fld1
  597.     fld     X
  598.     fyl2x
  599.   end;
  600.  
  601.   function LogN(Base, X: Extended): Extended;                { Log base N of X }
  602.   {&Frame-} {&Uses none}
  603.   asm
  604.     fld1
  605.     fld     X
  606.     fyl2x
  607.     fld1
  608.     fld     Base
  609.     fyl2x
  610.     fdivp   st(1),st
  611.   end;
  612.  
  613.  
  614.   { Exponential functions }
  615.  
  616.   { IntPower: Raise base to an integral power.  Fast. }
  617.   function IntPower(Base: Extended; Exponent: Integer): Extended;
  618.   {&Frame-} {&Uses none}
  619.   asm
  620.     mov     eax,Exponent
  621.     cdq
  622.     xor     eax,edx
  623.     sub     eax,edx          // eax := Abs(Exponent)
  624.     mov     ecx,eax
  625.     fld     Base
  626.     fld1
  627.     jecxz   @L1
  628.   @L0:
  629.     fmul    st,st(1)
  630.     loop    @L0
  631.   @L1:
  632.     cmp     dword ptr Exponent,0
  633.     jge     @L2
  634.     fld1
  635.     fdivrp                   // Result := 1 / Result
  636.   @L2:
  637.     fstp    st(1)
  638.   end;
  639.  
  640.   { Power: Raise base to any power.
  641.     For fractional exponents, or exponents > MaxInt, base must be > 0. }
  642.   function Power(Base, Exponent: Extended): Extended;
  643.   {&Frame-} {&Uses none}
  644.   asm
  645.     fld      Exponent
  646.     fld      Base
  647.     fabs
  648.     fyl2x
  649.     fld      st
  650.     frndint
  651.     fsub     st(1),st
  652.     fxch     st(1)
  653.     ftst
  654.     fstsw    ax
  655.     sahf
  656.     jb       @L0
  657.     f2xm1
  658.     jmp      @L1
  659.   @L0:
  660.     fchs
  661.     f2xm1
  662.     fld1
  663.     fadd     st,st(1)
  664.     fdivp    st(1),st
  665.     fchs
  666.   @L1:
  667.     fld1
  668.     faddp    st(1),st
  669.     fscale
  670.     fstp     st(1)
  671.   end;
  672.  
  673.  
  674.   { Miscellaneous Routines }
  675.  
  676.   { Frexp:  Separates the mantissa and exponent of X. }
  677.   procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
  678.   {&Frame-} {&Uses none}
  679.   asm
  680.     fld     X
  681.     fxtract
  682.     mov     eax,Mantissa
  683.     fstp    tbyte ptr [eax]
  684.     mov     eax,Exponent
  685.     fistp   word ptr [eax]
  686.     fwait
  687.   end;
  688.  
  689.   { Ldexp: returns X*2**P }
  690.   function Ldexp(X: Extended; P: Integer): Extended;
  691.   {&Frame-} {&Uses none}
  692.   asm
  693.     fild    P
  694.     fld     X
  695.     fscale
  696.     fstp    st(1)
  697.   end;
  698.  
  699.   { Ceil: Smallest integer >= X, |X| < MaxInt }
  700.   function Ceil(X: Extended):LongInt;
  701.   {&Frame-} {&Uses none}
  702.   asm
  703.     fld     x
  704.     fstcw   word ptr x      { save rounding control }
  705.     fwait
  706.     mov     ax, word ptr x
  707.     and     ax, 0f3ffh
  708.     or      ax, 00800h
  709.     xchg    ax, word ptr x
  710.     fldcw   word ptr x      { set rounding toward positive infinity }
  711.     frndint                 { round x }
  712.     fistp   dword ptr x+4   { store as LongInt }
  713.     fwait
  714.     xchg    ax, word ptr x
  715.     fldcw   word ptr x      { reset rounding control }
  716.     mov     eax,dword ptr x+4
  717.   end;
  718.  
  719.   { Floor: Largest integer <= X,  |X| < MaxInt }
  720.   function Floor(X: Extended): LongInt;
  721.   {&Frame-} {&Uses none}
  722.   asm
  723.     fld     x
  724.     fstcw   word ptr x      { save rounding control }
  725.     fwait
  726.     mov     ax, word ptr x
  727.     and     ax, 0f3ffh
  728.     or      ax, 00400h
  729.     xchg    ax, word ptr x
  730.     fldcw   word ptr x      { set rounding toward negative infinity }
  731.     frndint                 { round x }
  732.     fistp   dword ptr x+4   { store as LongInt }
  733.     fwait
  734.     xchg    ax, word ptr x
  735.     fldcw   word ptr x      { reset rounding control }
  736.     mov     eax,dword ptr x+4
  737.   end;
  738.  
  739.   { Poly: Evaluates a uniform polynomial of one variable at value X.
  740.       The coefficients are ordered in increasing powers of X:
  741.       Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
  742.   function Poly(X: Extended; const Coefficients: array of Double): Extended;
  743.   {&Frame-} {&Uses none}
  744.   asm
  745.     mov     ecx,Coefficients-4  // ecx = High(Coefficients)
  746.     mov     eax,ecx
  747.     shl     eax,3
  748.     fld     X
  749.     add     eax,Coefficients
  750.     fld     qword ptr [eax]
  751.     jecxz   @L1
  752.   @L0:
  753.     sub     eax,8
  754.     fmul    st,st(1)
  755.     fadd    qword ptr [eax]
  756.     loop    @L0
  757.   @L1:
  758.     fstp    st(1)
  759.   end;
  760.  
  761.   {-----------------------------------------------------------------------
  762.   Statistical functions.
  763.  
  764.   Common commercial spreadsheet macro names for these statistical and
  765.   financial functions are given in the comments preceding each function.
  766.   -----------------------------------------------------------------------}
  767.  
  768.   { Mean:  Arithmetic average of values.  (AVG):  SUM / N }
  769.   function Mean(const Data: array of Double): Extended;
  770.   {&Frame-} {&Uses none}
  771.   asm
  772.     mov     ecx,Data-4          // ecx = High(Data)
  773.     inc     ecx
  774.     mov     Data-4,ecx
  775.     mov     eax,Data
  776.     fldz
  777.   @L0:
  778.     fadd    qword ptr [eax]
  779.     add     eax,8
  780.     loop    @L0
  781.     fidiv   dword ptr Data-4
  782.   end;
  783.  
  784.   { Sum: Sum of values.  (SUM) }
  785.   function Sum(const Data: array of Double): Extended;
  786.   {&Frame-} {&Uses none}
  787.   asm
  788.     mov     ecx,Data-4
  789.     inc     ecx
  790.     mov     eax,Data
  791.     fldz
  792.   @L0:
  793.     fadd    qword ptr [eax]
  794.     add     eax,8
  795.     loop    @L0
  796.   end;
  797.  
  798.   function SumOfSquares(const Data: array of Double): Extended;
  799.   {&Frame-} {&Uses none}
  800.   asm
  801.     mov     ecx,Data-4
  802.     inc     ecx
  803.     mov     eax,Data
  804.     fldz
  805.   @L0:
  806.     fld     qword ptr [eax]
  807.     fmul    st,st
  808.     add     eax,8
  809.     faddp   st(1),st
  810.     loop    @L0
  811.   end;
  812.  
  813.   procedure SumsAndSquares(const Data: array of Double;
  814.                            var Sum, SumOfSquares: Extended);
  815.   {&Frame-} {&Uses none}
  816.   asm
  817.     mov     ecx,Data-4
  818.     inc     ecx
  819.     mov     eax,Data
  820.     fldz
  821.     fldz
  822.   @L0:
  823.     fld     qword ptr [eax]
  824.     fadd    st(2),st
  825.     fmul    st,st
  826.     add     eax,8
  827.     faddp   st(1),st
  828.     loop    @L0
  829.     mov     eax,SumOfSquares
  830.     fstp    tbyte ptr [eax]
  831.     mov     eax,Sum
  832.     fstp    tbyte ptr [eax]
  833.     fwait
  834.   end;
  835.  
  836.   { MinValue: Returns the smallest signed value in the data array (MIN) }
  837.   function MinValue(const Data: array of Double): Double;
  838.   {&Frame-} {&Uses ebx}
  839.   asm
  840.     mov     ebx,Data
  841.     mov     ecx,Data-4
  842.     fld     qword ptr [ebx]
  843.     jecxz   @L2
  844.   @L0:
  845.     add     ebx,8
  846.     fcom    qword ptr [ebx]
  847.     fstsw   ax
  848.     sahf
  849.     jbe     @L1
  850.     fstp    st
  851.     fld     qword ptr [ebx]
  852.   @L1:
  853.     loop    @L0
  854.   @L2:
  855.   end;
  856.  
  857.   { MaxValue: Returns the largest signed value in the data array (MAX) }
  858.   function MaxValue(const Data: array of Double): Double;
  859.   {&Frame-} {&Uses ebx}
  860.   asm
  861.     push    ebx
  862.     mov     ebx,Data
  863.     mov     ecx,Data-4
  864.     fld     qword ptr [ebx]
  865.     jecxz   @L2
  866.   @L0:
  867.     add     ebx,8
  868.     fcom    qword ptr [ebx]
  869.     fstsw   ax
  870.     sahf
  871.     jae     @L1
  872.     fstp    st
  873.     fld     qword ptr [ebx]
  874.   @L1:
  875.     loop    @L0
  876.   @L2:
  877.     pop     ebx
  878.   end;
  879.  
  880.   { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  881.   function StdDev(const Data: array of Double): Extended;
  882.   {&Frame-} {&Uses none}
  883.   asm
  884.     mov     eax,Data
  885.     mov     ecx,Data-4
  886.     inc     ecx
  887.     mov     Data-4,ecx
  888.     fldz
  889.   @L0:
  890.     fadd    qword ptr [eax]
  891.     add     eax,8
  892.     loop    @L0
  893.     fidiv   dword ptr Data-4 // xb = sx/N
  894.     mov     eax,Data
  895.     mov     ecx,Data-4
  896.     fldz
  897.   @L1:
  898.     fld     qword ptr [eax]
  899.     fsub    st,st(2)
  900.     fmul    st,st
  901.     faddp   st(1),st
  902.     add     eax,8
  903.     loop @L1
  904.     fstp    st(1)
  905.     dec     dword ptr Data-4
  906.     fidiv   dword ptr Data-4
  907.     fsqrt
  908.   end;
  909.  
  910.   { MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
  911.     calculating them separately.  Less accurate when the mean is very large
  912.     (> 10e7) or the variance is very small. }
  913.   procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  914.   {&Frame-} {&Uses none}
  915.   asm
  916.     mov     eax,Data
  917.     mov     ecx,Data-4
  918.     inc     ecx
  919.     mov     Data-4,ecx
  920.     fldz
  921.   @L0:
  922.     fadd    qword ptr [eax]
  923.     add     eax,8
  924.     loop    @L0
  925.     fidiv   dword ptr Data-4
  926.     mov     eax,Mean
  927.     fld     st
  928.     fstp    tbyte ptr [eax]
  929.     mov     eax,Data
  930.     mov     ecx,Data-4
  931.     fldz
  932.   @L1:
  933.     fld     qword ptr [eax]
  934.     fsub    st,st(2)
  935.     fmul    st,st
  936.     faddp   st(1),st
  937.     add     eax,8
  938.     loop @L1
  939.     fstp    st(1)
  940.     dec     dword ptr Data-4
  941.     fidiv   dword ptr Data-4
  942.     fsqrt
  943.     mov     eax,StdDev
  944.     fstp    tbyte ptr [eax]
  945.     fwait
  946.   end;
  947.  
  948.   { Population Standard Deviation (STDP): Sqrt(PopnVariance).
  949.     Used in some business and financial calculations. }
  950.   function PopnStdDev(const Data: array of Double): Extended;
  951.   {&Frame-} {&Uses none}
  952.   asm
  953.     mov     eax,Data
  954.     mov     ecx,Data-4
  955.     inc     ecx
  956.     mov     Data-4,ecx
  957.     fldz
  958.   @L0:
  959.     fadd    qword ptr [eax]
  960.     add     eax,8
  961.     loop    @L0
  962.     fidiv   dword ptr Data-4 // xb = sx/N
  963.     mov     eax,Data
  964.     mov     ecx,Data-4
  965.     fldz
  966.   @L1:
  967.     fld     qword ptr [eax]
  968.     fsub    st,st(2)
  969.     fmul    st,st
  970.     faddp   st(1),st
  971.     add     eax,8
  972.     loop @L1
  973.     fstp    st(1)
  974.     fidiv   dword ptr Data-4
  975.     fsqrt
  976.   end;
  977.  
  978.   { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
  979.   function Variance(const Data: array of Double): Extended;
  980.   {&Frame-} {&Uses none}
  981.   asm
  982.     mov     eax,Data
  983.     mov     ecx,Data-4
  984.     inc     ecx
  985.     mov     Data-4,ecx
  986.     fldz
  987.   @L0:
  988.     fadd    qword ptr [eax]
  989.     add     eax,8
  990.     loop    @L0
  991.     fidiv   dword ptr Data-4 // xb = sx/N
  992.     mov     eax,Data
  993.     mov     ecx,Data-4
  994.     fldz
  995.   @L1:
  996.     fld     qword ptr [eax]
  997.     fsub    st,st(2)
  998.     fmul    st,st
  999.     faddp   st(1),st
  1000.     add     eax,8
  1001.     loop @L1
  1002.     fstp    st(1)
  1003.     dec     dword ptr Data-4
  1004.     fidiv   dword ptr Data-4
  1005.   end;
  1006.  
  1007.   { Population Variance (VAR or VARP): TotalVariance/ N }
  1008.   function PopnVariance(const Data: array of Double): Extended;
  1009.   {&Frame-} {&Uses none}
  1010.   asm
  1011.     mov     eax,Data
  1012.     mov     ecx,Data-4
  1013.     inc     ecx
  1014.     mov     Data-4,ecx
  1015.     fldz
  1016.   @L0:
  1017.     fadd    qword ptr [eax]
  1018.     add     eax,8
  1019.     loop    @L0
  1020.     fidiv   dword ptr Data-4 // xb = sx/N
  1021.     mov     eax,Data
  1022.     mov     ecx,Data-4
  1023.     fldz
  1024.   @L1:
  1025.     fld     qword ptr [eax]
  1026.     fsub    st,st(2)
  1027.     fmul    st,st
  1028.     faddp   st(1),st
  1029.     add     eax,8
  1030.     loop @L1
  1031.     fstp    st(1)
  1032.     fidiv   dword ptr Data-4
  1033.     fwait
  1034.   end;
  1035.  
  1036.   { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
  1037.   function TotalVariance(const Data: array of Double): Extended;
  1038.   {&Frame-} {&Uses none}
  1039.   asm
  1040.     mov     eax,Data
  1041.     mov     ecx,Data-4
  1042.     inc     ecx
  1043.     mov     Data-4,ecx
  1044.     fldz
  1045.   @L0:
  1046.     fadd    qword ptr [eax]
  1047.     add     eax,8
  1048.     loop    @L0
  1049.     fidiv   dword ptr Data-4 // xb = sx/N
  1050.     mov     eax,Data
  1051.     mov     ecx,Data-4
  1052.     fldz
  1053.   @L1:
  1054.     fld     qword ptr [eax]
  1055.     fsub    st,st(2)
  1056.     fmul    st,st
  1057.     faddp   st(1),st
  1058.     add     eax,8
  1059.     loop @L1
  1060.     fstp    st(1)
  1061.   end;
  1062.  
  1063.   { Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
  1064.   function Norm(const Data: array of Double): Extended;
  1065.   {&Frame-} {&Uses none}
  1066.   asm
  1067.     mov     eax,Data
  1068.     fldz
  1069.     mov     ecx,Data-4
  1070.     inc     ecx
  1071.   @L0:
  1072.     fld     qword ptr [eax]
  1073.     fmul    st,st
  1074.     add     eax,8
  1075.     faddp   st(1),st
  1076.     loop    @L0
  1077.     fsqrt
  1078.   end;
  1079.  
  1080.   { MomentSkewKurtosis: Calculates the core factors of statistical analysis:
  1081.     the first four moments plus the coefficients of skewness and kurtosis.
  1082.     M1 is the Mean.  M2 is the Variance.
  1083.     Skew reflects symmetry of distribution: M3 / (M2**(3/2))
  1084.     Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
  1085.   procedure MomentSkewKurtosis(const Data: array of Double;
  1086.                                var M1, M2, M3, M4,
  1087.                                Skew, Kurtosis: Extended);
  1088.   {&Frame-} {&Uses none}
  1089.   asm
  1090.     mov     ecx,Data-4
  1091.     inc     ecx
  1092.     // calculate the mean
  1093.     mov     Data-4,ecx
  1094.     mov     eax,Data
  1095.     fldz
  1096.   @L0:
  1097.     fadd    qword ptr [eax]
  1098.     add     eax,8
  1099.     loop    @L0
  1100.     fidiv   dword ptr Data-4    // avg
  1101.     // now calculate first (absolute), second, third and fourth moments
  1102.     // of the deviation from the mean
  1103.     mov     ecx,Data-4
  1104.     mov     eax,Data
  1105.     fldz
  1106.     fldz
  1107.     fldz
  1108.     fldz                        // M1 M2 M3 M4 avg
  1109.   @L1:
  1110.     fld     qword ptr [eax]     // x M1 M2 M3 M4 avg
  1111.     fsub    st,st(5)            // dev M1 M2 M3 M4 avg
  1112.     fld     st                  // dev dev M1 M2 M3 M4 avg
  1113.     fabs                        // adev dev M1 M2 M3 M4 avg
  1114.     fadd    st(2),st            // adev dev M1+ M2 M3 M4 avg
  1115.     fmul    st,st               // dev^2 dev M1+ M2 M3 M4 avg
  1116.     fmul    st(1),st            // dev^2 dev^3 M1+ M2 M3 M4 avg
  1117.     fadd    st(3),st            // dev^2 dev^3 M1+ M2+ M3 M4 avg
  1118.     fmul    st,st               // dev^4 dev^3 M1+ M2+ M3 M4 avg
  1119.     faddp   st(5),st            // dev^3 M1+ M2+ M3 M4+ avg
  1120.     faddp   st(3),st            // M1+ M2+ M3+ M4+ avg
  1121.     add     eax,8
  1122.     loop    @L1
  1123.     mov     eax,M1
  1124.     fstp    tbyte ptr [eax]     // M2 M3 M4 avg
  1125.     mov     eax,M2
  1126.     fstp    tbyte ptr [eax]     // M3 M4 avg
  1127.     mov     eax,M3
  1128.     fstp    tbyte ptr [eax]     // M4 avg
  1129.     mov     eax,M4
  1130.     fstp    tbyte ptr [eax]     // avg
  1131.     fstp    st
  1132.     mov     eax,M2
  1133.     fld     tbyte ptr [eax]
  1134.     dec     dword ptr Data-4
  1135.     fidiv   dword ptr Data-4    // var
  1136.     ftst                        // var = 0?
  1137.     fstsw   ax
  1138.     sahf
  1139.     je      @Err
  1140.     fld     st                  // var var
  1141.     fsqrt                       // sdev var
  1142.     mov     eax,M3
  1143.     fld     tbyte ptr [eax]     // M3 sdev var
  1144.     inc     dword ptr Data-4
  1145.     fidiv   dword ptr Data-4    // M3/N sdev var
  1146.     fdiv    st,st(1)            // M3/(N*sdev) var
  1147.     fdiv    st,st(2)            // Skew=M3/(N*sdev*var) var
  1148.     mov     eax,Skew
  1149.     fstp    tbyte ptr [eax]     // var
  1150.     fmul    st,st               // var*var
  1151.     mov     eax,M4
  1152.     fld     tbyte ptr [eax]     // M4 var*var
  1153.     fidiv   dword ptr Data-4    // M4/N var*var
  1154.     fdivrp  st(1),st            // M4/(N*var*var)
  1155.     mov     dword ptr Data-4,3
  1156.     fisub   dword ptr Data-4    // Kurtosis = M4/(N*var*var) - 3
  1157.     mov     eax,Kurtosis
  1158.     fstp    tbyte ptr [eax]
  1159.     jmp     @Done
  1160.   @Err:
  1161.     // no skew/kurtosis when var=0
  1162.     fldz
  1163.     mov     eax,Skew
  1164.     fstp    tbyte ptr [eax]
  1165.     fldz
  1166.     mov     eax,Kurtosis
  1167.     fstp    tbyte ptr [eax]
  1168.   @Done:
  1169.   end;
  1170.  
  1171.   { RandG produces random numbers with Gaussian distribution about the mean.
  1172.     Useful for simulating data with sampling errors. }
  1173.   function RandG(Mean, StdDev: Extended): Extended;
  1174.   var
  1175.     v1,v2,rsq,fac : Extended;
  1176.   begin
  1177.     repeat
  1178.       v1 := 2.0 * random - 1.0;
  1179.       v2 := 2.0 * random - 1.0;
  1180.       rsq := sqr(v1) + sqr(v2);
  1181.     until (rsq < 1.0) and (rsq <> 0.0);
  1182.     fac := sqrt(-2.0*ln(rsq)/rsq);
  1183.     RandG := Mean + StdDev*v2*fac;
  1184.   end;
  1185.  
  1186.   {-----------------------------------------------------------------------
  1187.   Financial functions.  Standard set from Quattro Pro.
  1188.  
  1189.   Parameter conventions:
  1190.  
  1191.   From the point of view of A, amounts received by A are positive and
  1192.   amounts disbursed by A are negative (e.g. a borrower's loan repayments
  1193.   are regarded by the borrower as negative).
  1194.  
  1195.   Interest rates are per payment period.  11% annual percentage rate on a
  1196.   loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
  1197.  
  1198.   -----------------------------------------------------------------------}
  1199.  
  1200.   { Double Declining Balance (DDB) }
  1201.   function DoubleDecliningBalance(Cost, Salvage: Extended;
  1202.     Life, Period: Integer): Extended;
  1203.   VAR
  1204.     i : Integer;
  1205.     c,DDB : Extended;
  1206.   begin
  1207.     DDB := 0.0;
  1208.     for i := 1 TO Period do begin
  1209.       c := Cost*(1.0 - 2.0/Life);
  1210.       if c < Salvage then c := Salvage;
  1211.       DDB := Cost - c;
  1212.       Cost := c;
  1213.     end;
  1214.     DoubleDecliningBalance := DDB;
  1215.   end;
  1216.  
  1217.   { Future Value (FVAL) }
  1218.   function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  1219.     Extended; PaymentTime: TPaymentTime): Extended;
  1220.   VAR
  1221.     r, rN, rN1, s, p : Extended;
  1222.   begin
  1223.     r   := 1.0 + Rate;
  1224.     rN  := power(r,NPeriods);
  1225.     rN1 := r*rN;
  1226.     s   := -Payment*(rN1 - r)/(r - 1.0);
  1227.     if PaymentTime = ptEndOfPeriod then
  1228.       s := s/r;
  1229.     p   := -PresentValue*rN;
  1230.     FutureValue := s + p;
  1231.   end;
  1232.  
  1233.   { Interest Payment (IPAYMT)  }
  1234.   function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  1235.     FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1236.   var
  1237.     rf,r,pf,s,PayOf,PayOf0 : Extended;
  1238.  
  1239.     function FofPn(n: Integer) : Extended;
  1240.     begin
  1241.       FofPn := PresentValue*power(r,n);
  1242.     end;
  1243.  
  1244.     function FofAn(n: Integer) : Extended;
  1245.     begin
  1246.       result := pf*(power(r,n)-1.0)/rf;
  1247.       if PaymentTime = ptStartOfPeriod then
  1248.         result := r * result;
  1249.     end;
  1250.  
  1251.   begin
  1252.     rf  := Rate;
  1253.     r   := 1.0+rf;
  1254.     // calculate monthly payment pf
  1255.     s   := power(r,NPeriods);
  1256.     pf  := PresentValue*rf*(s/(s-1.0));
  1257.     // calculate outstanding principle
  1258.     PayOf := FofPn(Period) - FofAn(Period);
  1259.     // calculate previous outstanding principle
  1260.     PayOf0 := FofPn(Period-1) - FofAn(Period-1);
  1261.     InterestPayment := PayOf0 - PayOf;
  1262.   end;
  1263.  
  1264.   { Interest Rate (IRATE) }
  1265.   function InterestRate(NPeriods: Integer;
  1266.     Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1267.  
  1268.     function e(r: Extended) : Extended;
  1269.     var rN, s0, s1, s2 : Extended;
  1270.     begin
  1271.       rN := power(r,NPeriods);
  1272.       if PaymentTime = ptStartOfPeriod then begin
  1273.         s0 := rN*(PresentValue + Payment);
  1274.         s2 := 0.0;
  1275.       end
  1276.       else begin
  1277.         s0 := rN*PresentValue;
  1278.         s2 := Payment;
  1279.       end;
  1280.       if r = 1.0 then
  1281.         s1 := Payment*(NPeriods*rN/r - 1)
  1282.       else
  1283.         s1 := Payment*(rN-r)/(r-1);
  1284.       e := FutureValue + s0 + s1 + s2;
  1285.     end;
  1286.  
  1287.   var
  1288.     r0, r1, dr, e0, e1, de, eok : Extended;
  1289.     count : LongInt;
  1290.   begin
  1291.     if abs(FutureValue) > abs(PresentValue) then
  1292.       eok := abs(FutureValue)
  1293.     else
  1294.       eok := abs(PresentValue);
  1295.     r1 := 1.001;
  1296.     e1 := e(r1);
  1297.     count := 0;
  1298.     repeat
  1299.       e0 := e1;
  1300.       r0 := r1;
  1301.       de := ( e(r0+0.000001) - e0 ) / 0.000001;
  1302.       dr := -e0/de;
  1303.       if dr > 0.05 then dr := 0.05;
  1304.       if dr < -0.05 then dr := -0.05;
  1305.       e1 := e(r0+dr);
  1306.       while abs(e1) > abs(e0) do begin
  1307.         dr := 0.5*dr;
  1308.         e1 := e(r0+dr);
  1309.       end;
  1310.       r1 := r0 + dr;
  1311.       inc(count);
  1312.     until (abs(e1) < 1.0e-14*eok) or (count > 50);
  1313.     InterestRate := r1 - 1.0;
  1314.   end;
  1315.  
  1316.   { Discounted cash flow functions. }
  1317.  
  1318.   function Compound(R: Extended; N: Integer): Extended;
  1319.   { Return (1 + R)**N. }
  1320.   begin
  1321.     Result := IntPower(1.0 + R, N)
  1322.   end;
  1323.  
  1324.   function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
  1325.     var CompoundRN: Extended): Extended;
  1326.   { Set CompoundRN to Compound(R, N),
  1327.     return (1 + Rate*PaymentTime) * (1 - Compound(R, N)) / R }
  1328.   begin
  1329.     if R = 0.0 then
  1330.     begin
  1331.       CompoundRN := 1.0;
  1332.       Result := -N
  1333.     end
  1334.     else
  1335.     begin
  1336.       { 6.1E-5 approx= 2**-14 }
  1337.       if Abs(R) < 6.1E-5 then
  1338.       begin
  1339.         CompoundRN := Exp(N * LnXP1(R));
  1340.         Result := (CompoundRN - 1.0) / -R;
  1341.       end
  1342.       else
  1343.       begin
  1344.         CompoundRN := Compound(R, N);
  1345.         Result := (1 - CompoundRN) / R
  1346.       end;
  1347.   //    if PaymentTime = ptEndOfPeriod then Result := Result + R
  1348.       if PaymentTime = ptStartOfPeriod then Result := Result * (1 + R);
  1349.     end
  1350.   end;
  1351.  
  1352.   procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);
  1353.   { Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.
  1354.     Accumulate positive and negative terms separately. }
  1355.   var
  1356.     I: Integer;
  1357.     Neg, Pos, DNeg, DPos: Extended;
  1358.   begin
  1359.     Neg := 0.0;
  1360.     Pos := 0.0;
  1361.     DNeg := 0.0;
  1362.     DPos := 0.0;
  1363.     for I := Low(A) to High(A) do
  1364.     begin
  1365.       DNeg := X * DNeg + Neg;
  1366.       Neg := Neg * X;
  1367.       DPos := X * DPos + Pos;
  1368.       Pos := Pos * X;
  1369.       if A[I] >= 0.0 then Pos := Pos + A[I] else Neg := Neg + A[I]
  1370.     end;
  1371.     Poly.Neg := Neg;
  1372.     Poly.Pos := Pos;
  1373.     Poly.DNeg := DNeg * X;
  1374.     Poly.DPos := DPos * X;
  1375.   end;
  1376.  
  1377.   function RelSmall(X, Y: Extended): Boolean;
  1378.   { Returns True if X is small relative to Y }
  1379.   const
  1380.     C1: Double = 1E-15;
  1381.     C2: Double = 1E-12;
  1382.   begin
  1383.     Result := Abs(X) < (C1 + C2 * Abs(Y))
  1384.   end;
  1385.  
  1386.   function InternalRateOfReturn(Guess: Extended; const CashFlows: array of Double): Extended;
  1387.   {
  1388.   Use Newton's method to solve NPV = 0, where NPV is a polynomial in
  1389.   x = 1/(1+rate).  Split the coefficients into negative and postive sets:
  1390.     neg + pos = 0, so pos = -neg, so  -neg/pos = 1
  1391.   Then solve:
  1392.     log(-neg/pos) = 0
  1393.  
  1394.     Let  t = log(1/(1+r) = -LnXP1(r)
  1395.     then r = exp(-t) - 1
  1396.   Iterate on t, then use the last equation to compute r.
  1397.   }
  1398.   var
  1399.     T, Y: Extended;
  1400.     Poly: TPoly;
  1401.     K, Count: Integer;
  1402.  
  1403.     function ConditionP(const CashFlows: array of Double): Integer;
  1404.     { Guarantees existence and uniqueness of root.  The sign of payments
  1405.       must change exactly once, the net payout must be always > 0 for
  1406.       first portion, then each payment must be >= 0.
  1407.       Returns: 0 if condition not satisfied, > 0 if condition satisfied
  1408.       and this is the index of the first value considered a payback. }
  1409.     var
  1410.       X: Double;
  1411.       I, K: Integer;
  1412.     begin
  1413.       K := High(CashFlows);
  1414.       while (K >= 0) and (CashFlows[K] >= 0.0) do Dec(K);
  1415.       Inc(K);
  1416.       if K > 0 then begin
  1417.         X := 0.0;
  1418.         I := 0;
  1419.         while I < K do begin
  1420.           X := X + CashFlows[I];
  1421.           if X >= 0.0 then begin
  1422.             K := 0;
  1423.             Break
  1424.           end;
  1425.           Inc(I)
  1426.         end
  1427.       end;
  1428.       Result := K
  1429.     end;
  1430.  
  1431.   begin
  1432.     Result := 0;
  1433.     K := ConditionP(CashFlows);
  1434.     if K < 0 then ArgError;
  1435.     if K = 0 then
  1436.     begin
  1437.       if Guess <= -1.0 then ArgError;
  1438.       T := -LnXP1(Guess)
  1439.     end else
  1440.       T := 0.0;
  1441.     for Count := 1 to MaxIterations do
  1442.     begin
  1443.       PolyX(CashFlows, Exp(T), Poly);
  1444.       if Poly.Pos <= Poly.Neg then ArgError;
  1445.       if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then
  1446.       begin
  1447.         Result := -1.0;
  1448.         Exit;
  1449.       end;
  1450.       with Poly do
  1451.         Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);
  1452.       T := T - Y;
  1453.       if RelSmall(Y, T) then
  1454.       begin
  1455.         Result := Exp(-T) - 1.0;
  1456.         Exit;
  1457.       end
  1458.     end;
  1459.     ArgError;
  1460.   end;
  1461.  
  1462.   { Number of Periods (NPER) }
  1463.   function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  1464.     PaymentTime: TPaymentTime): Extended;
  1465.   { If Rate = 0 then nper := -(pv + fv) / pmt
  1466.     else t := pv + pmt * (1 + rate*pmtTime) / rate
  1467.          nper := LnXP1(-(pv + vf) / t) / LnXP1(rate)
  1468.   }
  1469.   var
  1470.     T: Extended;
  1471.   begin
  1472.     if Frac(Rate) = 0.0 then
  1473.       Result := -(PresentValue + FutureValue) / Payment
  1474.     else
  1475.     begin
  1476.       T := PresentValue + Payment * (1 + Integer(PaymentTime) * Rate) / Rate;
  1477.       Result := LnXP1(-(PresentValue + FutureValue) / T) / LnXP1(Rate)
  1478.     end
  1479.   end;
  1480.  
  1481.   { Net Present Value. (NPV) Needs array of cash flows. }
  1482.   function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  1483.     PaymentTime: TPaymentTime): Extended;
  1484.   var
  1485.     X, Y: Extended;
  1486.     I: Integer;
  1487.   begin
  1488.     if Rate <= -1.0 then ArgError;
  1489.     X := 1.0 / (1.0 + Rate);
  1490.     Y := CashFlows[High(CashFlows)];
  1491.     for I := High(CashFlows) - 1 downto Low(CashFlows) do
  1492.       Y := X * Y + CashFlows[I];
  1493.     if PaymentTime = ptStartOfPeriod then Y := X * Y;
  1494.     Result := Y;
  1495.   end;
  1496.  
  1497.   { Payment (PAYMT) }
  1498.   function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,
  1499.     FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):
  1500.     Extended;
  1501.   var
  1502.     Pmt, Rate1, T1, T2, T3: Extended;
  1503.     M: Integer;
  1504.  
  1505.     function Annuity(R: Extended; N: Integer): Extended;
  1506.     begin
  1507.       Result := (1 - Compound(R, N)) / R
  1508.     end;
  1509.  
  1510.   begin
  1511.     M := Period;
  1512.     Rate1 := 0.0;
  1513.     if PaymentTime = ptEndOfPeriod then
  1514.     begin
  1515.       Inc(M);
  1516.       Rate1 := Rate
  1517.     end;
  1518.     T1 := (1 + Rate1) * Annuity(Rate, NPeriods);
  1519.     T2 := (1 + Rate1) * Annuity(Rate, M);
  1520.     T3 := (1 + Rate1) * Annuity(Rate, NPeriods - M);
  1521.     Pmt := (FutureValue + PresentValue * Compound(Rate, NPeriods)) / T1;
  1522.     IntPmt := Rate *
  1523.               (FutureValue * Compound(Rate, M) - PresentValue * T2 * T3) / T1;
  1524.     Result := Pmt - IntPmt
  1525.   end;
  1526.  
  1527.   function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:
  1528.     Extended; PaymentTime: TPaymentTime): Extended;
  1529.   var
  1530.     Annuity, CompoundRN: Extended;
  1531.   begin
  1532.     Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1533.     if CompoundRN > 1.0E16 then
  1534.       Result := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)
  1535.     else
  1536.       Result := (PresentValue * CompoundRN + FutureValue) / Annuity
  1537.   end;
  1538.  
  1539.   { Period Payment (PPAYMT) }
  1540.   function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  1541.     PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1542.   var
  1543.     Junk: Extended;
  1544.   begin
  1545.     if (Period < 1) or (Period > NPeriods) then ArgError;
  1546.     Result := PaymentParts(Period, NPeriods, Rate, PresentValue, FutureValue,
  1547.                            PaymentTime, Junk);
  1548.   end;
  1549.  
  1550.   { Present Value (PVAL) }
  1551.   function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:
  1552.     Extended; PaymentTime: TPaymentTime): Extended;
  1553.   var
  1554.     Annuity, CompoundRN: Extended;
  1555.   begin
  1556.     Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1557.     if CompoundRN > 1.0E16 then
  1558.       Result := -(Payment / Rate * Integer(PaymentTime) * Payment)
  1559.     else
  1560.       Result := (Payment * Annuity - FutureValue) / CompoundRN
  1561.   end;
  1562.  
  1563.   { Straight Line depreciation (SLN) }
  1564.   function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  1565.   { Spreads depreciation linearly over life. }
  1566.   begin
  1567.     if Life < 1 then ArgError;
  1568.     Result := (Cost - Salvage) / Life
  1569.   end;
  1570.  
  1571.   { Sum-of-Years-Digits depreciation (SYD) }
  1572.   function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  1573.   { SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }
  1574.   var
  1575.     X1, X2: Extended;
  1576.   begin
  1577.     if (Period < 1) or (Life < Period) then ArgError;
  1578.     X1 := 2 * (Life - Period + 1);
  1579.     X2 := Life * (Life + 1);
  1580.     Result := (Cost - Salvage) * X1 / X2
  1581.   end;
  1582.  
  1583. end.
  1584.