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