home *** CD-ROM | disk | FTP | other *** search
/ PC Plus SuperCD (UK) 2000 March / pcp161b.iso / full / delphi / RUNIMAGE / DELPHI30 / SOURCE / RTL / SYS / MATH.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1997-08-03  |  42.5 KB  |  1,459 lines

  1.  
  2. {*******************************************************}
  3. {                                                       }
  4. {       Delphi Runtime Library                          }
  5. {       Math Unit                                       }
  6. {                                                       }
  7. {       Copyright (C) 1996,97 Borland International     }
  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 Delphi language or System unit. }
  16.  
  17. {$N+,S-}
  18.  
  19. interface
  20.  
  21. uses SysUtils;
  22.  
  23. const   { Ranges of the IEEE floating point types, including denormals }
  24.   MinSingle        =     1.5e-45;
  25.   MaxSingle        =     3.4e+38;
  26.   MinDouble        =     5.0e-324;
  27.   MaxDouble        =     1.7e+308;
  28.   MinExtended      =     3.4e-4932;
  29.   MaxExtended      =     1.1e+4932;
  30.   MinComp          =     -9.223372036854775807e+18;
  31.   MaxComp          =     9.223372036854775807e+18;
  32.  
  33. {-----------------------------------------------------------------------
  34. References:
  35.  
  36. 1) P.J. Plauger, "The Standard C Library", Prentice-Hall, 1992, Ch. 7.
  37. 2) W.J. Cody, Jr., and W. Waite, "Software Manual For the Elementary
  38.    Functions", Prentice-Hall, 1980.
  39. 3) Namir Shammas, "C/C++ Mathematical Algorithms for Scientists and Engineers",
  40.    McGraw-Hill, 1995, Ch 8.
  41. 4) H.T. Lau, "A Numerical Library in C for Scientists and Engineers",
  42.    CRC Press, 1994, Ch. 6.
  43. 5) "Pentium(tm) Processor User's Manual, Volume 3: Architecture
  44.    and Programming Manual", Intel, 1994
  45.  
  46. All angle parameters and results of trig functions are in radians.
  47.  
  48. Most of the following trig and log routines map directly to Intel 80387 FPU
  49. floating point machine instructions.  Input domains, output ranges, and
  50. error handling are determined largely by the FPU hardware.
  51. Routines coded in assembler favor the Pentium FPU pipeline architecture.
  52. -----------------------------------------------------------------------}
  53.  
  54. { Trigonometric functions }
  55. function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }
  56. function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }
  57.  
  58. { ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
  59.   IN: |Y| < 2^64, |X| < 2^64, X <> 0   OUT: [-PI..PI] radians }
  60. function ArcTan2(Y, X: Extended): Extended;
  61.  
  62. { SinCos is 2x faster than calling Sin and Cos separately for the same angle }
  63. procedure SinCos(Theta: Extended; var Sin, Cos: Extended) register;
  64. function Tan(X: Extended): Extended;
  65. function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }
  66. function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }
  67.  
  68. { Angle unit conversion routines }
  69. function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}
  70. function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
  71. function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
  72. function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }
  73. function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  74. function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  75.  
  76. { Hyperbolic functions and inverses }
  77. function Cosh(X: Extended): Extended;
  78. function Sinh(X: Extended): Extended;
  79. function Tanh(X: Extended): Extended;
  80. function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }
  81. function ArcSinh(X: Extended): Extended;
  82. function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }
  83.  
  84. { Logorithmic functions }
  85. function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }
  86. function Log10(X: Extended): Extended;                     { Log base 10 of X}
  87. function Log2(X: Extended): Extended;                      { Log base 2 of X }
  88. function LogN(Base, X: Extended): Extended;                { Log base N of X }
  89.  
  90. { Exponential functions }
  91.  
  92. { IntPower: Raise base to an integral power.  Fast. }
  93. function IntPower(Base: Extended; Exponent: Integer): Extended register;
  94.  
  95. { Power: Raise base to any power.
  96.   For fractional exponents, or |exponents| > MaxInt, base must be > 0. }
  97. function Power(Base, Exponent: Extended): Extended;
  98.  
  99.  
  100. { Miscellaneous Routines }
  101.  
  102. { Frexp:  Separates the mantissa and exponent of X. }
  103. procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer) register;
  104.  
  105. { Ldexp: returns X*2**P }
  106. function Ldexp(X: Extended; P: Integer): Extended register;
  107.  
  108. { Ceil: Smallest integer >= X, |X| < MaxInt }
  109. function Ceil(X: Extended):Integer;
  110.  
  111. { Floor: Largest integer <= X,  |X| < MaxInt }
  112. function Floor(X: Extended): Integer;
  113.  
  114. { Poly: Evaluates a uniform polynomial of one variable at value X.
  115.     The coefficients are ordered in increasing powers of X:
  116.     Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
  117. function Poly(X: Extended; const Coefficients: array of Double): Extended;
  118.  
  119. {-----------------------------------------------------------------------
  120. Statistical functions.
  121.  
  122. Common commercial spreadsheet macro names for these statistical and
  123. financial functions are given in the comments preceding each function.
  124. -----------------------------------------------------------------------}
  125.  
  126. { Mean:  Arithmetic average of values.  (AVG):  SUM / N }
  127. function Mean(const Data: array of Double): Extended;
  128.  
  129. { Sum: Sum of values.  (SUM) }
  130. function Sum(const Data: array of Double): Extended register;
  131. function SumInt(const Data: array of Integer): Integer register;
  132. function SumOfSquares(const Data: array of Double): Extended;
  133. procedure SumsAndSquares(const Data: array of Double;
  134.   var Sum, SumOfSquares: Extended) register;
  135.  
  136. { MinValue: Returns the smallest signed value in the data array (MIN) }
  137. function MinValue(const Data: array of Double): Double;
  138. function MinIntValue(const Data: array of Integer): Integer;
  139.  
  140. { MaxValue: Returns the largest signed value in the data array (MAX) }
  141. function MaxValue(const Data: array of Double): Double;
  142. function MaxIntValue(const Data: array of Integer): Integer;
  143.  
  144. { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  145. function StdDev(const Data: array of Double): Extended;
  146.  
  147. { MeanAndStdDev calculates Mean and StdDev in one call. }
  148. procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  149.  
  150. { Population Standard Deviation (STDP): Sqrt(PopnVariance).
  151.   Used in some business and financial calculations. }
  152. function PopnStdDev(const Data: array of Double): Extended;
  153.  
  154. { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
  155. function Variance(const Data: array of Double): Extended;
  156.  
  157. { Population Variance (VAR or VARP): TotalVariance/ N }
  158. function PopnVariance(const Data: array of Double): Extended;
  159.  
  160. { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
  161. function TotalVariance(const Data: array of Double): Extended;
  162.  
  163. { Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
  164. function Norm(const Data: array of Double): Extended;
  165.  
  166. { MomentSkewKurtosis: Calculates the core factors of statistical analysis:
  167.   the first four moments plus the coefficients of skewness and kurtosis.
  168.   M1 is the Mean.  M2 is the Variance.
  169.   Skew reflects symmetry of distribution: M3 / (M2**(3/2))
  170.   Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
  171. procedure MomentSkewKurtosis(const Data: array of Double;
  172.   var M1, M2, M3, M4, Skew, Kurtosis: Extended);
  173.  
  174. { RandG produces random numbers with Gaussian distribution about the mean.
  175.   Useful for simulating data with sampling errors. }
  176. function RandG(Mean, StdDev: Extended): Extended;
  177.  
  178. {-----------------------------------------------------------------------
  179. Financial functions.  Standard set from Quattro Pro.
  180.  
  181. Parameter conventions:
  182.  
  183. From the point of view of A, amounts received by A are positive and
  184. amounts disbursed by A are negative (e.g. a borrower's loan repayments
  185. are regarded by the borrower as negative).
  186.  
  187. Interest rates are per payment period.  11% annual percentage rate on a
  188. loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
  189.  
  190. -----------------------------------------------------------------------}
  191.  
  192. type
  193.   TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);
  194.  
  195. { Double Declining Balance (DDB) }
  196. function DoubleDecliningBalance(Cost, Salvage: Extended;
  197.   Life, Period: Integer): Extended;
  198.  
  199. { Future Value (FVAL) }
  200. function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  201.   Extended; PaymentTime: TPaymentTime): Extended;
  202.  
  203. { Interest Payment (IPAYMT)  }
  204. function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  205.   FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  206.  
  207. { Interest Rate (IRATE) }
  208. function InterestRate(NPeriods: Integer;
  209.   Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  210.  
  211. { Internal Rate of Return. (IRR) Needs array of cash flows. }
  212. function InternalRateOfReturn(Guess: Extended;
  213.   const CashFlows: array of Double): Extended;
  214.  
  215. { Number of Periods (NPER) }
  216. function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  217.   PaymentTime: TPaymentTime): Extended;
  218.  
  219. { Net Present Value. (NPV) Needs array of cash flows. }
  220. function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  221.   PaymentTime: TPaymentTime): Extended;
  222.  
  223. { Payment (PAYMT) }
  224. function Payment(Rate: Extended; NPeriods: Integer;
  225.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  226.  
  227. { Period Payment (PPAYMT) }
  228. function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  229.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  230.  
  231. { Present Value (PVAL) }
  232. function PresentValue(Rate: Extended; NPeriods: Integer;
  233.   Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  234.  
  235. { Straight Line depreciation (SLN) }
  236. function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  237.  
  238. { Sum-of-Years-Digits depreciation (SYD) }
  239. function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  240.  
  241. type
  242.   EInvalidArgument = class(EMathError) end;
  243.  
  244. implementation
  245.  
  246. function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
  247.   var CompoundRN: Extended): Extended; Forward;
  248. //function AnnuityF(R:Extended; N: Integer): Extended; Forward;
  249. function Compound(R: Extended; N: Integer): Extended; Forward;
  250. function RelSmall(X, Y: Extended): Boolean; Forward;
  251.  
  252. type
  253.   TPoly = record
  254.     Neg, Pos, DNeg, DPos: Extended
  255.   end;
  256.  
  257. const
  258.   MaxIterations = 15;
  259.  
  260. procedure ArgError(const Msg: string);
  261. begin
  262.   raise EInvalidArgument.Create(Msg);
  263. end;
  264.  
  265. function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180 }
  266. begin
  267.   Result := Degrees * (PI / 180);
  268. end;
  269.  
  270. function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
  271. begin
  272.   Result := Radians * (180 / PI);
  273. end;
  274.  
  275. function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
  276. begin
  277.   Result := Grads * (PI / 200);
  278. end;
  279.  
  280. function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI}
  281. begin
  282.   Result := Radians * (200 / PI);
  283. end;
  284.  
  285. function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  286. begin
  287.   Result := Cycles * (2 * PI);
  288. end;
  289.  
  290. function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  291. begin
  292.   Result := Radians / (2 * PI);
  293. end;
  294.  
  295. function LnXP1(X: Extended): Extended;
  296. { Return ln(1 + X).  Accurate for X near 0. }
  297. asm
  298.         FLDLN2
  299.         MOV     AX,WORD PTR X+8               { exponent }
  300.         FLD     X
  301.         CMP     AX,$3FFD                      { .4225 }
  302.         JB      @@1
  303.         FLD1
  304.         FADD
  305.         FYL2X
  306.         JMP     @@2
  307. @@1:
  308.         FYL2XP1
  309. @@2:
  310.         FWAIT
  311. end;
  312.  
  313. { Invariant: Y >= 0 & Result*X**Y = X**I.  Init Y = I and Result = 1. }
  314. {function IntPower(X: Extended; I: Integer): Extended;
  315. var
  316.   Y: Integer;
  317. begin
  318.   Y := Abs(I);
  319.   Result := 1.0;
  320.   while Y > 0 do begin
  321.     while not Odd(Y) do
  322.     begin
  323.       Y := Y shr 1;
  324.       X := X * X
  325.     end;
  326.     Dec(Y);
  327.     Result := Result * X
  328.   end;
  329.   if I < 0 then Result := 1.0 / Result
  330. end;
  331. }
  332. function IntPower(Base: Extended; Exponent: Integer): Extended;
  333. asm
  334.         mov     ecx, eax
  335.         cdq
  336.         fld1                      { Result := 1 }
  337.         xor     eax, edx
  338.         sub     eax, edx          { eax := Abs(Exponent) }
  339.         jz      @@3
  340.         fld     Base
  341.         jmp     @@2
  342. @@1:    fmul    ST, ST            { X := Base * Base }
  343. @@2:    shr     eax,1
  344.         jnc     @@1
  345.         fmul    ST(1),ST          { Result := Result * X }
  346.         jnz     @@1
  347.         fstp    st                { pop X from FPU stack }
  348.         cmp     ecx, 0
  349.         jge     @@3
  350.         fld1
  351.         fdivrp                    { Result := 1 / Result }
  352. @@3:
  353.         fwait
  354. end;
  355.  
  356. function Compound(R: Extended; N: Integer): Extended;
  357. { Return (1 + R)**N. }
  358. begin
  359.   Result := IntPower(1.0 + R, N)
  360. end;
  361.  
  362. (*
  363. function AnnuityF(R:Extended; N: Integer): Extended;
  364. {Return (((1+R)**N)-1)/R, the annuity growth factor}
  365. begin
  366.   { 6.1E-5 approx= 2**-14 }
  367.     if ( ABS(R)<6.1E-5 ) then
  368.         AnnuityF := N*(1+(N-1)*R/2) {linear approximation for small R}
  369.     else
  370.       AnnuityF := (Exp(N*LNXP1(R))-1)/R;
  371. end; {AnnuityF}
  372. *)
  373.  
  374. function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
  375.   var CompoundRN: Extended): Extended;
  376. { Set CompoundRN to Compound(R, N),
  377.     return (1+Rate*PaymentTime)*(Compound(R,N)-1)/R;
  378. }
  379. begin
  380.   if R = 0.0 then
  381.   begin
  382.     CompoundRN := 1.0;
  383.     Result := N;
  384.   end
  385.   else
  386.   begin
  387.     { 6.1E-5 approx= 2**-14 }
  388.     if Abs(R) < 6.1E-5 then
  389.     begin
  390.       CompoundRN := Exp(N * LnXP1(R));
  391.       Result := N*(1+(N-1)*R/2);
  392.     end
  393.     else
  394.     begin
  395.       CompoundRN := Compound(R, N);
  396.       Result := (CompoundRN-1) / R
  397.     end;
  398.     if PaymentTime = ptStartOfPeriod then
  399.       Result := Result * (1 + R);
  400.   end;
  401. end; {Annuity2}
  402.  
  403.  
  404. procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);
  405. { Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.
  406.   Accumulate positive and negative terms separately. }
  407. var
  408.   I: Integer;
  409.   Neg, Pos, DNeg, DPos: Extended;
  410. begin
  411.   Neg := 0.0;
  412.   Pos := 0.0;
  413.   DNeg := 0.0;
  414.   DPos := 0.0;
  415.     for I := High(A) downto  Low(A) do
  416.   begin
  417.     DNeg := X * DNeg + Neg;
  418.     Neg := Neg * X;
  419.     DPos := X * DPos + Pos;
  420.     Pos := Pos * X;
  421.     if A[I] >= 0.0 then
  422.       Pos := Pos + A[I]
  423.     else
  424.       Neg := Neg + A[I]
  425.   end;
  426.   Poly.Neg := Neg;
  427.   Poly.Pos := Pos;
  428.   Poly.DNeg := DNeg * X;
  429.   Poly.DPos := DPos * X;
  430. end; {PolyX}
  431.  
  432.  
  433. function RelSmall(X, Y: Extended): Boolean;
  434. { Returns True if X is small relative to Y }
  435. const
  436.   C1: Double = 1E-15;
  437.   C2: Double = 1E-12;
  438. begin
  439.   Result := Abs(X) < (C1 + C2 * Abs(Y))
  440. end;
  441.  
  442. { Math functions. }
  443.  
  444. function ArcCos(X: Extended): Extended;
  445. begin
  446.   Result := ArcTan2(Sqrt(1 - X*X), X);
  447. end;
  448.  
  449. function ArcSin(X: Extended): Extended;
  450. begin
  451.   Result := ArcTan2(X, Sqrt(1 - X*X))
  452. end;
  453.  
  454. function ArcTan2(Y, X: Extended): Extended;
  455. asm
  456.         FLD     Y
  457.         FLD     X
  458.         FPATAN
  459.         FWAIT
  460. end;
  461.  
  462. function Tan(X: Extended): Extended;
  463. {  Tan := Sin(X) / Cos(X) }
  464. asm
  465.         FLD    X
  466.         FPTAN
  467.         FSTP   ST(0)      { FPTAN pushes 1.0 after result }
  468.         FWAIT
  469. end;
  470.  
  471. function CoTan(X: Extended): Extended;
  472. { CoTan := Cos(X) / Sin(X) = 1 / Tan(X) }
  473. asm
  474.         FLD   X
  475.         FPTAN
  476.         FDIVRP
  477.         FWAIT
  478. end;
  479.  
  480. function Hypot(X, Y: Extended): Extended;
  481. { formula: Sqrt(X*X + Y*Y)
  482.   implemented as:  |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision
  483. var
  484.   Temp: Extended;
  485. begin
  486.   X := Abs(X);
  487.   Y := Abs(Y);
  488.   if X > Y then
  489.   begin
  490.     Temp := X;
  491.     X := Y;
  492.     Y := Temp;
  493.   end;
  494.   if X = 0 then
  495.     Result := Y
  496.   else         // Y > X, X <> 0, so Y > 0
  497.     Result := Y * Sqrt(1 + Sqr(X/Y));
  498. end;
  499. }
  500. asm
  501.         FLD     Y
  502.         FABS
  503.         FLD     X
  504.         FABS
  505.         FCOM
  506.         FNSTSW  AX
  507.         TEST    AH,$45
  508.         JNZ      @@1        // if ST > ST(1) then swap
  509.         FXCH    ST(1)      // put larger number in ST(1)
  510. @@1:    FLDZ
  511.         FCOMP
  512.         FNSTSW  AX
  513.         TEST    AH,$40     // if ST = 0, return ST(1)
  514.         JZ      @@2
  515.         FSTP    ST         // eat ST(0)
  516.         JMP     @@3
  517. @@2:    FDIV    ST,ST(1)   // ST := ST / ST(1)
  518.         FMUL    ST,ST      // ST := ST * ST
  519.         FLD1
  520.         FADD               // ST := ST + 1
  521.         FSQRT              // ST := Sqrt(ST)
  522.         FMUL               // ST(1) := ST * ST(1); Pop ST
  523. @@3:    FWAIT
  524. end;
  525.  
  526.  
  527. procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
  528. asm
  529.         FLD     Theta
  530.         FSINCOS
  531.         FSTP    tbyte ptr [edx]    // Cos
  532.         FSTP    tbyte ptr [eax]    // Sin
  533.         FWAIT
  534. end;
  535.  
  536. { Extract exponent and mantissa from X }
  537. procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
  538. { Mantissa ptr in EAX, Exponent ptr in EDX }
  539. asm
  540.         FLD     X
  541.         PUSH    EAX
  542.         MOV     dword ptr [edx], 0    { if X = 0, return 0 }
  543.  
  544.         FTST
  545.         FSTSW   AX
  546.         FWAIT
  547.         SAHF
  548.         JZ      @@Done
  549.  
  550.         FXTRACT                 // ST(1) = exponent, (pushed) ST = fraction
  551.         FXCH
  552.  
  553. // The FXTRACT instruction normalizes the fraction 1 bit higher than
  554. // wanted for the definition of frexp() so we need to tweak the result
  555. // by scaling the fraction down and incrementing the exponent.
  556.  
  557.         FISTP   dword ptr [edx]
  558.         FLD1
  559.         FCHS
  560.         FXCH
  561.         FSCALE                  // scale fraction
  562.         INC     dword ptr [edx] // exponent biased to match
  563.         FSTP ST(1)              // discard -1, leave fraction as TOS
  564.  
  565. @@Done:
  566.         POP     EAX
  567.         FSTP    tbyte ptr [eax]
  568.         FWAIT
  569. end;
  570.  
  571. function Ldexp(X: Extended; P: Integer): Extended;
  572.   { Result := X * (2^P) }
  573. asm
  574.         PUSH    EAX
  575.         FILD    dword ptr [ESP]
  576.         FLD     X
  577.         FSCALE
  578.         POP     EAX
  579.         FSTP    ST(1)
  580.         FWAIT
  581. end;
  582.  
  583. function Ceil(X: Extended): Integer;
  584. begin
  585.   Result := Trunc(X);
  586.   if Frac(X) > 0 then
  587.     Inc(Result);
  588. end;
  589.  
  590. function Floor(X: Extended): Integer;
  591. begin
  592.   Result := Trunc(X);
  593.   if Frac(X) < 0 then
  594.     Dec(Result);
  595. end;
  596.  
  597. { Conversion of bases:  Log.b(X) = Log.a(X) / Log.a(b)  }
  598.  
  599. function Log10(X: Extended): Extended;
  600.   { Log.10(X) := Log.2(X) * Log.10(2) }
  601. asm
  602.         FLDLG2     { Log base ten of 2 }
  603.         FLD     X
  604.         FYL2X
  605.         FWAIT
  606. end;
  607.  
  608. function Log2(X: Extended): Extended;
  609. asm
  610.         FLD1
  611.         FLD     X
  612.         FYL2X
  613.         FWAIT
  614. end;
  615.  
  616. function LogN(Base, X: Extended): Extended;
  617. { Log.N(X) := Log.2(X) / Log.2(N) }
  618. asm
  619.         FLD1
  620.         FLD     X
  621.         FYL2X
  622.         FLD1
  623.         FLD     Base
  624.         FYL2X
  625.         FDIV
  626.         FWAIT
  627. end;
  628.  
  629. function Poly(X: Extended; const Coefficients: array of Double): Extended;
  630. { Horner's method }
  631. var
  632.   I: Integer;
  633. begin
  634.   Result := Coefficients[High(Coefficients)];
  635.   for I := High(Coefficients)-1 downto Low(Coefficients) do
  636.     Result := Result * X + Coefficients[I];
  637. end;
  638.  
  639. function Power(Base, Exponent: Extended): Extended;
  640. begin
  641.   if Exponent = 0.0 then
  642.     Result := 1.0               { n**0 = 1 }
  643.   else if (Base = 0.0) and (Exponent > 0.0) then
  644.     Result := 0.0               { 0**n = 0, n > 0 }
  645.   else if (Frac(Exponent) = 0.0) and (Abs(Exponent) <= MaxInt) then
  646.     Result := IntPower(Base, Trunc(Exponent))
  647.   else
  648.     Result := Exp(Exponent * Ln(Base))
  649. end;
  650.  
  651. { Hyperbolic functions }
  652.  
  653. function CoshSinh(X: Extended; Factor: Double): Extended;
  654. begin
  655.   Result := Exp(X) / 2;
  656.   Result := Result + Factor / Result;
  657. end;
  658.  
  659. function Cosh(X: Extended): Extended;
  660. begin
  661.   Result := CoshSinh(X, 0.25)
  662. end;
  663.  
  664. function Sinh(X: Extended): Extended;
  665. begin
  666.   Result := CoshSinh(X, -0.25)
  667. end;
  668.  
  669. const
  670.   MaxTanhDomain = 5678.22249441322; // Ln(MaxExtended)/2
  671.  
  672. function Tanh(X: Extended): Extended;
  673. begin
  674.   if X > MaxTanhDomain then
  675.     Result := 1.0
  676.   else if X < -MaxTanhDomain then
  677.     Result := -1.0
  678.   else
  679.   begin
  680.     Result := Exp(X);
  681.     Result := Result * Result;
  682.     Result := (Result - 1.0) / (Result + 1.0)
  683.   end;
  684. end;
  685.  
  686. function ArcCosh(X: Extended): Extended;
  687. begin
  688.   if X <= 1.0 then
  689.     Result := 0.0
  690.   else if X > 1.0e10 then
  691.     Result := Ln(2) + Ln(X)
  692.   else
  693.     Result := Ln(X + Sqrt((X - 1.0) * (X + 1.0)));
  694. end;
  695.  
  696. function ArcSinh(X: Extended): Extended;
  697. var
  698.   Neg: Boolean;
  699. begin
  700.   if X = 0 then
  701.     Result := 0
  702.   else
  703.   begin
  704.     Neg := (X < 0);
  705.     X := Abs(X);
  706.     if X > 1.0e10 then
  707.       Result := Ln(2) + Ln(X)
  708.     else
  709.     begin
  710.       Result := X*X;
  711.       Result := LnXP1(X + Result / (1 + Sqrt(1 + Result)));
  712.     end;
  713.     if Neg then Result := -Result;
  714.   end;
  715. end;
  716.  
  717. function ArcTanh(X: Extended): Extended;
  718. var
  719.   Neg: Boolean;
  720. begin
  721.   if X = 0 then
  722.     Result := 0
  723.   else
  724.   begin
  725.     Neg := (X < 0);
  726.     X := Abs(X);
  727.     if X >= 1 then
  728.       Result := MaxExtended
  729.     else
  730.       Result := 0.5 * LnXP1((2.0 * X) / (1.0 - X));
  731.     if Neg then Result := -Result;
  732.   end;
  733. end;
  734.  
  735. { Statistical functions }
  736.  
  737. function Mean(const Data: array of Double): Extended;
  738. begin
  739.   Result := SUM(Data) / (High(Data) - Low(Data) + 1)
  740. end;
  741.  
  742. function MinValue(const Data: array of Double): Double;
  743. var
  744.   I: Integer;
  745. begin
  746.   Result := Data[Low(Data)];
  747.   for I := Low(Data) + 1 to High(Data) do
  748.     if Result > Data[I] then
  749.       Result := Data[I];
  750. end;
  751.  
  752. function MinIntValue(const Data: array of Integer): Integer;
  753. var
  754.   I: Integer;
  755. begin
  756.   Result := Data[Low(Data)];
  757.   for I := Low(Data) + 1 to High(Data) do
  758.     if Result > Data[I] then
  759.       Result := Data[I];
  760. end;
  761.  
  762. function MaxValue(const Data: array of Double): Double;
  763. var
  764.   I: Integer;
  765. begin
  766.   Result := Data[Low(Data)];
  767.   for I := Low(Data) + 1 to High(Data) do
  768.     if Result < Data[I] then
  769.       Result := Data[I];
  770. end;
  771.  
  772. function MaxIntValue(const Data: array of Integer): Integer;
  773. var
  774.   I: Integer;
  775. begin
  776.   Result := Data[Low(Data)];
  777.   for I := Low(Data) + 1 to High(Data) do
  778.     if Result < Data[I] then
  779.       Result := Data[I];
  780. end;
  781.  
  782. procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  783. var
  784.   S: Extended;
  785.   N,I: Integer;
  786. begin
  787.   N := High(Data)- Low(Data) + 1;
  788.   if N = 1 then
  789.   begin
  790.     Mean := Data[0];
  791.     StdDev := Data[0];
  792.     Exit;
  793.   end;
  794.   Mean := Sum(Data) / N;
  795.   S := 0;               // sum differences from the mean, for greater accuracy
  796.   for I := Low(Data) to High(Data) do
  797.     S := S + Sqr(Mean - Data[I]);
  798.   StdDev := Sqrt(S / (N - 1));
  799. end;
  800.  
  801. procedure MomentSkewKurtosis(const Data: array of Double;
  802.   var M1, M2, M3, M4, Skew, Kurtosis: Extended);
  803. var
  804.   Sum, SumSquares, SumCubes, SumQuads, OverN, Accum, M1Sqr, S2N, S3N: Extended;
  805.   I: Integer;
  806. begin
  807.   OverN := 1 / (High(Data) - Low(Data) + 1);
  808.   Sum := 0;
  809.   SumSquares := 0;
  810.   SumCubes := 0;
  811.   SumQuads := 0;
  812.   for I := Low(Data) to High(Data) do
  813.   begin
  814.     Sum := Sum + Data[I];
  815.     Accum := Sqr(Data[I]);
  816.     SumSquares := SumSquares + Accum;
  817.     Accum := Accum*Data[I];
  818.     SumCubes := SumCubes + Accum;
  819.     SumQuads := SumQuads + Accum*Data[I];
  820.   end;
  821.   M1 := Sum * OverN;
  822.   M1Sqr := Sqr(M1);
  823.   S2N := SumSquares * OverN;
  824.   S3N := SumCubes * OverN;
  825.   M2 := S2N - M1Sqr;
  826.   M3 := S3N - (M1 * 3 * S2N) + 2*M1Sqr*M1;
  827.   M4 := (SumQuads * OverN) - (M1 * 4 * S3N) + (M1Sqr*6*S2N - 3*Sqr(M1Sqr));
  828.   Skew := M3 * Power(M2, -3/2);   // = M3 / Power(M2, 3/2)
  829.   Kurtosis := M4 / Sqr(M2);
  830. end;
  831.  
  832. function Norm(const Data: array of Double): Extended;
  833. begin
  834.   Result := Sqrt(SumOfSquares(Data));
  835. end;
  836.  
  837. function PopnStdDev(const Data: array of Double): Extended;
  838. begin
  839.   Result := Sqrt(PopnVariance(Data))
  840. end;
  841.  
  842. function PopnVariance(const Data: array of Double): Extended;
  843. begin
  844.   Result := TotalVariance(Data) / (High(Data) - Low(Data) + 1)
  845. end;
  846.  
  847. function RandG(Mean, StdDev: Extended): Extended;
  848. { Marsaglia-Bray algorithm }
  849. var
  850.   U1, S2: Extended;
  851. begin
  852.   repeat
  853.     U1 := 2*Random - 1;
  854.     S2 := Sqr(U1) + Sqr(2*Random-1);
  855.   until S2 < 1;
  856.   Result := Sqrt(-2*Ln(S2)/S2) * U1 * StdDev + Mean;
  857. end;
  858.  
  859. function StdDev(const Data: array of Double): Extended;
  860. begin
  861.   Result := Sqrt(Variance(Data))
  862. end;
  863.  
  864. procedure RaiseOverflowError; forward;
  865.  
  866. function SumInt(const Data: array of Integer): Integer;
  867. {var
  868.   I: Integer;
  869. begin
  870.   Result := 0;
  871.   for I := Low(Data) to High(Data) do
  872.     Result := Result + Data[I]
  873. end; }
  874. asm  // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1
  875.      // loop unrolled 4 times, 5 clocks per loop, 1.2 clocks per datum
  876.       PUSH EBX
  877.       MOV  ECX, EAX         // ecx = ptr to data
  878.       MOV  EBX, EDX
  879.       XOR  EAX, EAX
  880.       AND  EDX, not 3
  881.       AND  EBX, 3
  882.       SHL  EDX, 2
  883.       JMP  @Vector.Pointer[EBX*4]
  884. @Vector:
  885.       DD @@1
  886.       DD @@2
  887.       DD @@3
  888.       DD @@4
  889. @@4:
  890.       ADD  EAX, [ECX+12+EDX]
  891.       JO   RaiseOverflowError
  892. @@3:
  893.       ADD  EAX, [ECX+8+EDX]
  894.       JO   RaiseOverflowError
  895. @@2:
  896.       ADD  EAX, [ECX+4+EDX]
  897.       JO   RaiseOverflowError
  898. @@1:
  899.       ADD  EAX, [ECX+EDX]
  900.       JO   RaiseOverflowError
  901.       SUB  EDX,16
  902.       JNS  @@4
  903.       POP  EBX
  904. end;
  905.  
  906. procedure RaiseOverflowError;
  907. begin
  908.   raise EIntOverflow.Create(SIntOverflow);
  909. end;
  910.  
  911. function SUM(const Data: array of Double): Extended;
  912. {var
  913.   I: Integer;
  914. begin
  915.   Result := 0.0;
  916.   for I := Low(Data) to High(Data) do
  917.     Result := Result + Data[I]
  918. end; }
  919. asm  // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1
  920.      // Uses 4 accumulators to minimize read-after-write delays and loop overhead
  921.      // 5 clocks per loop, 4 items per loop = 1.2 clocks per item
  922.        FLDZ
  923.        MOV      ECX, EDX
  924.        FLD      ST(0)
  925.        AND      EDX, not 3
  926.        FLD      ST(0)
  927.        AND      ECX, 3
  928.        FLD      ST(0)
  929.        SHL      EDX, 3      // count * sizeof(Double) = count * 8
  930.        JMP      @Vector.Pointer[ECX*4]
  931. @Vector:
  932.        DD @@1
  933.        DD @@2
  934.        DD @@3
  935.        DD @@4
  936. @@4:   FADD     qword ptr [EAX+EDX+24]    // 1
  937.        FXCH     ST(3)                     // 0
  938. @@3:   FADD     qword ptr [EAX+EDX+16]    // 1
  939.        FXCH     ST(2)                     // 0
  940. @@2:   FADD     qword ptr [EAX+EDX+8]     // 1
  941.        FXCH     ST(1)                     // 0
  942. @@1:   FADD     qword ptr [EAX+EDX]       // 1
  943.        FXCH     ST(2)                     // 0
  944.        SUB      EDX, 32
  945.        JNS      @@4
  946.        FADDP    ST(3),ST                  // ST(3) := ST + ST(3); Pop ST
  947.        FADD                               // ST(1) := ST + ST(1); Pop ST
  948.        FADD                               // ST(1) := ST + ST(1); Pop ST
  949.        FWAIT
  950. end;
  951.  
  952. function SumOfSquares(const Data: array of Double): Extended;
  953. var
  954.   I: Integer;
  955. begin
  956.   Result := 0.0;
  957.   for I := Low(Data) to High(Data) do
  958.     Result := Result + Sqr(Data[I]);
  959. end;
  960.  
  961. procedure SumsAndSquares(const Data: array of Double; var Sum, SumOfSquares: Extended);
  962. {var
  963.   I: Integer;
  964. begin
  965.   Sum := 0;
  966.   SumOfSquares := 0;
  967.   for I := Low(Data) to High(Data) do
  968.   begin
  969.     Sum := Sum + Data[I];
  970.     SumOfSquares := SumOfSquares + Data[I]*Data[I];
  971.   end;
  972. end;  }
  973. asm  // IN:  EAX = ptr to Data
  974.      //      EDX = High(Data) = Count - 1
  975.      //      ECX = ptr to Sum
  976.      // Est. 17 clocks per loop, 4 items per loop = 4.5 clocks per data item
  977.        FLDZ                 // init Sum accumulator
  978.        PUSH     ECX
  979.        MOV      ECX, EDX
  980.        FLD      ST(0)       // init Sqr1 accum.
  981.        AND      EDX, not 3
  982.        FLD      ST(0)       // init Sqr2 accum.
  983.        AND      ECX, 3
  984.        FLD      ST(0)       // init/simulate last data item left in ST
  985.        SHL      EDX, 3      // count * sizeof(Double) = count * 8
  986.        JMP      @Vector.Pointer[ECX*4]
  987. @Vector:
  988.        DD @@1
  989.        DD @@2
  990.        DD @@3
  991.        DD @@4
  992. @@4:   FADD                            // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
  993.        FLD     qword ptr [EAX+EDX+24]  // Load Data1
  994.        FADD    ST(3),ST                // Sum := Sum + Data1
  995.        FMUL    ST,ST                   // Data1 := Sqr(Data1)
  996. @@3:   FLD     qword ptr [EAX+EDX+16]  // Load Data2
  997.        FADD    ST(4),ST                // Sum := Sum + Data2
  998.        FMUL    ST,ST                   // Data2 := Sqr(Data2)
  999.        FXCH                            // Move Sqr(Data1) into ST(0)
  1000.        FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data1); Pop Data1
  1001. @@2:   FLD     qword ptr [EAX+EDX+8]   // Load Data3
  1002.        FADD    ST(4),ST                // Sum := Sum + Data3
  1003.        FMUL    ST,ST                   // Data3 := Sqr(Data3)
  1004.        FXCH                            // Move Sqr(Data2) into ST(0)
  1005.        FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data2); Pop Data2
  1006. @@1:   FLD     qword ptr [EAX+EDX]     // Load Data4
  1007.        FADD    ST(4),ST                // Sum := Sum + Data4
  1008.        FMUL    ST,ST                   // Sqr(Data4)
  1009.        FXCH                            // Move Sqr(Data3) into ST(0)
  1010.        FADDP   ST(3),ST                // Sqr1 := Sqr1 + Sqr(Data3); Pop Data3
  1011.        SUB     EDX,32
  1012.        JNS     @@4
  1013.        FADD                         // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
  1014.        POP     ECX
  1015.        FADD                         // Sqr1 := Sqr2 + Sqr1; Pop Sqr2
  1016.        FXCH                         // Move Sum1 into ST(0)
  1017.        MOV     EAX, SumOfSquares
  1018.        FSTP    tbyte ptr [ECX]      // Sum := Sum1; Pop Sum1
  1019.        FSTP    tbyte ptr [EAX]      // SumOfSquares := Sum1; Pop Sum1
  1020.        FWAIT
  1021. end;
  1022.  
  1023. function TotalVariance(const Data: array of Double): Extended;
  1024. var
  1025.   Sum, SumSquares: Extended;
  1026. begin
  1027.   SumsAndSquares(Data, Sum, SumSquares);
  1028.   Result := SumSquares - Sqr(Sum)/(High(Data) - Low(Data) + 1);
  1029. end;
  1030.  
  1031. function Variance(const Data: array of Double): Extended;
  1032. begin
  1033.   Result := TotalVariance(Data) / (High(Data) - Low(Data))
  1034. end;
  1035.  
  1036.  
  1037. { Depreciation functions. }
  1038.  
  1039. function DoubleDecliningBalance(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  1040. { dv := cost * (1 - 2/life)**(period - 1)
  1041.   DDB = (2/life) * dv
  1042.   if DDB > dv - salvage then DDB := dv - salvage
  1043.   if DDB < 0 then DDB := 0
  1044. }
  1045. var
  1046.   DepreciatedVal, Factor: Extended;
  1047. begin
  1048.   Result := 0;
  1049.     if (Period < 1) or (Life < Period) or (Life < 1) or (Cost <= Salvage) then
  1050.     Exit;
  1051.  
  1052.     {depreciate everything in period 1 if life is only one or two periods}
  1053.     if ( Life <= 2 ) then
  1054.   begin
  1055.         if ( Period = 1 ) then
  1056.       DoubleDecliningBalance:=Cost-Salvage
  1057.         else
  1058.             DoubleDecliningBalance:=0; {all depreciation occurred in first period}
  1059.         exit;
  1060.   end;
  1061.   Factor := 2.0 / Life;
  1062.  
  1063.   DepreciatedVal := Cost * IntPower((1.0 - Factor), Period - 1);
  1064.        {DepreciatedVal is Cost-(sum of previous depreciation results)}
  1065.  
  1066.   Result := Factor * DepreciatedVal;
  1067.        {Nominal computed depreciation for this period.  The rest of the
  1068.             function applies limits to this nominal value. }
  1069.  
  1070.     {Only depreciate until total depreciation equals cost-salvage.}
  1071.   if Result > DepreciatedVal - Salvage then
  1072.     Result := DepreciatedVal - Salvage;
  1073.  
  1074.     {No more depreciation after salvage value is reached.  This is mostly a nit.
  1075.      If Result is negative at this point, it's very close to zero.}
  1076.   if Result < 0.0 then
  1077.     Result := 0.0;
  1078. end;
  1079.  
  1080. function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  1081. { Spreads depreciation linearly over life. }
  1082. begin
  1083.   if Life < 1 then ArgError('SLNDepreciation');
  1084.   Result := (Cost - Salvage) / Life
  1085. end;
  1086.  
  1087. function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  1088. { SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }
  1089. { Note: life*(life+1)/2 = 1+2+3+...+life "sum of years"
  1090.                 The depreciation factor varies from life/sum_of_years in first period = 1
  1091.                                                                              downto  1/sum_of_years in last period = life.
  1092.                 Total depreciation over life is cost-salvage.}
  1093. var
  1094.   X1, X2: Extended;
  1095. begin
  1096.   Result := 0;
  1097.     if (Period < 1) or (Life < Period) or (Cost <= Salvage) then Exit;
  1098.   X1 := 2 * (Life - Period + 1);
  1099.   X2 := Life * (Life + 1);
  1100.   Result := (Cost - Salvage) * X1 / X2
  1101. end;
  1102.  
  1103. { Discounted cash flow functions. }
  1104.  
  1105. function InternalRateOfReturn(Guess: Extended; const CashFlows: array of Double): Extended;
  1106. {
  1107. Use Newton's method to solve NPV = 0, where NPV is a polynomial in
  1108. x = 1/(1+rate).  Split the coefficients into negative and postive sets:
  1109.   neg + pos = 0, so pos = -neg, so  -neg/pos = 1
  1110. Then solve:
  1111.   log(-neg/pos) = 0
  1112.  
  1113.   Let  t = log(1/(1+r) = -LnXP1(r)
  1114.   then r = exp(-t) - 1
  1115. Iterate on t, then use the last equation to compute r.
  1116. }
  1117. var
  1118.   T, Y: Extended;
  1119.   Poly: TPoly;
  1120.   K, Count: Integer;
  1121.  
  1122.   function ConditionP(const CashFlows: array of Double): Integer;
  1123.   { Guarantees existence and uniqueness of root.  The sign of payments
  1124.     must change exactly once, the net payout must be always > 0 for
  1125.     first portion, then each payment must be >= 0.
  1126.     Returns: 0 if condition not satisfied, > 0 if condition satisfied
  1127.     and this is the index of the first value considered a payback. }
  1128.   var
  1129.     X: Double;
  1130.     I, K: Integer;
  1131.   begin
  1132.     K := High(CashFlows);
  1133.     while (K >= 0) and (CashFlows[K] >= 0.0) do Dec(K);
  1134.     Inc(K);
  1135.     if K > 0 then
  1136.     begin
  1137.       X := 0.0;
  1138.       I := 0;
  1139.       while I < K do begin
  1140.           X := X + CashFlows[I];
  1141.         if X >= 0.0 then
  1142.         begin
  1143.              K := 0;
  1144.             Break
  1145.           end;
  1146.           Inc(I)
  1147.       end
  1148.     end;
  1149.     ConditionP := K
  1150.   end;
  1151.  
  1152. begin
  1153.   InternalRateOfReturn := 0;
  1154.   K := ConditionP(CashFlows);
  1155.   if K < 0 then ArgError('InternalRateOfReturn');
  1156.   if K = 0 then
  1157.   begin
  1158.     if Guess <= -1.0 then ArgError('InternalRateOfReturn');
  1159.     T := -LnXP1(Guess)
  1160.   end else
  1161.     T := 0.0;
  1162.   for Count := 1 to MaxIterations do
  1163.   begin
  1164.     PolyX(CashFlows, Exp(T), Poly);
  1165.     if Poly.Pos <= Poly.Neg then ArgError('InternalRateOfReturn');
  1166.     if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then
  1167.     begin
  1168.       InternalRateOfReturn := -1.0;
  1169.       Exit;
  1170.     end;
  1171.     with Poly do
  1172.       Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);
  1173.     T := T - Y;
  1174.     if RelSmall(Y, T) then
  1175.     begin
  1176.       InternalRateOfReturn := Exp(-T) - 1.0;
  1177.       Exit;
  1178.     end
  1179.   end;
  1180.   ArgError('InternalRateOfReturn');
  1181. end;
  1182.  
  1183. function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  1184.   PaymentTime: TPaymentTime): Extended;
  1185. { Caution: The sign of NPV is reversed from what would be expected for standard
  1186.    cash flows!}
  1187. var
  1188.   rr: Extended;
  1189.   I: Integer;
  1190. begin
  1191.   if Rate <= -1.0 then ArgError('NetPresentValue');
  1192.   rr := 1/(1+Rate);
  1193.   result := 0;
  1194.   for I := High(CashFlows) downto Low(CashFlows) do
  1195.     result := rr * result + CashFlows[I];
  1196.   if PaymentTime = ptEndOfPeriod then result := rr * result;
  1197. end;
  1198.  
  1199. { Annuity functions. }
  1200.  
  1201. {---------------
  1202. From the point of view of A, amounts received by A are positive and
  1203. amounts disbursed by A are negative (e.g. a borrower's loan repayments
  1204. are regarded by the borrower as negative).
  1205.  
  1206. Given interest rate r, number of periods n:
  1207.   compound(r, n) = (1 + r)**n               "Compounding growth factor"
  1208.   annuity(r, n) = (compound(r, n)-1) / r   "Annuity growth factor"
  1209.  
  1210. Given future value fv, periodic payment pmt, present value pv and type
  1211. of payment (start, 1 , or end of period, 0) pmtTime, financial variables satisfy:
  1212.  
  1213.   fv = -pmt*(1 + r*pmtTime)*annuity(r, n) - pv*compound(r, n)
  1214.  
  1215. For fv, pv, pmt:
  1216.  
  1217.   C := compound(r, n)
  1218.   A := (1 + r*pmtTime)*annuity(r, n)
  1219.   Compute both at once in Annuity2.
  1220.  
  1221.   if C > 1E16 then A = C/r, so:
  1222.     fv := meaningless
  1223.     pv := -pmt*(pmtTime+1/r)
  1224.     pmt := -pv*r/(1 + r*pmtTime)
  1225.   else
  1226.     fv := -pmt(1+r*pmtTime)*A - pv*C
  1227.     pv := (-pmt(1+r*pmtTime)*A - fv)/C
  1228.     pmt := (-pv*C-fv)/((1+r*pmtTime)*A)
  1229. ---------------}
  1230.  
  1231. function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,
  1232.   FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):
  1233.   Extended;
  1234. var
  1235.         Crn:extended; { =Compound(Rate,NPeriods) }
  1236.         Crp:extended; { =Compound(Rate,Period-1) }
  1237.         Arn:extended; { =AnnuityF(Rate,NPeriods) }
  1238.  
  1239. begin
  1240.   if Rate <= -1.0 then ArgError('PaymentParts');
  1241.   Crp:=Compound(Rate,Period-1);
  1242.   Arn:=Annuity2(Rate,NPeriods,PaymentTime,Crn);
  1243.   IntPmt:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;
  1244.   PaymentParts:=(-FutureValue-PresentValue)*Crp/Arn;
  1245. end;
  1246.  
  1247. function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  1248.   Extended; PaymentTime: TPaymentTime): Extended;
  1249. var
  1250.   Annuity, CompoundRN: Extended;
  1251. begin
  1252.   if Rate <= -1.0 then ArgError('FutureValue');
  1253.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1254.   if CompoundRN > 1.0E16 then ArgError('FutureValue');
  1255.   FutureValue := -Payment * Annuity - PresentValue * CompoundRN
  1256. end;
  1257.  
  1258. function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  1259.   FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1260. var
  1261.   Crp:extended; { compound(rate,period-1)}
  1262.   Crn:extended; { compound(rate,nperiods)}
  1263.   Arn:extended; { annuityf(rate,nperiods)}
  1264. begin
  1265.   if (Rate <= -1.0)
  1266.     or (Period < 1) or (Period > NPeriods) then ArgError('InterestPayment');
  1267.     Crp:=Compound(Rate,Period-1);
  1268.     Arn:=Annuity2(Rate,Nperiods,PaymentTime,Crn);
  1269.     InterestPayment:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;
  1270. end;
  1271.  
  1272. function InterestRate(NPeriods: Integer;
  1273.   Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1274. {
  1275. Given:
  1276.   First and last payments are non-zero and of opposite signs.
  1277.   Number of periods N >= 2.
  1278. Convert data into cash flow of first, N-1 payments, last with
  1279. first < 0, payment > 0, last > 0.
  1280. Compute the IRR of this cash flow:
  1281.   0 = first + pmt*x + pmt*x**2 + ... + pmt*x**(N-1) + last*x**N
  1282. where x = 1/(1 + rate).
  1283. Substitute x = exp(t) and apply Newton's method to
  1284.   f(t) = log(pmt*x + ... + last*x**N) / -first
  1285. which has a unique root given the above hypotheses.
  1286. }
  1287. var
  1288.   X, Y, Z, First, Pmt, Last, T, ET, EnT, ET1: Extended;
  1289.   Count: Integer;
  1290.   Reverse: Boolean;
  1291.  
  1292.   function LostPrecision(X: Extended): Boolean;
  1293.   asm
  1294.         XOR     EAX, EAX
  1295.         MOV     BX,WORD PTR X+8
  1296.         INC     EAX
  1297.         AND     EBX, $7FF0
  1298.         JZ      @@1
  1299.         CMP     EBX, $7FF0
  1300.         JE      @@1
  1301.         XOR     EAX,EAX
  1302.   @@1:
  1303.   end;
  1304.  
  1305. begin
  1306.   Result := 0;
  1307.   if NPeriods <= 0 then ArgError('InterestRate');
  1308.   Pmt := Payment;
  1309.   if PaymentTime = ptEndOfPeriod then
  1310.   begin
  1311.     X := PresentValue;
  1312.     Y := FutureValue + Payment
  1313.   end
  1314.   else
  1315.   begin
  1316.     X := PresentValue + Payment;
  1317.     Y := FutureValue
  1318.   end;
  1319.   First := X;
  1320.   Last := Y;
  1321.   Reverse := False;
  1322.   if First * Payment > 0.0 then
  1323.   begin
  1324.     Reverse := True;
  1325.     T := First;
  1326.     First := Last;
  1327.     Last := T
  1328.   end;
  1329.   if first > 0.0 then
  1330.   begin
  1331.     First := -First;
  1332.     Pmt := -Pmt;
  1333.     Last := -Last
  1334.   end;
  1335.   if (First = 0.0) or (Last < 0.0) then ArgError('InterestRate');
  1336.   T := 0.0;                     { Guess at solution }
  1337.   for Count := 1 to MaxIterations do
  1338.   begin
  1339.     EnT := Exp(NPeriods * T);
  1340.     if {LostPrecision(EnT)} ent=(ent+1) then
  1341.     begin
  1342.       Result := -Pmt / First;
  1343.       if Reverse then
  1344.         Result := Exp(-LnXP1(Result)) - 1.0;
  1345.       Exit;
  1346.     end;
  1347.     ET := Exp(T);
  1348.     ET1 := ET - 1.0;
  1349.     if ET1 = 0.0 then
  1350.     begin
  1351.       X := NPeriods;
  1352.       Y := X * (X - 1.0) / 2.0
  1353.     end
  1354.     else
  1355.     begin
  1356.       X := ET * (Exp((NPeriods - 1) * T)-1.0) / ET1;
  1357.       Y := (NPeriods * EnT - ET - X * ET) / ET1
  1358.     end;
  1359.     Z := Pmt * X + Last * EnT;
  1360.     Y := Ln(Z / -First) / ((Pmt * Y + Last * NPeriods *EnT) / Z);
  1361.     T := T - Y;
  1362.     if RelSmall(Y, T) then
  1363.     begin
  1364.       if not Reverse then T := -T;
  1365.       InterestRate := Exp(T)-1.0;
  1366.       Exit;
  1367.     end
  1368.   end;
  1369.   ArgError('InterestRate');
  1370. end;
  1371.  
  1372. function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  1373.   PaymentTime: TPaymentTime): Extended;
  1374.  
  1375. { If Rate = 0 then nper := -(pv + fv) / pmt
  1376.   else cf := pv + pmt * (1 + rate*pmtTime) / rate
  1377.        nper := LnXP1(-(pv + fv) / cf) / LnXP1(rate) }
  1378.  
  1379. var
  1380.   PVRPP: Extended; { =PV*Rate+Payment } {"initial cash flow"}
  1381.   T:     Extended;
  1382.  
  1383. begin
  1384.  
  1385.   if Rate <= -1.0 then ArgError('NumberOfPeriods');
  1386.  
  1387. {whenever both Payment and PaymentTime are given together, the PaymentTime has the effect
  1388.  of modifying the effective Payment by the interest accrued on the Payment}
  1389.  
  1390.   if ( PaymentTime=ptStartOfPeriod ) then
  1391.     Payment:=Payment*(1+Rate);
  1392.  
  1393. {if the payment exactly matches the interest accrued periodically on the
  1394.  presentvalue, then an infinite number of payments are going to be
  1395.  required to effect a change from presentvalue to futurevalue.  The
  1396.  following catches that specific error where payment is exactly equal,
  1397.  but opposite in sign to the interest on the present value.  If PVRPP
  1398.  ("initial cash flow") is simply close to zero, the computation will
  1399.  be numerically unstable, but not as likely to cause an error.}
  1400.  
  1401.   PVRPP:=PresentValue*Rate+Payment;
  1402.   if PVRPP=0 then ArgError('NumberOfPeriods');
  1403.  
  1404.   { 6.1E-5 approx= 2**-14 }
  1405.   if ( ABS(Rate)<6.1E-5 ) then
  1406.     Result:=-(PresentValue+FutureValue)/PVRPP
  1407.   else
  1408.   begin
  1409.  
  1410. {starting with the initial cash flow, each compounding period cash flow
  1411.  should result in the current value approaching the final value.  The
  1412.  following test combines a number of simultaneous conditions to ensure
  1413.  reasonableness of the cashflow before computing the NPER.}
  1414.  
  1415.     T:= -(PresentValue+FutureValue)*Rate/PVRPP;
  1416.     if  T<=-1.0  then ArgError('NumberOfPeriods');
  1417.     Result := LnXP1(T) / LnXP1(Rate)
  1418.   end;
  1419.   NumberOfPeriods:=Result;
  1420. end;
  1421.  
  1422. function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:
  1423.   Extended; PaymentTime: TPaymentTime): Extended;
  1424. var
  1425.   Annuity, CompoundRN: Extended;
  1426. begin
  1427.   if Rate <= -1.0 then ArgError('Payment');
  1428.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1429.   if CompoundRN > 1.0E16 then
  1430.     Payment := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)
  1431.   else
  1432.     Payment := (-PresentValue * CompoundRN - FutureValue) / Annuity
  1433. end;
  1434.  
  1435. function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  1436.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1437. var
  1438.   Junk: Extended;
  1439. begin
  1440.   if (Rate <= -1.0) or (Period < 1) or (Period > NPeriods) then ArgError('PeriodPayment');
  1441.   PeriodPayment := PaymentParts(Period, NPeriods, Rate, PresentValue,
  1442.        FutureValue, PaymentTime, Junk);
  1443. end;
  1444.  
  1445. function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:
  1446.   Extended; PaymentTime: TPaymentTime): Extended;
  1447. var
  1448.   Annuity, CompoundRN: Extended;
  1449. begin
  1450.   if Rate <= -1.0 then ArgError('PresentValue');
  1451.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1452.   if CompoundRN > 1.0E16 then
  1453.     PresentValue := -(Payment / Rate * Integer(PaymentTime) * Payment)
  1454.   else
  1455.     PresentValue := (-Payment * Annuity - FutureValue) / CompoundRN
  1456. end;
  1457.  
  1458. end.
  1459.