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