| blob | 9ade6a6f43f9493d67e76126a755ee1049cee2f6 |
1 {=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
3 KKKKK KKKKK OOOOOOOOO LLLLL
4 KKKKK KKKKK OOOOOOOOOOOOO LLLLL
5 KKKKK KKKKK OOOOO OOOOO LLLLL
6 KKKKK KKKKK OOOOO OOOOO LLLLL
7 KKKKKKKKKK OOOOO OOOOO LLLLL
8 KKKKK KKKKK OOOOO OOOOO LLLLL
9 KKKKK KKKKK OOOOO OOOOO LLLLL
10 KKKKK KKKKK OOOOOOOOOOOOO LLLLLLLLLLLLL
11 KKKKK KKKKK OOOOOOOOO LLLLLLLLLLLLL
13 Key Objects Library (C) 2000 by Kladov Vladimir.
15 mailto: bonanzas@xcl.cjb.net
16 Home: http://kol.nm.ru
17 http://xcl.cjb.net
18 http://xcl.nm.ru
20 =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-}
21 {
22 This code is grabbed from standard math.pas unit,
23 provided by Borland Delphi. This unit is for working with
24 engineering (mathematical) functions. The main difference
25 is that err unit specially designed to handle exceptions
26 for KOL is used instead of SysUtils. This allows to make
27 size of the executable smaller for about 5K. though this
28 value is insignificant for project made with VCL, it can
29 be more than 15% of executable file size made with KOL.
30 }
32 {*******************************************************}
33 { }
34 { Borland Delphi Runtime Library }
35 { Math Unit }
36 { }
37 { Copyright (C) 1996,99 Inprise Corporation }
38 { }
39 {*******************************************************}
43 { This unit contains high-performance arithmetic, trigonometric, logorithmic,
44 statistical and financial calculation routines which supplement the math
45 routines that are part of the Delphi language or System unit. }
47 {$N+,S-}
49 {$I KOLDEF.INC}
51 interface
65 {-----------------------------------------------------------------------
66 References:
68 1) P.J. Plauger, "The Standard C Library", Prentice-Hall, 1992, Ch. 7.
69 2) W.J. Cody, Jr., and W. Waite, "Software Manual For the Elementary
70 Functions", Prentice-Hall, 1980.
71 3) Namir Shammas, "C/C++ Mathematical Algorithms for Scientists and Engineers",
72 McGraw-Hill, 1995, Ch 8.
73 4) H.T. Lau, "A Numerical Library in C for Scientists and Engineers",
74 CRC Press, 1994, Ch. 6.
75 5) "Pentium(tm) Processor User's Manual, Volume 3: Architecture
76 and Programming Manual", Intel, 1994
78 All angle parameters and results of trig functions are in radians.
80 Most of the following trig and log routines map directly to Intel 80387 FPU
81 floating point machine instructions. Input domains, output ranges, and
82 error handling are determined largely by the FPU hardware.
83 Routines coded in assembler favor the Pentium FPU pipeline architecture.
84 -----------------------------------------------------------------------}
92 { Trigonometric functions }
96 { ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
97 IN: |Y| < 2^64, |X| < 2^64, X <> 0 OUT: [-PI..PI] radians }
100 { SinCos is 2x faster than calling Sin and Cos separately for the same angle }
106 { Angle unit conversion routines }
114 { Hyperbolic functions and inverses }
122 { Logorithmic functions }
128 { Exponential functions }
130 { IntPower: Raise base to an integral power. Fast. }
131 //function IntPower(Base: Extended; Exponent: Integer): Extended register;
132 // -- already defined in kol.pas
134 { Power: Raise base to any power.
135 For fractional exponents, or |exponents| > MaxInt, base must be > 0. }
139 { Miscellaneous Routines }
141 { Frexp: Separates the mantissa and exponent of X. }
144 { Ldexp: returns X*2**P }
147 { Ceil: Smallest integer >= X, |X| < MaxInt }
150 { Floor: Largest integer <= X, |X| < MaxInt }
153 { Poly: Evaluates a uniform polynomial of one variable at value X.
154 The coefficients are ordered in increasing powers of X:
155 Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
158 {-----------------------------------------------------------------------
159 Statistical functions.
161 Common commercial spreadsheet macro names for these statistical and
162 financial functions are given in the comments preceding each function.
163 -----------------------------------------------------------------------}
165 { Mean: Arithmetic average of values. (AVG): SUM / N }
168 { Sum: Sum of values. (SUM) }
175 { MinValue: Returns the smallest signed value in the data array (MIN) }
180 {$IFDEF _D4orHigher}
181 overload;
186 {$ENDIF}
188 { MaxValue: Returns the largest signed value in the data array (MAX) }
193 {$IFDEF _D4orHigher}
194 overload;
199 {$ENDIF}
201 { Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
204 { MeanAndStdDev calculates Mean and StdDev in one call. }
207 { Population Standard Deviation (STDP): Sqrt(PopnVariance).
208 Used in some business and financial calculations. }
211 { Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
214 { Population Variance (VAR or VARP): TotalVariance/ N }
217 { Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
220 { Norm: The Euclidean L2-norm. Sqrt(SumOfSquares) }
223 { MomentSkewKurtosis: Calculates the core factors of statistical analysis:
224 the first four moments plus the coefficients of skewness and kurtosis.
225 M1 is the Mean. M2 is the Variance.
226 Skew reflects symmetry of distribution: M3 / (M2**(3/2))
227 Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
231 { RandG produces random numbers with Gaussian distribution about the mean.
232 Useful for simulating data with sampling errors. }
235 {-----------------------------------------------------------------------
236 Financial functions. Standard set from Quattro Pro.
238 Parameter conventions:
240 From the point of view of A, amounts received by A are positive and
241 amounts disbursed by A are negative (e.g. a borrower's loan repayments
242 are regarded by the borrower as negative).
244 Interest rates are per payment period. 11% annual percentage rate on a
245 loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
247 -----------------------------------------------------------------------}
249 type
252 { Double Declining Balance (DDB) }
256 { Future Value (FVAL) }
260 { Interest Payment (IPAYMT) }
264 { Interest Rate (IRATE) }
268 { Internal Rate of Return. (IRR) Needs array of cash flows. }
272 { Number of Periods (NPER) }
276 { Net Present Value. (NPV) Needs array of cash flows. }
280 { Payment (PAYMT) }
284 { Period Payment (PPAYMT) }
288 { Present Value (PVAL) }
292 { Straight Line depreciation (SLN) }
295 { Sum-of-Years-Digits depreciation (SYD) }
298 {type
299 EInvalidArgument = class(EMathError) end;}
301 implementation
303 {$IFNDEF _D2orD3}
305 {$ENDIF}
308 begin
316 begin
324 begin
332 begin
340 begin
351 type
356 const
360 begin
365 begin
370 begin
375 begin
380 begin
385 begin
390 begin
395 { Return ln(1 + X). Accurate for X near 0. }
396 asm
397 FLDLN2
399 FLD X
402 FLD1
403 FADD
404 FYL2X
406 @@1:
407 FYL2XP1
408 @@2:
409 FWAIT
410 end;
412 { Invariant: Y >= 0 & Result*X**Y = X**I. Init Y = I and Result = 1. }
413 {function IntPower(X: Extended; I: Integer): Extended;
414 var
415 Y: Integer;
416 begin
417 Y := Abs(I);
418 Result := 1.0;
419 while Y > 0 do begin
420 while not Odd(Y) do
421 begin
422 Y := Y shr 1;
423 X := X * X
424 end;
425 Dec(Y);
426 Result := Result * X
427 end;
428 if I < 0 then Result := 1.0 / Result
429 end;
430 }
431 (* -- already defined in kol.pas
432 function IntPower(Base: Extended; Exponent: Integer): Extended;
433 asm
434 mov ecx, eax
435 cdq
436 fld1 { Result := 1 }
437 xor eax, edx
438 sub eax, edx { eax := Abs(Exponent) }
439 jz @@3
440 fld Base
441 jmp @@2
442 @@1: fmul ST, ST { X := Base * Base }
443 @@2: shr eax,1
444 jnc @@1
445 fmul ST(1),ST { Result := Result * X }
446 jnz @@1
447 fstp st { pop X from FPU stack }
448 cmp ecx, 0
449 jge @@3
450 fld1
451 fdivrp { Result := 1 / Result }
452 @@3:
453 fwait
454 end;
455 *)
457 function Compound(R: Extended; N: Integer): Extended;
458 { Return (1 + R)**N. }
459 begin
460 Result := IntPower(1.0 + R, N)
461 end;
463 function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
464 var CompoundRN: Extended): Extended;
465 { Set CompoundRN to Compound(R, N),
466 return (1+Rate*PaymentTime)*(Compound(R,N)-1)/R;
467 }
468 begin
469 if R = 0.0 then
470 begin
471 CompoundRN := 1.0;
472 Result := N;
473 end
474 else
475 begin
476 { 6.1E-5 approx= 2**-14 }
477 if EAbs(R) < 6.1E-5 then
478 begin
479 CompoundRN := Exp(N * LnXP1(R));
480 Result := N*(1+(N-1)*R/2);
481 end
482 else
483 begin
484 CompoundRN := Compound(R, N);
485 Result := (CompoundRN-1) / R
486 end;
487 if PaymentTime = ptStartOfPeriod then
488 Result := Result * (1 + R);
489 end;
490 end; {Annuity2}
493 procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);
494 { Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.
495 Accumulate positive and negative terms separately. }
496 var
497 I: Integer;
498 Neg, Pos, DNeg, DPos: Extended;
499 begin
500 Neg := 0.0;
501 Pos := 0.0;
502 DNeg := 0.0;
503 DPos := 0.0;
504 for I := High(A) downto Low(A) do
505 begin
506 DNeg := X * DNeg + Neg;
507 Neg := Neg * X;
508 DPos := X * DPos + Pos;
509 Pos := Pos * X;
510 if A[I] >= 0.0 then
511 Pos := Pos + A[I]
512 else
513 Neg := Neg + A[I]
514 end;
515 Poly.Neg := Neg;
516 Poly.Pos := Pos;
517 Poly.DNeg := DNeg * X;
518 Poly.DPos := DPos * X;
519 end; {PolyX}
522 function RelSmall(X, Y: Extended): Boolean;
523 { Returns True if X is small relative to Y }
524 const
525 C1: Double = 1E-15;
526 C2: Double = 1E-12;
527 begin
528 Result := EAbs(X) < (C1 + C2 * EAbs(Y))
529 end;
531 { Math functions. }
533 function ArcCos(X: Extended): Extended;
534 begin
535 Result := ArcTan2(Sqrt(1 - X*X), X);
536 end;
538 function ArcSin(X: Extended): Extended;
539 begin
540 Result := ArcTan2(X, Sqrt(1 - X*X))
541 end;
543 function ArcTan2(Y, X: Extended): Extended;
544 asm
545 FLD Y
546 FLD X
547 FPATAN
548 FWAIT
549 end;
551 function Tan(X: Extended): Extended;
552 { Tan := Sin(X) / Cos(X) }
553 asm
554 FLD X
555 FPTAN
556 FSTP ST(0) { FPTAN pushes 1.0 after result }
557 FWAIT
558 end;
560 function CoTan(X: Extended): Extended;
561 { CoTan := Cos(X) / Sin(X) = 1 / Tan(X) }
562 asm
563 FLD X
564 FPTAN
565 FDIVRP
566 FWAIT
567 end;
569 function Hypot(X, Y: Extended): Extended;
570 { formula: Sqrt(X*X + Y*Y)
571 implemented as: |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision
572 var
573 Temp: Extended;
574 begin
575 X := Abs(X);
576 Y := Abs(Y);
577 if X > Y then
578 begin
579 Temp := X;
580 X := Y;
581 Y := Temp;
582 end;
583 if X = 0 then
584 Result := Y
585 else // Y > X, X <> 0, so Y > 0
586 Result := Y * Sqrt(1 + Sqr(X/Y));
587 end;
588 }
589 asm
590 FLD Y
591 FABS
592 FLD X
593 FABS
594 FCOM
595 FNSTSW AX
596 TEST AH,$45
597 JNZ @@1 // if ST > ST(1) then swap
598 FXCH ST(1) // put larger number in ST(1)
599 @@1: FLDZ
600 FCOMP
601 FNSTSW AX
602 TEST AH,$40 // if ST = 0, return ST(1)
603 JZ @@2
604 FSTP ST // eat ST(0)
605 JMP @@3
606 @@2: FDIV ST,ST(1) // ST := ST / ST(1)
607 FMUL ST,ST // ST := ST * ST
608 FLD1
609 FADD // ST := ST + 1
610 FSQRT // ST := Sqrt(ST)
611 FMUL // ST(1) := ST * ST(1); Pop ST
612 @@3: FWAIT
613 end;
616 procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
617 asm
618 FLD Theta
619 FSINCOS
620 FSTP tbyte ptr [edx] // Cos
621 FSTP tbyte ptr [eax] // Sin
622 FWAIT
623 end;
625 { Extract exponent and mantissa from X }
626 procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
627 { Mantissa ptr in EAX, Exponent ptr in EDX }
628 asm
629 FLD X
630 PUSH EAX
631 MOV dword ptr [edx], 0 { if X = 0, return 0 }
633 FTST
634 FSTSW AX
635 FWAIT
636 SAHF
637 JZ @@Done
639 FXTRACT // ST(1) = exponent, (pushed) ST = fraction
640 FXCH
642 // The FXTRACT instruction normalizes the fraction 1 bit higher than
643 // wanted for the definition of frexp() so we need to tweak the result
644 // by scaling the fraction down and incrementing the exponent.
646 FISTP dword ptr [edx]
647 FLD1
648 FCHS
649 FXCH
650 FSCALE // scale fraction
651 INC dword ptr [edx] // exponent biased to match
652 FSTP ST(1) // discard -1, leave fraction as TOS
654 @@Done:
655 POP EAX
656 FSTP tbyte ptr [eax]
657 FWAIT
658 end;
660 function Ldexp(X: Extended; P: Integer): Extended;
661 { Result := X * (2^P) }
662 asm
663 PUSH EAX
664 FILD dword ptr [ESP]
665 FLD X
666 FSCALE
667 POP EAX
668 FSTP ST(1)
669 FWAIT
670 end;
672 function Ceil(X: Extended): Integer;
673 begin
674 Result := Integer(Trunc(X));
675 if Frac(X) > 0 then
676 Inc(Result);
677 end;
679 function Floor(X: Extended): Integer;
680 begin
681 Result := Integer(Trunc(X));
682 if Frac(X) < 0 then
683 Dec(Result);
684 end;
686 { Conversion of bases: Log.b(X) = Log.a(X) / Log.a(b) }
688 function Log10(X: Extended): Extended;
689 { Log.10(X) := Log.2(X) * Log.10(2) }
690 asm
691 FLDLG2 { Log base ten of 2 }
692 FLD X
693 FYL2X
694 FWAIT
695 end;
697 function Log2(X: Extended): Extended;
698 asm
699 FLD1
700 FLD X
701 FYL2X
702 FWAIT
703 end;
705 function LogN(Base, X: Extended): Extended;
706 { Log.N(X) := Log.2(X) / Log.2(N) }
707 asm
708 FLD1
709 FLD X
710 FYL2X
711 FLD1
712 FLD Base
713 FYL2X
714 FDIV
715 FWAIT
716 end;
718 function Poly(X: Extended; const Coefficients: array of Double): Extended;
719 { Horner's method }
720 var
721 I: Integer;
722 begin
723 Result := Coefficients[High(Coefficients)];
724 for I := High(Coefficients)-1 downto Low(Coefficients) do
725 Result := Result * X + Coefficients[I];
726 end;
728 function Power(Base, Exponent: Extended): Extended;
729 begin
730 if Exponent = 0.0 then
731 Result := 1.0 { n**0 = 1 }
732 else if (Base = 0.0) and (Exponent > 0.0) then
733 Result := 0.0 { 0**n = 0, n > 0 }
734 else if (Frac(Exponent) = 0.0) and (EAbs(Exponent) <= MaxInt) then
735 Result := IntPower(Base, Integer(Trunc(Exponent)))
736 else
737 Result := Exp(Exponent * Ln(Base))
738 end;
740 { Hyperbolic functions }
742 function CoshSinh(X: Extended; Factor: Double): Extended;
743 begin
744 Result := Exp(X) / 2;
745 Result := Result + Factor / Result;
746 end;
748 function Cosh(X: Extended): Extended;
749 begin
750 Result := CoshSinh(X, 0.25)
751 end;
753 function Sinh(X: Extended): Extended;
754 begin
755 Result := CoshSinh(X, -0.25)
756 end;
758 const
759 MaxTanhDomain = 5678.22249441322; // Ln(MaxExtended)/2
761 function Tanh(X: Extended): Extended;
762 begin
763 if X > MaxTanhDomain then
764 Result := 1.0
765 else if X < -MaxTanhDomain then
766 Result := -1.0
767 else
768 begin
769 Result := Exp(X);
770 Result := Result * Result;
771 Result := (Result - 1.0) / (Result + 1.0)
772 end;
773 end;
775 function ArcCosh(X: Extended): Extended;
776 begin
777 if X <= 1.0 then
778 Result := 0.0
779 else if X > 1.0e10 then
780 Result := Ln(2) + Ln(X)
781 else
782 Result := Ln(X + Sqrt((X - 1.0) * (X + 1.0)));
783 end;
785 function ArcSinh(X: Extended): Extended;
786 var
787 Neg: Boolean;
788 begin
789 if X = 0 then
790 Result := 0
791 else
792 begin
793 Neg := (X < 0);
794 X := EAbs(X);
795 if X > 1.0e10 then
796 Result := Ln(2) + Ln(X)
797 else
798 begin
799 Result := X*X;
800 Result := LnXP1(X + Result / (1 + Sqrt(1 + Result)));
801 end;
802 if Neg then Result := -Result;
803 end;
804 end;
806 function ArcTanh(X: Extended): Extended;
807 var
808 Neg: Boolean;
809 begin
810 if X = 0 then
811 Result := 0
812 else
813 begin
814 Neg := (X < 0);
815 X := EAbs(X);
816 if X >= 1 then
817 Result := MaxExtended
818 else
819 Result := 0.5 * LnXP1((2.0 * X) / (1.0 - X));
820 if Neg then Result := -Result;
821 end;
822 end;
824 { Statistical functions }
826 function Mean(const Data: array of Double): Extended;
827 begin
828 Result := SUM(Data) / (High(Data) - Low(Data) + 1)
829 end;
831 function MinValue(const Data: array of Double): Double;
832 var
833 I: Integer;
834 begin
835 Result := Data[Low(Data)];
836 for I := Low(Data) + 1 to High(Data) do
837 if Result > Data[I] then
838 Result := Data[I];
839 end;
841 function MinIntValue(const Data: array of Integer): Integer;
842 var
843 I: Integer;
844 begin
845 Result := Data[Low(Data)];
846 for I := Low(Data) + 1 to High(Data) do
847 if Result > Data[I] then
848 Result := Data[I];
849 end;
851 function Min(A,B: Integer): Integer;
852 begin
853 if A < B then
854 Result := A
855 else
856 Result := B;
857 end;
859 {$IFDEF _D4orHigher}
860 function Min(A,B: I64): I64;
861 begin
862 if Cmp64( A, B ) < 0 then
863 Result := A
864 else
865 Result := B;
866 end;
868 function Min(A,B: Single): Single;
869 begin
870 if A < B then
871 Result := A
872 else
873 Result := B;
874 end;
876 function Min(A,B: Double): Double;
877 begin
878 if A < B then
879 Result := A
880 else
881 Result := B;
882 end;
884 function Min(A,B: Extended): Extended;
885 begin
886 if A < B then
887 Result := A
888 else
889 Result := B;
890 end;
891 {$ENDIF}
893 function MaxValue(const Data: array of Double): Double;
894 var
895 I: Integer;
896 begin
897 Result := Data[Low(Data)];
898 for I := Low(Data) + 1 to High(Data) do
899 if Result < Data[I] then
900 Result := Data[I];
901 end;
903 function MaxIntValue(const Data: array of Integer): Integer;
904 var
905 I: Integer;
906 begin
907 Result := Data[Low(Data)];
908 for I := Low(Data) + 1 to High(Data) do
909 if Result < Data[I] then
910 Result := Data[I];
911 end;
913 function Max(A,B: Integer): Integer;
914 begin
915 if A > B then
916 Result := A
917 else
918 Result := B;
919 end;
921 {$IFDEF _D4orHigher}
922 function Max(A,B: I64): I64;
923 begin
924 if Cmp64( A, B ) > 0 then
925 Result := A
926 else
927 Result := B;
928 end;
930 function Max(A,B: Single): Single;
931 begin
932 if A > B then
933 Result := A
934 else
935 Result := B;
936 end;
938 function Max(A,B: Double): Double;
939 begin
940 if A > B then
941 Result := A
942 else
943 Result := B;
944 end;
946 function Max(A,B: Extended): Extended;
947 begin
948 if A > B then
949 Result := A
950 else
951 Result := B;
952 end;
953 {$ENDIF}
955 procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
956 var
957 S: Extended;
958 N,I: Integer;
959 begin
960 N := High(Data)- Low(Data) + 1;
961 if N = 1 then
962 begin
963 Mean := Data[0];
964 StdDev := Data[0];
965 Exit;
966 end;
967 Mean := Sum(Data) / N;
968 S := 0; // sum differences from the mean, for greater accuracy
969 for I := Low(Data) to High(Data) do
970 S := S + Sqr(Mean - Data[I]);
971 StdDev := Sqrt(S / (N - 1));
972 end;
974 procedure MomentSkewKurtosis(const Data: array of Double;
975 var M1, M2, M3, M4, Skew, Kurtosis: Extended);
976 var
977 Sum, SumSquares, SumCubes, SumQuads, OverN, Accum, M1Sqr, S2N, S3N: Extended;
978 I: Integer;
979 begin
980 OverN := 1 / (High(Data) - Low(Data) + 1);
981 Sum := 0;
982 SumSquares := 0;
983 SumCubes := 0;
984 SumQuads := 0;
985 for I := Low(Data) to High(Data) do
986 begin
987 Sum := Sum + Data[I];
988 Accum := Sqr(Data[I]);
989 SumSquares := SumSquares + Accum;
990 Accum := Accum*Data[I];
991 SumCubes := SumCubes + Accum;
992 SumQuads := SumQuads + Accum*Data[I];
993 end;
994 M1 := Sum * OverN;
995 M1Sqr := Sqr(M1);
996 S2N := SumSquares * OverN;
997 S3N := SumCubes * OverN;
998 M2 := S2N - M1Sqr;
999 M3 := S3N - (M1 * 3 * S2N) + 2*M1Sqr*M1;
1000 M4 := (SumQuads * OverN) - (M1 * 4 * S3N) + (M1Sqr*6*S2N - 3*Sqr(M1Sqr));
1001 Skew := M3 * Power(M2, -3/2); // = M3 / Power(M2, 3/2)
1002 Kurtosis := M4 / Sqr(M2);
1003 end;
1005 function Norm(const Data: array of Double): Extended;
1006 begin
1007 Result := Sqrt(SumOfSquares(Data));
1008 end;
1010 function PopnStdDev(const Data: array of Double): Extended;
1011 begin
1012 Result := Sqrt(PopnVariance(Data))
1013 end;
1015 function PopnVariance(const Data: array of Double): Extended;
1016 begin
1017 Result := TotalVariance(Data) / (High(Data) - Low(Data) + 1)
1018 end;
1020 function RandG(Mean, StdDev: Extended): Extended;
1021 { Marsaglia-Bray algorithm }
1022 var
1023 U1, S2: Extended;
1024 begin
1025 repeat
1026 U1 := 2*Random - 1;
1027 S2 := Sqr(U1) + Sqr(2*Random-1);
1028 until S2 < 1;
1029 Result := Sqrt(-2*Ln(S2)/S2) * U1 * StdDev + Mean;
1030 end;
1032 function StdDev(const Data: array of Double): Extended;
1033 begin
1034 Result := Sqrt(Variance(Data))
1035 end;
1037 procedure RaiseOverflowError; forward;
1039 function SumInt(const Data: array of Integer): Integer;
1040 {var
1041 I: Integer;
1042 begin
1043 Result := 0;
1044 for I := Low(Data) to High(Data) do
1045 Result := Result + Data[I]
1046 end; }
1047 asm // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1
1048 // loop unrolled 4 times, 5 clocks per loop, 1.2 clocks per datum
1049 PUSH EBX
1050 MOV ECX, EAX // ecx = ptr to data
1051 MOV EBX, EDX
1052 XOR EAX, EAX
1053 AND EDX, not 3
1054 AND EBX, 3
1055 SHL EDX, 2
1056 JMP @Vector.Pointer[EBX*4]
1057 @Vector:
1058 DD @@1
1059 DD @@2
1060 DD @@3
1061 DD @@4
1062 @@4:
1063 ADD EAX, [ECX+12+EDX]
1064 JO RaiseOverflowError
1065 @@3:
1066 ADD EAX, [ECX+8+EDX]
1067 JO RaiseOverflowError
1068 @@2:
1069 ADD EAX, [ECX+4+EDX]
1070 JO RaiseOverflowError
1071 @@1:
1072 ADD EAX, [ECX+EDX]
1073 JO RaiseOverflowError
1074 SUB EDX,16
1075 JNS @@4
1076 POP EBX
1077 end;
1079 procedure RaiseOverflowError;
1080 begin
1081 raise Exception.Create(e_IntOverflow, SIntOverflow);
1082 end;
1084 function SUM(const Data: array of Double): Extended;
1085 {var
1086 I: Integer;
1087 begin
1088 Result := 0.0;
1089 for I := Low(Data) to High(Data) do
1090 Result := Result + Data[I]
1091 end; }
1092 asm // IN: EAX = ptr to Data, EDX = High(Data) = Count - 1
1093 // Uses 4 accumulators to minimize read-after-write delays and loop overhead
1094 // 5 clocks per loop, 4 items per loop = 1.2 clocks per item
1095 FLDZ
1096 MOV ECX, EDX
1097 FLD ST(0)
1098 AND EDX, not 3
1099 FLD ST(0)
1100 AND ECX, 3
1101 FLD ST(0)
1102 SHL EDX, 3 // count * sizeof(Double) = count * 8
1103 JMP @Vector.Pointer[ECX*4]
1104 @Vector:
1105 DD @@1
1106 DD @@2
1107 DD @@3
1108 DD @@4
1109 @@4: FADD qword ptr [EAX+EDX+24] // 1
1110 FXCH ST(3) // 0
1111 @@3: FADD qword ptr [EAX+EDX+16] // 1
1112 FXCH ST(2) // 0
1113 @@2: FADD qword ptr [EAX+EDX+8] // 1
1114 FXCH ST(1) // 0
1115 @@1: FADD qword ptr [EAX+EDX] // 1
1116 FXCH ST(2) // 0
1117 SUB EDX, 32
1118 JNS @@4
1119 FADDP ST(3),ST // ST(3) := ST + ST(3); Pop ST
1120 FADD // ST(1) := ST + ST(1); Pop ST
1121 FADD // ST(1) := ST + ST(1); Pop ST
1122 FWAIT
1123 end;
1125 function SumOfSquares(const Data: array of Double): Extended;
1126 var
1127 I: Integer;
1128 begin
1129 Result := 0.0;
1130 for I := Low(Data) to High(Data) do
1131 Result := Result + Sqr(Data[I]);
1132 end;
1134 procedure SumsAndSquares(const Data: array of Double; var Sum, SumOfSquares: Extended);
1135 {var
1136 I: Integer;
1137 begin
1138 Sum := 0;
1139 SumOfSquares := 0;
1140 for I := Low(Data) to High(Data) do
1141 begin
1142 Sum := Sum + Data[I];
1143 SumOfSquares := SumOfSquares + Data[I]*Data[I];
1144 end;
1145 end; }
1146 asm // IN: EAX = ptr to Data
1147 // EDX = High(Data) = Count - 1
1148 // ECX = ptr to Sum
1149 // Est. 17 clocks per loop, 4 items per loop = 4.5 clocks per data item
1150 FLDZ // init Sum accumulator
1151 PUSH ECX
1152 MOV ECX, EDX
1153 FLD ST(0) // init Sqr1 accum.
1154 AND EDX, not 3
1155 FLD ST(0) // init Sqr2 accum.
1156 AND ECX, 3
1157 FLD ST(0) // init/simulate last data item left in ST
1158 SHL EDX, 3 // count * sizeof(Double) = count * 8
1159 JMP @Vector.Pointer[ECX*4]
1160 @Vector:
1161 DD @@1
1162 DD @@2
1163 DD @@3
1164 DD @@4
1165 @@4: FADD // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
1166 FLD qword ptr [EAX+EDX+24] // Load Data1
1167 FADD ST(3),ST // Sum := Sum + Data1
1168 FMUL ST,ST // Data1 := Sqr(Data1)
1169 @@3: FLD qword ptr [EAX+EDX+16] // Load Data2
1170 FADD ST(4),ST // Sum := Sum + Data2
1171 FMUL ST,ST // Data2 := Sqr(Data2)
1172 FXCH // Move Sqr(Data1) into ST(0)
1173 FADDP ST(3),ST // Sqr1 := Sqr1 + Sqr(Data1); Pop Data1
1174 @@2: FLD qword ptr [EAX+EDX+8] // Load Data3
1175 FADD ST(4),ST // Sum := Sum + Data3
1176 FMUL ST,ST // Data3 := Sqr(Data3)
1177 FXCH // Move Sqr(Data2) into ST(0)
1178 FADDP ST(3),ST // Sqr1 := Sqr1 + Sqr(Data2); Pop Data2
1179 @@1: FLD qword ptr [EAX+EDX] // Load Data4
1180 FADD ST(4),ST // Sum := Sum + Data4
1181 FMUL ST,ST // Sqr(Data4)
1182 FXCH // Move Sqr(Data3) into ST(0)
1183 FADDP ST(3),ST // Sqr1 := Sqr1 + Sqr(Data3); Pop Data3
1184 SUB EDX,32
1185 JNS @@4
1186 FADD // Sqr2 := Sqr2 + Sqr(Data4); Pop Data4
1187 POP ECX
1188 FADD // Sqr1 := Sqr2 + Sqr1; Pop Sqr2
1189 FXCH // Move Sum1 into ST(0)
1190 MOV EAX, SumOfSquares
1191 FSTP tbyte ptr [ECX] // Sum := Sum1; Pop Sum1
1192 FSTP tbyte ptr [EAX] // SumOfSquares := Sum1; Pop Sum1
1193 FWAIT
1194 end;
1196 function TotalVariance(const Data: array of Double): Extended;
1197 var
1198 Sum, SumSquares: Extended;
1199 begin
1200 SumsAndSquares(Data, Sum, SumSquares);
1201 Result := SumSquares - Sqr(Sum)/(High(Data) - Low(Data) + 1);
1202 end;
1204 function Variance(const Data: array of Double): Extended;
1205 begin
1206 Result := TotalVariance(Data) / (High(Data) - Low(Data))
1207 end;
1210 { Depreciation functions. }
1212 function DoubleDecliningBalance(Cost, Salvage: Extended; Life, Period: Integer): Extended;
1213 { dv := cost * (1 - 2/life)**(period - 1)
1214 DDB = (2/life) * dv
1215 if DDB > dv - salvage then DDB := dv - salvage
1216 if DDB < 0 then DDB := 0
1217 }
1218 var
1219 DepreciatedVal, Factor: Extended;
1220 begin
1221 Result := 0;
1222 if (Period < 1) or (Life < Period) or (Life < 1) or (Cost <= Salvage) then
1223 Exit;
1225 {depreciate everything in period 1 if life is only one or two periods}
1226 if ( Life <= 2 ) then
1227 begin
1228 if ( Period = 1 ) then
1229 DoubleDecliningBalance:=Cost-Salvage
1230 else
1231 DoubleDecliningBalance:=0; {all depreciation occurred in first period}
1232 exit;
1233 end;
1234 Factor := 2.0 / Life;
1236 DepreciatedVal := Cost * IntPower((1.0 - Factor), Period - 1);
1237 {DepreciatedVal is Cost-(sum of previous depreciation results)}
1239 Result := Factor * DepreciatedVal;
1240 {Nominal computed depreciation for this period. The rest of the
1241 function applies limits to this nominal value. }
1243 {Only depreciate until total depreciation equals cost-salvage.}
1244 if Result > DepreciatedVal - Salvage then
1245 Result := DepreciatedVal - Salvage;
1247 {No more depreciation after salvage value is reached. This is mostly a nit.
1248 If Result is negative at this point, it's very close to zero.}
1249 if Result < 0.0 then
1250 Result := 0.0;
1251 end;
1253 function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
1254 { Spreads depreciation linearly over life. }
1255 begin
1256 if Life < 1 then ArgError('SLNDepreciation');
1257 Result := (Cost - Salvage) / Life
1258 end;
1260 function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
1261 { SYD = (cost - salvage) * (life - period + 1) / (life*(life + 1)/2) }
1262 { Note: life*(life+1)/2 = 1+2+3+...+life "sum of years"
1263 The depreciation factor varies from life/sum_of_years in first period = 1
1264 downto 1/sum_of_years in last period = life.
1265 Total depreciation over life is cost-salvage.}
1266 var
1267 X1, X2: Extended;
1268 begin
1269 Result := 0;
1270 if (Period < 1) or (Life < Period) or (Cost <= Salvage) then Exit;
1271 X1 := 2 * (Life - Period + 1);
1272 X2 := Life * (Life + 1);
1273 Result := (Cost - Salvage) * X1 / X2
1274 end;
1276 { Discounted cash flow functions. }
1278 function InternalRateOfReturn(Guess: Extended; const CashFlows: array of Double): Extended;
1279 {
1280 Use Newton's method to solve NPV = 0, where NPV is a polynomial in
1281 x = 1/(1+rate). Split the coefficients into negative and postive sets:
1282 neg + pos = 0, so pos = -neg, so -neg/pos = 1
1283 Then solve:
1284 log(-neg/pos) = 0
1286 Let t = log(1/(1+r) = -LnXP1(r)
1287 then r = exp(-t) - 1
1288 Iterate on t, then use the last equation to compute r.
1289 }
1290 var
1291 T, Y: Extended;
1292 Poly: TPoly;
1293 K, Count: Integer;
1295 function ConditionP(const CashFlows: array of Double): Integer;
1296 { Guarantees existence and uniqueness of root. The sign of payments
1297 must change exactly once, the net payout must be always > 0 for
1298 first portion, then each payment must be >= 0.
1299 Returns: 0 if condition not satisfied, > 0 if condition satisfied
1300 and this is the index of the first value considered a payback. }
1301 var
1302 X: Double;
1303 I, K: Integer;
1304 begin
1305 K := High(CashFlows);
1306 while (K >= 0) and (CashFlows[K] >= 0.0) do Dec(K);
1307 Inc(K);
1308 if K > 0 then
1309 begin
1310 X := 0.0;
1311 I := 0;
1312 while I < K do begin
1313 X := X + CashFlows[I];
1314 if X >= 0.0 then
1315 begin
1316 K := 0;
1317 Break
1318 end;
1319 Inc(I)
1320 end
1321 end;
1322 ConditionP := K
1323 end;
1325 begin
1326 InternalRateOfReturn := 0;
1327 K := ConditionP(CashFlows);
1328 if K < 0 then ArgError('InternalRateOfReturn');
1329 if K = 0 then
1330 begin
1331 if Guess <= -1.0 then ArgError('InternalRateOfReturn');
1332 T := -LnXP1(Guess)
1333 end else
1334 T := 0.0;
1335 for Count := 1 to MaxIterations do
1336 begin
1337 PolyX(CashFlows, Exp(T), Poly);
1338 if Poly.Pos <= Poly.Neg then ArgError('InternalRateOfReturn');
1339 if (Poly.Neg >= 0.0) or (Poly.Pos <= 0.0) then
1340 begin
1341 InternalRateOfReturn := -1.0;
1342 Exit;
1343 end;
1344 with Poly do
1345 Y := Ln(-Neg / Pos) / (DNeg / Neg - DPos / Pos);
1346 T := T - Y;
1347 if RelSmall(Y, T) then
1348 begin
1349 InternalRateOfReturn := Exp(-T) - 1.0;
1350 Exit;
1351 end
1352 end;
1353 ArgError('InternalRateOfReturn');
1354 end;
1356 function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
1357 PaymentTime: TPaymentTime): Extended;
1358 { Caution: The sign of NPV is reversed from what would be expected for standard
1359 cash flows!}
1360 var
1361 rr: Extended;
1362 I: Integer;
1363 begin
1364 if Rate <= -1.0 then ArgError('NetPresentValue');
1365 rr := 1/(1+Rate);
1366 result := 0;
1367 for I := High(CashFlows) downto Low(CashFlows) do
1368 result := rr * result + CashFlows[I];
1369 if PaymentTime = ptEndOfPeriod then result := rr * result;
1370 end;
1372 { Annuity functions. }
1374 {---------------
1375 From the point of view of A, amounts received by A are positive and
1376 amounts disbursed by A are negative (e.g. a borrower's loan repayments
1377 are regarded by the borrower as negative).
1379 Given interest rate r, number of periods n:
1380 compound(r, n) = (1 + r)**n "Compounding growth factor"
1381 annuity(r, n) = (compound(r, n)-1) / r "Annuity growth factor"
1383 Given future value fv, periodic payment pmt, present value pv and type
1384 of payment (start, 1 , or end of period, 0) pmtTime, financial variables satisfy:
1386 fv = -pmt*(1 + r*pmtTime)*annuity(r, n) - pv*compound(r, n)
1388 For fv, pv, pmt:
1390 C := compound(r, n)
1391 A := (1 + r*pmtTime)*annuity(r, n)
1392 Compute both at once in Annuity2.
1394 if C > 1E16 then A = C/r, so:
1395 fv := meaningless
1396 pv := -pmt*(pmtTime+1/r)
1397 pmt := -pv*r/(1 + r*pmtTime)
1398 else
1399 fv := -pmt(1+r*pmtTime)*A - pv*C
1400 pv := (-pmt(1+r*pmtTime)*A - fv)/C
1401 pmt := (-pv*C-fv)/((1+r*pmtTime)*A)
1402 ---------------}
1404 function PaymentParts(Period, NPeriods: Integer; Rate, PresentValue,
1405 FutureValue: Extended; PaymentTime: TPaymentTime; var IntPmt: Extended):
1406 Extended;
1407 var
1408 Crn:extended; { =Compound(Rate,NPeriods) }
1409 Crp:extended; { =Compound(Rate,Period-1) }
1410 Arn:extended; { =AnnuityF(Rate,NPeriods) }
1412 begin
1413 if Rate <= -1.0 then ArgError('PaymentParts');
1414 Crp:=Compound(Rate,Period-1);
1415 Arn:=Annuity2(Rate,NPeriods,PaymentTime,Crn);
1416 IntPmt:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;
1417 PaymentParts:=(-FutureValue-PresentValue)*Crp/Arn;
1418 end;
1420 function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
1421 Extended; PaymentTime: TPaymentTime): Extended;
1422 var
1423 Annuity, CompoundRN: Extended;
1424 begin
1425 if Rate <= -1.0 then ArgError('FutureValue');
1426 Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
1427 if CompoundRN > 1.0E16 then ArgError('FutureValue');
1428 FutureValue := -Payment * Annuity - PresentValue * CompoundRN
1429 end;
1431 function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
1432 FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
1433 var
1434 Crp:extended; { compound(rate,period-1)}
1435 Crn:extended; { compound(rate,nperiods)}
1436 Arn:extended; { annuityf(rate,nperiods)}
1437 begin
1438 if (Rate <= -1.0)
1439 or (Period < 1) or (Period > NPeriods) then ArgError('InterestPayment');
1440 Crp:=Compound(Rate,Period-1);
1441 Arn:=Annuity2(Rate,Nperiods,PaymentTime,Crn);
1442 InterestPayment:=(FutureValue*(Crp-1)-PresentValue*(Crn-Crp))/Arn;
1443 end;
1445 function InterestRate(NPeriods: Integer;
1446 Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
1447 {
1448 Given:
1449 First and last payments are non-zero and of opposite signs.
1450 Number of periods N >= 2.
1451 Convert data into cash flow of first, N-1 payments, last with
1452 first < 0, payment > 0, last > 0.
1453 Compute the IRR of this cash flow:
1454 0 = first + pmt*x + pmt*x**2 + ... + pmt*x**(N-1) + last*x**N
1455 where x = 1/(1 + rate).
1456 Substitute x = exp(t) and apply Newton's method to
1457 f(t) = log(pmt*x + ... + last*x**N) / -first
1458 which has a unique root given the above hypotheses.
1459 }
1460 var
1461 X, Y, Z, First, Pmt, Last, T, ET, EnT, ET1: Extended;
1462 Count: Integer;
1463 Reverse: Boolean;
1465 function LostPrecision(X: Extended): Boolean;
1466 asm
1467 XOR EAX, EAX
1468 MOV BX,WORD PTR X+8
1469 INC EAX
1470 AND EBX, $7FF0
1471 JZ @@1
1472 CMP EBX, $7FF0
1473 JE @@1
1474 XOR EAX,EAX
1475 @@1:
1476 end;
1478 begin
1479 Result := 0;
1480 if NPeriods <= 0 then ArgError('InterestRate');
1481 Pmt := Payment;
1482 if PaymentTime = ptEndOfPeriod then
1483 begin
1484 X := PresentValue;
1485 Y := FutureValue + Payment
1486 end
1487 else
1488 begin
1489 X := PresentValue + Payment;
1490 Y := FutureValue
1491 end;
1492 First := X;
1493 Last := Y;
1494 Reverse := False;
1495 if First * Payment > 0.0 then
1496 begin
1497 Reverse := True;
1498 T := First;
1499 First := Last;
1500 Last := T
1501 end;
1502 if first > 0.0 then
1503 begin
1504 First := -First;
1505 Pmt := -Pmt;
1506 Last := -Last
1507 end;
1508 if (First = 0.0) or (Last < 0.0) then ArgError('InterestRate');
1509 T := 0.0; { Guess at solution }
1510 for Count := 1 to MaxIterations do
1511 begin
1512 EnT := Exp(NPeriods * T);
1513 if {LostPrecision(EnT)} ent=(ent+1) then
1514 begin
1515 Result := -Pmt / First;
1516 if Reverse then
1517 Result := Exp(-LnXP1(Result)) - 1.0;
1518 Exit;
1519 end;
1520 ET := Exp(T);
1521 ET1 := ET - 1.0;
1522 if ET1 = 0.0 then
1523 begin
1524 X := NPeriods;
1525 Y := X * (X - 1.0) / 2.0
1526 end
1527 else
1528 begin
1529 X := ET * (Exp((NPeriods - 1) * T)-1.0) / ET1;
1530 Y := (NPeriods * EnT - ET - X * ET) / ET1
1531 end;
1532 Z := Pmt * X + Last * EnT;
1533 Y := Ln(Z / -First) / ((Pmt * Y + Last * NPeriods *EnT) / Z);
1534 T := T - Y;
1535 if RelSmall(Y, T) then
1536 begin
1537 if not Reverse then T := -T;
1538 InterestRate := Exp(T)-1.0;
1539 Exit;
1540 end
1541 end;
1542 ArgError('InterestRate');
1543 end;
1545 function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
1546 PaymentTime: TPaymentTime): Extended;
1548 { If Rate = 0 then nper := -(pv + fv) / pmt
1549 else cf := pv + pmt * (1 + rate*pmtTime) / rate
1550 nper := LnXP1(-(pv + fv) / cf) / LnXP1(rate) }
1552 var
1553 PVRPP: Extended; { =PV*Rate+Payment } {"initial cash flow"}
1554 T: Extended;
1556 begin
1558 if Rate <= -1.0 then ArgError('NumberOfPeriods');
1560 {whenever both Payment and PaymentTime are given together, the PaymentTime has the effect
1561 of modifying the effective Payment by the interest accrued on the Payment}
1563 if ( PaymentTime=ptStartOfPeriod ) then
1564 Payment:=Payment*(1+Rate);
1566 {if the payment exactly matches the interest accrued periodically on the
1567 presentvalue, then an infinite number of payments are going to be
1568 required to effect a change from presentvalue to futurevalue. The
1569 following catches that specific error where payment is exactly equal,
1570 but opposite in sign to the interest on the present value. If PVRPP
1571 ("initial cash flow") is simply close to zero, the computation will
1572 be numerically unstable, but not as likely to cause an error.}
1574 PVRPP:=PresentValue*Rate+Payment;
1575 if PVRPP=0 then ArgError('NumberOfPeriods');
1577 { 6.1E-5 approx= 2**-14 }
1578 if ( EAbs(Rate)<6.1E-5 ) then
1579 Result:=-(PresentValue+FutureValue)/PVRPP
1580 else
1581 begin
1583 {starting with the initial cash flow, each compounding period cash flow
1584 should result in the current value approaching the final value. The
1585 following test combines a number of simultaneous conditions to ensure
1586 reasonableness of the cashflow before computing the NPER.}
1588 T:= -(PresentValue+FutureValue)*Rate/PVRPP;
1589 if T<=-1.0 then ArgError('NumberOfPeriods');
1590 Result := LnXP1(T) / LnXP1(Rate)
1591 end;
1592 NumberOfPeriods:=Result;
1593 end;
1595 function Payment(Rate: Extended; NPeriods: Integer; PresentValue, FutureValue:
1596 Extended; PaymentTime: TPaymentTime): Extended;
1597 var
1598 Annuity, CompoundRN: Extended;
1599 begin
1600 if Rate <= -1.0 then ArgError('Payment');
1601 Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
1602 if CompoundRN > 1.0E16 then
1603 Payment := -PresentValue * Rate / (1 + Integer(PaymentTime) * Rate)
1604 else
1605 Payment := (-PresentValue * CompoundRN - FutureValue) / Annuity
1606 end;
1608 function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
1609 PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
1610 var
1611 Junk: Extended;
1612 begin
1613 if (Rate <= -1.0) or (Period < 1) or (Period > NPeriods) then ArgError('PeriodPayment');
1614 PeriodPayment := PaymentParts(Period, NPeriods, Rate, PresentValue,
1615 FutureValue, PaymentTime, Junk);
1616 end;
1618 function PresentValue(Rate: Extended; NPeriods: Integer; Payment, FutureValue:
1619 Extended; PaymentTime: TPaymentTime): Extended;
1620 var
1621 Annuity, CompoundRN: Extended;
1622 begin
1623 if Rate <= -1.0 then ArgError('PresentValue');
1624 Annuity := Annuity2(Rate, NPeriods, PaymentTime, CompoundRN);
1625 if CompoundRN > 1.0E16 then
1626 PresentValue := -(Payment / Rate * Integer(PaymentTime) * Payment)
1627 else
1628 PresentValue := (-Payment * Annuity - FutureValue) / CompoundRN
1629 end;
1631 end.