home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 1999 February / DPPCPRO0299.ISO / February / Delphi / Install / DATA.Z / MATH.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1996-06-09  |  37.3 KB  |  1,311 lines
  1.  
  2. {*******************************************************}
  3. {                                                       }
  4. {       Delphi Runtime Library                          }
  5. {       Math Unit                                       }
  6. {                                                       }
  7. {       Copyright (C) 1996 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 SumOfSquares(const Data: array of Double): Extended;
  132. procedure SumsAndSquares(const Data: array of Double;
  133.   var Sum, SumOfSquares: Extended) register;
  134.  
  135. { MinValue: Returns the smallest signed value in the data array (MIN) }
  136. function MinValue(const Data: array of Double): Double;
  137.  
  138. { MaxValue: Returns the largest signed value in the data array (MAX) }
  139. function MaxValue(const Data: array of Double): Double;
  140.  
  141. { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
  142. function StdDev(const Data: array of Double): Extended;
  143.  
  144. { MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
  145.   calculating them separately.  Less accurate when the mean is very large
  146.   (> 10e7) or the variance is very small. }
  147. procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  148.  
  149. { Population Standard Deviation (STDP): Sqrt(PopnVariance).
  150.   Used in some business and financial calculations. }
  151. function PopnStdDev(const Data: array of Double): Extended;
  152.  
  153. { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
  154. function Variance(const Data: array of Double): Extended;
  155.  
  156. { Population Variance (VAR or VARP): TotalVariance/ N }
  157. function PopnVariance(const Data: array of Double): Extended;
  158.  
  159. { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
  160. function TotalVariance(const Data: array of Double): Extended;
  161.  
  162. { Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
  163. function Norm(const Data: array of Double): Extended;
  164.  
  165. { MomentSkewKurtosis: Calculates the core factors of statistical analysis:
  166.   the first four moments plus the coefficients of skewness and kurtosis.
  167.   M1 is the Mean.  M2 is the Variance.
  168.   Skew reflects symmetry of distribution: M3 / (M2**(3/2))
  169.   Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
  170. procedure MomentSkewKurtosis(const Data: array of Double;
  171.   var M1, M2, M3, M4, Skew, Kurtosis: Extended);
  172.  
  173. { RandG produces random numbers with Gaussian distribution about the mean.
  174.   Useful for simulating data with sampling errors. }
  175. function RandG(Mean, StdDev: Extended): Extended;
  176.  
  177. {-----------------------------------------------------------------------
  178. Financial functions.  Standard set from Quattro Pro.
  179.  
  180. Parameter conventions:
  181.  
  182. From the point of view of A, amounts received by A are positive and
  183. amounts disbursed by A are negative (e.g. a borrower's loan repayments
  184. are regarded by the borrower as negative).
  185.  
  186. Interest rates are per payment period.  11% annual percentage rate on a
  187. loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
  188.  
  189. -----------------------------------------------------------------------}
  190.  
  191. type
  192.   TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);
  193.  
  194. { Double Declining Balance (DDB) }
  195. function DoubleDecliningBalance(Cost, Salvage: Extended;
  196.   Life, Period: Integer): Extended;
  197.  
  198. { Future Value (FVAL) }
  199. function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  200.   Extended; PaymentTime: TPaymentTime): Extended;
  201.  
  202. { Interest Payment (IPAYMT)  }
  203. function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  204.   FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  205.  
  206. { Interest Rate (IRATE) }
  207. function InterestRate(NPeriods: Integer;
  208.   Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  209.  
  210. { Internal Rate of Return. (IRR) Needs array of cash flows. }
  211. function InternalRateOfReturn(Guess: Extended;
  212.   const CashFlows: array of Double): Extended;
  213.  
  214. { Number of Periods (NPER) }
  215. function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  216.   PaymentTime: TPaymentTime): Extended;
  217.  
  218. { Net Present Value. (NPV) Needs array of cash flows. }
  219. function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  220.   PaymentTime: TPaymentTime): Extended;
  221.  
  222. { Payment (PAYMT) }
  223. function Payment(Rate: Extended; NPeriods: Integer;
  224.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  225.  
  226. { Period Payment (PPAYMT) }
  227. function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  228.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  229.  
  230. { Present Value (PVAL) }
  231. function PresentValue(Rate: Extended; NPeriods: Integer;
  232.   Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  233.  
  234. { Straight Line depreciation (SLN) }
  235. function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  236.  
  237. { Sum-of-Years-Digits depreciation (SYD) }
  238. function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  239.  
  240. type
  241.   EInvalidArgument = class(EMathError) end;
  242.  
  243. implementation
  244.  
  245. function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
  246.   var CompoundRN: Extended): Extended; Forward;
  247. function Compound(R: Extended; N: Integer): Extended; Forward;
  248. function RelSmall(X, Y: Extended): Boolean; Forward;
  249.  
  250. type
  251.   TPoly = record
  252.     Neg, Pos, DNeg, DPos: Extended
  253.   end;
  254.  
  255. const
  256.   MaxIterations = 15;
  257.  
  258. procedure ArgError;
  259. begin
  260.   raise EInvalidArgument.Create('');
  261. end;
  262.  
  263. function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180 }
  264. begin
  265.   Result := Degrees * (PI / 180);
  266. end;
  267.  
  268. function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
  269. begin
  270.   Result := Radians * (180 / PI);
  271. end;
  272.  
  273. function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
  274. begin
  275.   Result := Grads * (PI / 200);
  276. end;
  277.  
  278. function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI}
  279. begin
  280.   Result := Radians * (200 / PI);
  281. end;
  282.  
  283. function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
  284. begin
  285.   Result := Cycles * (2 * PI);
  286. end;
  287.  
  288. function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
  289. begin
  290.   Result := Radians / (2 * PI);
  291. end;
  292.  
  293. function LnXP1(X: Extended): Extended;
  294. { Return ln(1 + X).  Accurate for X near 0. }
  295. asm
  296.         FLDLN2
  297.         MOV     AX,WORD PTR X+8               { exponent }
  298.         FLD     X
  299.         CMP     AX,$3FFD                      { .4225 }
  300.         JB      @@1
  301.         FLD1
  302.         FADD
  303.         FYL2X
  304.         JMP     @@2
  305. @@1:
  306.         FYL2XP1
  307. @@2:
  308.         FWAIT
  309. end;
  310.  
  311. { Invariant: Y >= 0 & Result*X**Y = X**I.  Init Y = I and Result = 1. }
  312. {function IntPower(X: Extended; I: Integer): Extended;
  313. var
  314.   Y: Integer;
  315. begin
  316.   Y := Abs(I);
  317.   Result := 1.0;
  318.   while Y > 0 do begin
  319.     while not Odd(Y) do
  320.     begin
  321.       Y := Y shr 1;
  322.       X := X * X
  323.     end;
  324.     Dec(Y);
  325.     Result := Result * X
  326.   end;
  327.   if I < 0 then Result := 1.0 / Result
  328. end;
  329. }
  330. function IntPower(Base: Extended; Exponent: Integer): Extended;
  331. asm
  332.         mov     ecx, eax
  333.         cdq
  334.         fld1                      { Result := 1 }
  335.         xor     eax, edx
  336.         sub     eax, edx          { eax := Abs(Exponent) }
  337.         jz      @@3
  338.         fld     Base
  339.         jmp     @@2
  340. @@1:    fmul    ST, ST            { X := Base * Base }
  341. @@2:    shr     eax,1
  342.         jnc     @@1
  343.         fmul    ST(1),ST          { Result := Result * X }
  344.         jnz     @@1
  345.         fstp    st                { pop X from FPU stack }
  346.         cmp     ecx, 0
  347.         jge     @@3
  348.         fld1
  349.         fdivrp                    { Result := 1 / Result }
  350. @@3:
  351.         fwait
  352. end;
  353.  
  354. function Compound(R: Extended; N: Integer): Extended;
  355. { Return (1 + R)**N. }
  356. begin
  357.   Result := IntPower(1.0 + R, N)
  358. end;
  359.  
  360. function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
  361.   var CompoundRN: Extended): Extended;
  362. { Set CompoundRN to Compound(R, N),
  363.   return (1 + Rate*PaymentTime) * (1 - Compound(R, N)) / R }
  364. begin
  365.   if R = 0.0 then
  366.   begin
  367.     CompoundRN := 1.0;
  368.     Result := -N
  369.   end
  370.   else
  371.   begin
  372.     { 6.1E-5 approx= 2**-14 }
  373.     if Abs(R) < 6.1E-5 then
  374.     begin
  375.       CompoundRN := Exp(N * LnXP1(R));
  376.       Result := (CompoundRN - 1.0) / -R;
  377.     end
  378.     else
  379.     begin
  380.       CompoundRN := Compound(R, N);
  381.       Result := (1 - CompoundRN) / R
  382.     end;
  383. //    if PaymentTime = ptEndOfPeriod then Result := Result + R
  384.     if PaymentTime = ptStartOfPeriod then Result := Result * (1 + R);
  385.   end
  386. end;
  387.  
  388. procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);
  389. { Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.
  390.   Accumulate positive and negative terms separately. }
  391. var
  392.   I: Integer;
  393.   Neg, Pos, DNeg, DPos: Extended;
  394. begin
  395.   Neg := 0.0;
  396.   Pos := 0.0;
  397.   DNeg := 0.0;
  398.   DPos := 0.0;
  399.   for I := Low(A) to High(A) do
  400.   begin
  401.     DNeg := X * DNeg + Neg;
  402.     Neg := Neg * X;
  403.     DPos := X * DPos + Pos;
  404.     Pos := Pos * X;
  405.     if A[I] >= 0.0 then Pos := Pos + A[I] else Neg := Neg + A[I]
  406.   end;
  407.   Poly.Neg := Neg;
  408.   Poly.Pos := Pos;
  409.   Poly.DNeg := DNeg * X;
  410.   Poly.DPos := DPos * X;
  411. end;
  412.  
  413. function RelSmall(X, Y: Extended): Boolean;
  414. { Returns True if X is small relative to Y }
  415. const
  416.   C1: Double = 1E-15;
  417.   C2: Double = 1E-12;
  418. begin
  419.   Result := Abs(X) < (C1 + C2 * Abs(Y))
  420. end;
  421.  
  422. { Math functions. }
  423.  
  424. function ArcCos(X: Extended): Extended;
  425. begin
  426.   Result := ArcTan2(Sqrt(1 - X*X), X);
  427. end;
  428.  
  429. function ArcSin(X: Extended): Extended;
  430. begin
  431.   Result := ArcTan2(X, Sqrt(1 - X*X))
  432. end;
  433.  
  434. function ArcTan2(Y, X: Extended): Extended;
  435. asm
  436.         FLD     Y
  437.         FLD     X
  438.         FPATAN
  439.         FWAIT
  440. end;
  441.  
  442. function Tan(X: Extended): Extended;
  443. {  Tan := Sin(X) / Cos(X) }
  444. asm
  445.         FLD    X
  446.         FPTAN
  447.         FSTP   ST(0)      { FPTAN pushes 1.0 after result }
  448.         FWAIT
  449. end;
  450.  
  451. function CoTan(X: Extended): Extended;
  452. { CoTan := Cos(X) / Sin(X) = 1 / Tan(X) }
  453. asm
  454.         FLD   X
  455.         FPTAN
  456.         FDIVRP
  457.         FWAIT
  458. end;
  459.  
  460. function Hypot(X, Y: Extended): Extended;
  461. { formula: Sqrt(X*X + Y*Y)
  462.   implemented as:  |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision
  463. var
  464.   Temp: Extended;
  465. begin
  466.   X := Abs(X);
  467.   Y := Abs(Y);
  468.   if X > Y then
  469.   begin
  470.     Temp := X;
  471.     X := Y;
  472.     Y := Temp;
  473.   end;
  474.   if X = 0 then
  475.     Result := Y
  476.   else         // Y > X, X <> 0, so Y > 0
  477.     Result := Y * Sqrt(1 + Sqr(X/Y));
  478. end;
  479. }
  480. asm
  481.         FLD     Y
  482.         FABS
  483.         FLD     X
  484.         FABS
  485.         FCOM
  486.         FNSTSW  AX
  487.         TEST    AH,$F
  488.         JNZ      @@1        // if ST > ST(1) then swap
  489.         FXCH    ST(1)      // put larger number in ST(1)
  490. @@1:    FLDZ
  491.         FCOMP
  492.         FNSTSW  AX
  493.         TEST    AH,$4      // if ST = 0, return ST(1)
  494.         JZ      @@2
  495.         FSTP    ST         // eat ST(0)
  496.         JMP     @@3
  497. @@2:    FDIV    ST,ST(1)   // ST := ST / ST(1)
  498.         FMUL    ST,ST      // ST := ST * ST
  499.         FLD1
  500.         FADD               // ST := ST + 1
  501.         FSQRT              // ST := Sqrt(ST)
  502.         FMUL               // ST(1) := ST * ST(1); Pop ST
  503. @@3:    FWAIT
  504. end;
  505.  
  506.  
  507. procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
  508. asm
  509.         FLD     Theta
  510.         FSINCOS
  511.         FSTP    tbyte ptr [edx]    // Cos
  512.         FSTP    tbyte ptr [eax]    // Sin
  513.         FWAIT
  514. end;
  515.  
  516. { Extract exponent and mantissa from X }
  517. procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
  518. { Mantissa ptr in EAX, Exponent ptr in EDX }
  519. asm
  520.         FLD     X
  521.         PUSH    EAX
  522.         MOV     dword ptr [edx], 0    { if X = 0, return 0 }
  523.  
  524.         FTST
  525.         FSTSW   AX
  526.         FWAIT
  527.         SAHF
  528.         JZ      @@Done
  529.  
  530.         FXTRACT                 // ST(1) = exponent, (pushed) ST = fraction
  531.         FXCH
  532.  
  533. // The FXTRACT instruction normalizes the fraction 1 bit higher than
  534. // wanted for the definition of frexp() so we need to tweak the result
  535. // by scaling the fraction down and incrementing the exponent.
  536.  
  537.         FISTP   dword ptr [edx]
  538.         FLD1
  539.         FCHS
  540.         FXCH
  541.         FSCALE                  // scale fraction
  542.         INC     dword ptr [edx] // exponent biased to match
  543.         FSTP ST(1)              // discard -1, leave fraction as TOS
  544.  
  545. @@Done:
  546.         POP     EAX
  547.         FSTP    tbyte ptr [eax]
  548.         FWAIT
  549. end;
  550.  
  551. function Ldexp(X: Extended; P: Integer): Extended;
  552.   { Result := X * (2^P) }
  553. asm
  554.         PUSH    EAX
  555.         FILD    dword ptr [ESP]
  556.         FLD     X
  557.         FSCALE
  558.         POP     EAX
  559.         FSTP    ST(1)
  560.         FWAIT
  561. end;
  562.  
  563. function Ceil(X: Extended): Integer;
  564. begin
  565.   Result := Trunc(X);
  566.   if Frac(X) > 0 then
  567.     Inc(Result);
  568. end;
  569.  
  570. function Floor(X: Extended): Integer;
  571. begin
  572.   Result := Trunc(X);
  573.   if Frac(X) < 0 then
  574.     Dec(Result);
  575. end;
  576.  
  577. { Conversion of bases:  Log.b(X) = Log.a(X) / Log.a(b)  }
  578.  
  579. function Log10(X: Extended): Extended;
  580.   { Log.10(X) := Log.2(X) * Log.10(2) }
  581. asm
  582.         FLDLG2     { Log base ten of 2 }
  583.         FLD     X
  584.         FYL2X
  585.         FWAIT
  586. end;
  587.  
  588. function Log2(X: Extended): Extended;
  589. asm
  590.         FLD1
  591.         FLD     X
  592.         FYL2X
  593.         FWAIT
  594. end;
  595.  
  596. function LogN(Base, X: Extended): Extended;
  597. { Log.N(X) := Log.2(X) / Log.2(N) }
  598. asm
  599.         FLD1
  600.         FLD     X
  601.         FYL2X
  602.         FLD1
  603.         FLD     Base
  604.         FYL2X
  605.         FDIV
  606.         FWAIT
  607. end;
  608.  
  609. function Poly(X: Extended; const Coefficients: array of Double): Extended;
  610. { Horner's method }
  611. var
  612.   I: Integer;
  613. begin
  614.   Result := Coefficients[High(Coefficients)];
  615.   for I := High(Coefficients)-1 downto Low(Coefficients) do
  616.     Result := Result * X + Coefficients[I];
  617. end;
  618.  
  619. function Power(Base, Exponent: Extended): Extended;
  620. begin
  621.   { Make special cases of Exponent = 0 and Exponent = integer.
  622.     Error if Base < 0 and Exponent not integer. }
  623.   if Exponent = 0.0 then
  624.     Result := 1.0               { By fiat, 0**0 = 1 }
  625.   else if (Frac(Exponent) = 0.0) and (Exponent < MaxInt) then
  626.     Result := IntPower(Base, Trunc(Exponent))
  627.   else
  628.   begin
  629.     if Base < 0.0 then ArgError;
  630.     Result := Exp(Exponent * Ln(Base))
  631.   end
  632. end;
  633.  
  634. { Hyperbolic functions }
  635.  
  636. function CoshSinh(X: Extended; Factor: Double): Extended;
  637. begin
  638.   Result := Exp(X) / 2;
  639.   Result := Result + Factor / Result;
  640. end;
  641.  
  642. function Cosh(X: Extended): Extended;
  643. begin
  644.   Result := CoshSinh(X, 0.25)
  645. end;
  646.  
  647. function Sinh(X: Extended): Extended;
  648. begin
  649.   Result := CoshSinh(X, -0.25)
  650. end;
  651.  
  652. const
  653.   MaxTanhDomain = 5678.22249441322; // Ln(MaxExtended)/2
  654.  
  655. function Tanh(X: Extended): Extended;
  656. begin
  657.   if X > MaxTanhDomain then
  658.     Result := 1.0
  659.   else if X < -MaxTanhDomain then
  660.     Result := -1.0
  661.   else
  662.   begin
  663.     Result := Exp(X);
  664.     Result := Result * Result;
  665.     Result := (Result - 1.0) / (Result + 1.0)
  666.   end;
  667. end;
  668.  
  669. function ArcCosh(X: Extended): Extended;
  670. begin
  671.   if X <= 1.0 then
  672.     Result := 0.0
  673.   else if X > 1.0e10 then
  674.     Result := Ln(2) + Ln(X)
  675.   else
  676.     Result := Ln(X + Sqrt((X - 1.0) * (X + 1.0)));
  677. end;
  678.  
  679. function ArcSinh(X: Extended): Extended;
  680. var
  681.   Neg: Boolean;
  682. begin
  683.   if X = 0 then
  684.     Result := 0
  685.   else
  686.   begin
  687.     Neg := (X < 0);
  688.     X := Abs(X);
  689.     if X > 1.0e10 then
  690.       Result := Ln(2) + Ln(X)
  691.     else
  692.     begin
  693.       Result := X*X;
  694.       Result := LnXP1(X + Result / (1 + Sqrt(1 + Result)));
  695.     end;
  696.     if Neg then Result := -Result;
  697.   end;
  698. end;
  699.  
  700. function ArcTanh(X: Extended): Extended;
  701. var
  702.   Neg: Boolean;
  703. begin
  704.   if X = 0 then
  705.     Result := 0
  706.   else
  707.   begin
  708.     Neg := (X < 0);
  709.     X := Abs(X);
  710.     if X >= 1 then
  711.       Result := MaxExtended
  712.     else
  713.       Result := 0.5 * LnXP1((2.0 * X) / (1.0 - X));
  714.     if Neg then Result := -Result;
  715.   end;
  716. end;
  717.  
  718. { Statistical functions }
  719.  
  720. function Mean(const Data: array of Double): Extended;
  721. var
  722.   I: Integer;
  723. begin
  724.   Result := SUM(Data) / (High(Data) - Low(Data) + 1)
  725. end;
  726.  
  727. function MinValue(const Data: array of Double): Double;
  728. var
  729.   I: Integer;
  730. begin
  731.   Result := Data[Low(Data)];
  732.   for I := Low(Data) + 1 to High(Data) do
  733.     if Result > Data[I] then
  734.       Result := Data[I];
  735. end;
  736.  
  737. function MaxValue(const Data: array of Double): Double;
  738. var
  739.   I: Integer;
  740. begin
  741.   Result := Data[Low(Data)];
  742.   for I := Low(Data) + 1 to High(Data) do
  743.     if Result < Data[I] then
  744.       Result := Data[I];
  745. end;
  746.  
  747. procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
  748. var
  749.   Sums, Squares: Extended;
  750.   N: Integer;
  751. begin
  752.   SumsAndSquares(Data, Sums, Squares);
  753.   N := High(Data)- Low(Data) + 1;
  754.   Mean := Sums / N;
  755.   StdDev := Sqrt((Squares - Sqr(Sums) / N) / (N - 1));
  756. end;
  757.  
  758. procedure MomentSkewKurtosis(const Data: array of Double;
  759.   var M1, M2, M3, M4, Skew, Kurtosis: Extended);
  760. var
  761.   Sum, SumSquares, SumCubes, SumQuads, OverN, Accum, M1Sqr, S2N, S3N: Extended;
  762.   I: Integer;
  763. begin
  764.   OverN := 1 / (High(Data) - Low(Data) + 1);
  765.   Sum := 0;
  766.   SumSquares := 0;
  767.   SumCubes := 0;
  768.   SumQuads := 0;
  769.   for I := Low(Data) to High(Data) do
  770.   begin
  771.     Sum := Sum + Data[I];
  772.     Accum := Sqr(Data[I]);
  773.     SumSquares := SumSquares + Accum;
  774.     Accum := Accum*Data[I];
  775.     SumCubes := SumCubes + Accum;
  776.     SumQuads := SumQuads + Accum*Data[I];
  777.   end;
  778.   M1 := Sum * OverN;
  779.   M1Sqr := Sqr(M1);
  780.   S2N := SumSquares * OverN;
  781.   S3N := SumCubes * OverN;
  782.   M2 := S2N - M1Sqr;
  783.   M3 := S3N - (M1 * 3 * S2N) + 2*M1Sqr*M1;
  784.   M4 := (SumQuads * OverN) - (M1 * 4 * S3N) + (M1Sqr*6*S2N - 3*Sqr(M1Sqr));
  785.   Skew := M3 * Power(M2, -3/2);   // = M3 / Power(M2, 3/2)
  786.   Kurtosis := M4 / Sqr(M2);
  787. end;
  788.  
  789. function Norm(const Data: array of Double): Extended;
  790. begin
  791.   Result := Sqrt(SumOfSquares(Data));
  792. end;
  793.  
  794. function PopnStdDev(const Data: array of Double): Extended;
  795. begin
  796.   Result := Sqrt(PopnVariance(Data))
  797. end;
  798.  
  799. function PopnVariance(const Data: array of Double): Extended;
  800. begin
  801.   Result := TotalVariance(Data) / (High(Data) - Low(Data) + 1)
  802. end;
  803.  
  804. function RandG(Mean, StdDev: Extended): Extended;
  805. { Marsaglia-Bray algorithm }
  806. var
  807.   U1, S2: Extended;
  808. begin
  809.   repeat
  810.     U1 := 2*Random - 1;
  811.     S2 := Sqr(U1) + Sqr(2*Random-1);
  812.   until S2 < 1;
  813.   Result := Sqrt(-2*Ln(S2)/S2) * U1 * StdDev + Mean;
  814. end;
  815.  
  816. function StdDev(const Data: array of Double): Extended;
  817. begin
  818.   Result := Sqrt(Variance(Data))
  819. end;
  820.  
  821. function SUM(const Data: array of Double): Extended;
  822. {var
  823.   I: Integer;
  824. begin
  825.   Result := 0.0;
  826.   for I := Low(Data) to High(Data) do
  827.     Result := Result + Data[I]
  828. end; }
  829. asm  // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1
  830.      // Uses 2 accumulators to minimize read-after-write delays
  831.      // Est. 4 clocks per loop, 2 items per loop = 2 clocks per datum
  832.        FLD      qword ptr [EAX]    // Init first accumulator
  833.        INC      EDX
  834.        ADD      EAX,8
  835.        SHR      EDX,1
  836.        JC       @@OddCount
  837.        JZ       @@End
  838.        JMP      @@1
  839. @@OddCount:
  840.        FADD     qword ptr [EAX]
  841.        ADD      EAX,8
  842. @@1:   FLDZ                        // Init second accumulator
  843. @@2:   FADD     qword ptr [EAX]    // 1
  844.        FXCH     ST(1)              // 0
  845.        FADD     qword ptr [EAX+8]  // 1
  846.        FXCH     ST(1)              // 0
  847.        ADD      EAX,16             // 1
  848.        DEC      EDX                // 0
  849.        JNZ      @@2                // 1
  850.        FADD                        // ST(1) := ST + ST(1); Pop ST
  851. @@End:
  852.        FWAIT
  853. end;
  854.  
  855. function SumOfSquares(const Data: array of Double): Extended;
  856. var
  857.   I: Integer;
  858. begin
  859.   Result := 0.0;
  860.   for I := Low(Data) to High(Data) do
  861.     Result := Result + Sqr(Data[I]);
  862. end;
  863.  
  864. procedure SumsAndSquares(const Data: array of Double; var Sum, SumOfSquares: Extended);
  865. {var
  866.   I: Integer;
  867. begin
  868.   Sum := 0;
  869.   SumOfSquares := 0;
  870.   for I := Low(Data) to High(Data) do
  871.   begin
  872.     Sum := Sum + Data[I];
  873.     SumOfSquares := SumOfSquares + Data[I]*Data[I];
  874.   end;
  875. end;  }
  876. asm  // IN:  EAX = ptr to Data
  877.      //      EDX = High(Data) = Count - 1
  878.      //      ECX = ptr to Sum
  879.      // Est. 18 clocks per loop, 4 items per loop = 4.5 clocks per data item
  880.        PUSH    EBX
  881.        INC     EDX
  882.        FLDZ
  883.        MOV     EBX,EDX
  884.        FLD     ST(0)
  885.        AND     EBX,3
  886.        JZ      @@1
  887.        JMP     @@P1
  888. @@P0:  FADD
  889. @@P1:  FLD     qword ptr [EAX]
  890.        FADD    ST(2),ST
  891.        FMUL    ST,ST
  892.        ADD     EAX,8
  893.        DEC     EBX
  894.        JNZ      @@P0
  895.        FADD
  896. @@1:   SHR     EDX,2
  897.        JZ      @@End
  898.        FLDZ                         // FPU stack: Sqr2, Sqr1, Sum
  899.        JMP     @@3
  900. @@2:   FADD                         // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
  901. @@3:   FLD     qword ptr [EAX]      // Load Data1
  902.        FADD    ST(3),ST             // Sum := Sum + Data1
  903.        FMUL    ST,ST                // Data1 := Sqr(Data1)
  904.        FLD     qword ptr [EAX+8]    // Load Data2
  905.        FADD    ST(4),ST             // Sum := Sum + Data2
  906.        FMUL    ST,ST                // Data2 := Sqr(Data2)
  907.        FXCH                         // Move Sqr(Data1) into ST(0)
  908.        FADDP   ST(3),ST             // Sqr1 := Sqr1 + Sqr(Data1); Pop Data1
  909.        FLD     qword ptr [EAX+16]   // Load Data3
  910.        FADD    ST(4),ST             // Sum := Sum + Data3
  911.        FMUL    ST,ST                // Data3 := Sqr(Data3)
  912.        FXCH                         // Move Sqr(Data2) into ST(0)
  913.        FADDP   ST(3),ST             // Sqr1 := Sqr1 + Sqr(Data2); Pop Data2
  914.        FLD     qword ptr [EAX+24]   // Load Data4
  915.        FADD    ST(4),ST             // Sum := Sum + Data4
  916.        FMUL    ST,ST                // Sqr(Data4)
  917.        FXCH                         // Move Sqr(Data3) into ST(0)
  918.        FADDP   ST(3),ST             // Sqr1 := Sqr1 + Sqr(Data3); Pop Data3
  919.        ADD     EAX,32
  920.        DEC     EDX
  921.        JNZ     @@2
  922.        FADD                         // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
  923.        FADD                         // Sqr1 := Sqr2 + Sqr1; Pop Sqr2
  924. @@End:
  925.        MOV     EAX, SumOfSquares
  926.        FXCH                         // Move Sum1 into ST(0)
  927.        FSTP    tbyte ptr [ECX]      // Sum := Sum1; Pop Sum1
  928.        FSTP    tbyte ptr [EAX]      // SumOfSquares := Sum1; Pop Sum1
  929.        POP     EBX
  930.        FWAIT
  931. end;
  932.  
  933. function TotalVariance(const Data: array of Double): Extended;
  934. var
  935.   Sum, SumSquares: Extended;
  936. begin
  937.   SumsAndSquares(Data, Sum, SumSquares);
  938.   Result := SumSquares - Sqr(Sum)/(High(Data) - Low(Data) + 1);
  939. end;
  940.  
  941. function Variance(const Data: array of Double): Extended;
  942. begin
  943.   Result := TotalVariance(Data) / (High(Data) - Low(Data))
  944. end;
  945.  
  946.  
  947. { Depreciation functions. }
  948.  
  949. function DoubleDecliningBalance(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  950. { dv := cost * (1 - 2/life)**(period - 1)
  951.   DDB = (2/life) * dv
  952.   if DDB > dv - salvage then DDB := dv - salvage
  953.   if DDB < 0 then DDB := 0
  954. }
  955. var
  956.   DepreciatedVal, Factor: Extended;
  957. begin
  958.   if (Period < 1) or (Life < Period) or (Life <= 2) or (Cost < Salvage) then
  959.     ArgError;
  960.   Factor := 2.0 / Life;
  961.   DepreciatedVal := IntPower(Cost * (1.0 - Factor), Period - 1);
  962.   Result := Factor * DepreciatedVal;
  963.   if Result > DepreciatedVal - Salvage then
  964.     Result := DepreciatedVal - Salvage;
  965.   if Result < 0.0 then
  966.     Result := 0.0
  967. end;
  968.  
  969. function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
  970. { Spreads depreciation linearly over life. }
  971. begin
  972.   if Life < 1 then ArgError;
  973.   Result := (Cost - Salvage) / Life
  974. end;
  975.  
  976. function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
  977. { SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }
  978. var
  979.   X1, X2: Extended;
  980. begin
  981.   if (Period < 1) or (Life < Period) then ArgError;
  982.   X1 := 2 * (Life - Period + 1);
  983.   X2 := Life * (Life + 1);
  984.   Result := (Cost - Salvage) * X1 / X2
  985. end;
  986.  
  987. { Discounted cash flow functions. }
  988.  
  989. function InternalRateOfReturn(Guess: Extended; const CashFlows: array of Double): Extended;
  990. {
  991. Use Newton's method to solve NPV = 0, where NPV is a polynomial in
  992. x = 1/(1+rate).  Split the coefficients into negative and postive sets:
  993.   neg + pos = 0, so pos = -neg, so  -neg/pos = 1
  994. Then solve:
  995.   log(-neg/pos) = 0
  996.  
  997.   Let  t = log(1/(1+r) = -LnXP1(r)
  998.   then r = exp(-t) - 1
  999. Iterate on t, then use the last equation to compute r.
  1000. }
  1001. var
  1002.   T, Y: Extended;
  1003.   Poly: TPoly;
  1004.   K, Count: Integer;
  1005.  
  1006.   function ConditionP(const CashFlows: array of Double): Integer;
  1007.   { Guarantees existence and uniqueness of root.  The sign of payments
  1008.     must change exactly once, the net payout must be always > 0 for
  1009.     first portion, then each payment must be >= 0.
  1010.     Returns: 0 if condition not satisfied, > 0 if condition satisfied
  1011.     and this is the index of the first value considered a payback. }
  1012.   var
  1013.     X: Double;
  1014.     I, K: Integer;
  1015.   begin
  1016.     K := High(CashFlows);
  1017.     while (K >= 0) and (CashFlows[K] >= 0.0) do Dec(K);
  1018.     Inc(K);
  1019.     if K > 0 then begin
  1020.       X := 0.0;
  1021.       I := 0;
  1022.       while I < K do begin
  1023.         X := X + CashFlows[I];
  1024.         if X >= 0.0 then begin
  1025.           K := 0;
  1026.           Break
  1027.         end;
  1028.         Inc(I)
  1029.       end
  1030.     end;
  1031.     Result := K
  1032.   end;
  1033.  
  1034. begin
  1035.   Result := 0;
  1036.   K := ConditionP(CashFlows);
  1037.   if K < 0 then ArgError;
  1038.   if K = 0 then
  1039.   begin
  1040.     if Guess <= -1.0 then ArgError;
  1041.     T := -LnXP1(Guess)
  1042.   end else
  1043.     T := 0.0;
  1044.   for Count := 1 to MaxIterations do
  1045.   begin
  1046.     PolyX(CashFlows, Exp(T), Poly);
  1047.     if Poly.Pos <= Poly.Neg then ArgError;
  1048.     if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then
  1049.     begin
  1050.       Result := -1.0;
  1051.       Exit;
  1052.     end;
  1053.     with Poly do
  1054.       Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);
  1055.     T := T - Y;
  1056.     if RelSmall(Y, T) then
  1057.     begin
  1058.       Result := Exp(-T) - 1.0;
  1059.       Exit;
  1060.     end
  1061.   end;
  1062.   ArgError;
  1063. end;
  1064.  
  1065. function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  1066.   PaymentTime: TPaymentTime): Extended;
  1067. var
  1068.   X, Y: Extended;
  1069.   I: Integer;
  1070. begin
  1071.   if Rate <= -1.0 then ArgError;
  1072.   X := 1.0 / (1.0 + Rate);
  1073.   Y := CashFlows[High(CashFlows)];
  1074.   for I := High(CashFlows) - 1 downto Low(CashFlows) do
  1075.     Y := X * Y + CashFlows[I];
  1076.   if PaymentTime = ptStartOfPeriod then Y := X * Y;
  1077.   Result := Y;
  1078. end;
  1079.  
  1080. { Annuity functions. }
  1081.  
  1082. {---------------
  1083. From the point of view of A, amounts received by A are positive and
  1084. amounts disbursed by A are negative (e.g. a borrower's loan repayments
  1085. are regarded by the borrower as negative).
  1086.  
  1087. Given interest rate r, number of periods n:
  1088.   compound(r, n) = (1 + r)**n
  1089.   annuity(r, n) = (1 - compound(r, n)) / r
  1090.  
  1091. Given future value fv, periodic payment pmt, present value pv and type
  1092. of payment (start or end of period) pmtTime, financial variables satisfy:
  1093.  
  1094.   fv = pmt*(1 * r*pmtTime)*annuity(r, n) - pv*compound(r, n)
  1095.  
  1096. For fv, pv, pmt:
  1097.  
  1098.   C := compound(r, n)
  1099.   A := (1 + r*pmtTime)*annuity(r, n)
  1100.   Compute both at once in Annuity2.
  1101.  
  1102.   if C > 1E16 then A = -C/r, so:
  1103.     fv := meaningless
  1104.     pv := -(pmt/r + pmt*pmtTime)
  1105.     pmt := -pv*r/(1 + r*pmtTime)
  1106.   else
  1107.     fv := pmt*A - pv*C
  1108.     pv := (pmt*A - fv)/C
  1109.     pmt := (pv*C + fv)/A
  1110. ---------------}
  1111.  
  1112. function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,
  1113.   FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):
  1114.   Extended;
  1115. var
  1116.   Pmt, Rate1, T1, T2, T3: Extended;
  1117.   M: Integer;
  1118.  
  1119.   function Annuity(R: Extended; N: Integer): Extended;
  1120.   begin
  1121.     Result := (1 - Compound(R, N)) / R
  1122.   end;
  1123.  
  1124. begin
  1125.   M := Period;
  1126.   Rate1 := 0.0;
  1127.   if PaymentTime = ptEndOfPeriod then
  1128.   begin
  1129.     Inc(M);
  1130.     Rate1 := Rate
  1131.   end;
  1132.   T1 := (1 + Rate1) * Annuity(Rate, NPeriods);
  1133.   T2 := (1 + Rate1) * Annuity(Rate, M);
  1134.   T3 := (1 + Rate1) * Annuity(Rate, NPeriods - M);
  1135.   Pmt := (FutureValue + PresentValue * Compound(Rate, NPeriods)) / T1;
  1136.   IntPmt := Rate *
  1137.             (FutureValue * Compound(Rate, M) - PresentValue * T2 * T3) / T1;
  1138.   Result := Pmt - IntPmt
  1139. end;
  1140.  
  1141. function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  1142.   Extended; PaymentTime: TPaymentTime): Extended;
  1143. var
  1144.   Annuity, CompoundRN: Extended;
  1145. begin
  1146.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1147.   if CompoundRN > 1.0E16 then ArgError;
  1148.   Result := Payment * Annuity - PresentValue * CompoundRN
  1149. end;
  1150.  
  1151. function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  1152.   FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1153. begin
  1154.   if (Period < 1) or (Period > NPeriods) then ArgError;
  1155.   PaymentParts(Period, NPeriods, Rate, PresentValue, FutureValue, PaymentTime, Result);
  1156. end;
  1157.  
  1158. function InterestRate(NPeriods: Integer;
  1159.   Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1160. {
  1161. Given:
  1162.   First and last payments are non-zero and of opposite signs.
  1163.   Number of periods N >= 2.
  1164. Convert data into cash flow of first, N-1 payments, last with
  1165. first < 0, payment > 0, last > 0.
  1166. Compute the IRR of this cash flow:
  1167.   0 = first + pmt*x + pmt*x**2 + ... + pmt*x**(N-1) + last*x**N
  1168. where x = 1/(1 + rate).
  1169. Substitute x = exp(t) and apply Newton's method to
  1170.   f(t) = log(pmt*x + ... + last*x**N) / -first
  1171. which has a unique root given the above hypotheses.
  1172. }
  1173. var
  1174.   X, Y, Z, First, Pmt, Last, T, ET, EnT, ET1: Extended;
  1175.   Count: Integer;
  1176.   Reverse: Boolean;
  1177.  
  1178.   function LostPrecision(X: Extended): Boolean;
  1179.   asm
  1180.         XOR     EAX, EAX
  1181.         MOV     BX,WORD PTR X+8
  1182.         INC     EAX
  1183.         AND     EBX, $7FF0
  1184.         JZ      @@1
  1185.         CMP     EBX, $7FF0
  1186.         JE      @@1
  1187.         XOR     EAX,EAX
  1188.   @@1:
  1189.   end;
  1190.  
  1191. begin
  1192.   Result := 0;
  1193.   if NPeriods <= 0 then ArgError;
  1194.   Pmt := Payment;
  1195.   if PaymentTime = ptStartOfPeriod then
  1196.   begin
  1197.     X := PresentValue;
  1198.     Y := FutureValue + Payment
  1199.   end
  1200.   else
  1201.   begin
  1202.     X := PresentValue + Payment;
  1203.     Y := FutureValue
  1204.   end;
  1205.   First := X;
  1206.   Last := Y;
  1207.   Reverse := False;
  1208.   if First * Payment > 0.0 then
  1209.   begin
  1210.     Reverse := True;
  1211.     T := First;
  1212.     First := Last;
  1213.     Last := T
  1214.   end;
  1215.   if first > 0.0 then
  1216.   begin
  1217.     First := -First;
  1218.     Pmt := -Pmt;
  1219.     Last := -Last
  1220.   end;
  1221.   if (First = 0.0) or (Last < 0.0) then ArgError;
  1222.   T := 0.0;                     { Guess at solution }
  1223.   for Count := 1 to MaxIterations do
  1224.   begin
  1225.     EnT := Exp(NPeriods * T);
  1226.     if LostPrecision(EnT) then
  1227.     begin
  1228.       Result := -Pmt / First;
  1229.       if Reverse then
  1230.         Result := Exp(-LnXP1(Result)) - 1.0;
  1231.       Exit;
  1232.     end;
  1233.     ET := Exp(T);
  1234.     ET1 := ET - 1.0;
  1235.     if ET1 = 0.0 then
  1236.     begin
  1237.       X := NPeriods;
  1238.       Y := X * (X - 1.0) / 2.0
  1239.     end
  1240.     else
  1241.     begin
  1242.       X := ET * (Exp((NPeriods - 1) * T)-1.0) / ET1;
  1243.       Y := (NPeriods * EnT - ET - X * ET) / ET1
  1244.     end;
  1245.     Z := Pmt * X + Last * EnT;
  1246.     Y := Ln(Z / -First) / ((Pmt * Y + Last * NPeriods *EnT) / Z);
  1247.     T := T - Y;
  1248.     if RelSmall(Y, T) then
  1249.     begin
  1250.       if not Reverse then T := -T;
  1251.       Result := Exp(T)-1.0;
  1252.       Exit;
  1253.     end
  1254.   end;
  1255.   ArgError;
  1256. end;
  1257.  
  1258. function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  1259.   PaymentTime: TPaymentTime): Extended;
  1260. { If Rate = 0 then nper := -(pv + fv) / pmt
  1261.   else t := pv + pmt * (1 + rate*pmtTime) / rate
  1262.        nper := LnXP1(-(pv + vf) / t) / LnXP1(rate)
  1263. }
  1264. var
  1265.   T: Extended;
  1266. begin
  1267.   if Frac(Rate) = 0.0 then
  1268.     Result := -(PresentValue + FutureValue) / Payment
  1269.   else
  1270.   begin
  1271.     T := PresentValue + Payment * (1 + Integer(PaymentTime) * Rate) / Rate;
  1272.     Result := LnXP1(-(PresentValue + FutureValue) / T) / LnXP1(Rate)
  1273.   end
  1274. end;
  1275.  
  1276. function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:
  1277.   Extended; PaymentTime: TPaymentTime): Extended;
  1278. var
  1279.   Annuity, CompoundRN: Extended;
  1280. begin
  1281.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1282.   if CompoundRN > 1.0E16 then
  1283.     Result := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)
  1284.   else
  1285.     Result := (PresentValue * CompoundRN + FutureValue) / Annuity
  1286. end;
  1287.  
  1288. function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  1289.   PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
  1290. var
  1291.   Junk: Extended;
  1292. begin
  1293.   if (Period < 1) or (Period > NPeriods) then ArgError;
  1294.   Result := PaymentParts(Period, NPeriods, Rate, PresentValue, FutureValue,
  1295.                          PaymentTime, Junk);
  1296. end;
  1297.  
  1298. function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:
  1299.   Extended; PaymentTime: TPaymentTime): Extended;
  1300. var
  1301.   Annuity, CompoundRN: Extended;
  1302. begin
  1303.   Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
  1304.   if CompoundRN > 1.0E16 then
  1305.     Result := -(Payment / Rate * Integer(PaymentTime) * Payment)
  1306.   else
  1307.     Result := (Payment * Annuity - FutureValue) / CompoundRN
  1308. end;
  1309.  
  1310. end.
  1311.