-
Notifications
You must be signed in to change notification settings - Fork 31
/
kolmath.pas
1631 lines (1426 loc) · 46.4 KB
/
kolmath.pas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
{=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
KKKKK KKKKK OOOOOOOOO LLLLL
KKKKK KKKKK OOOOOOOOOOOOO LLLLL
KKKKK KKKKK OOOOO OOOOO LLLLL
KKKKK KKKKK OOOOO OOOOO LLLLL
KKKKKKKKKK OOOOO OOOOO LLLLL
KKKKK KKKKK OOOOO OOOOO LLLLL
KKKKK KKKKK OOOOO OOOOO LLLLL
KKKKK KKKKK OOOOOOOOOOOOO LLLLLLLLLLLLL
KKKKK KKKKK OOOOOOOOO LLLLLLLLLLLLL
Key Objects Library (C) 2000 by Kladov Vladimir.
mailto: [email protected]
Home: http://kol.nm.ru
http://xcl.cjb.net
http://xcl.nm.ru
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-}
{
This code is grabbed from standard math.pas unit,
provided by Borland Delphi. This unit is for working with
engineering (mathematical) functions. The main difference
is that err unit specially designed to handle exceptions
for KOL is used instead of SysUtils. This allows to make
size of the executable smaller for about 5K. though this
value is insignificant for project made with VCL, it can
be more than 15% of executable file size made with KOL.
}
{*******************************************************}
{ }
{ Borland Delphi Runtime Library }
{ Math Unit }
{ }
{ Copyright (C) 1996,99 Inprise Corporation }
{ }
{*******************************************************}
unit kolmath;
{ This unit contains high-performance arithmetic, trigonometric, logorithmic,
statistical and financial calculation routines which supplement the math
routines that are part of the Delphi language or System unit. }
{$N+,S-}
{$I KOLDEF.INC}
interface
uses err, kol;
const { Ranges of the IEEE floating point types, including denormals }
MinSingle = 1.5e-45;
MaxSingle = 3.4e+38;
MinDouble = 5.0e-324;
MaxDouble = 1.7e+308;
MinExtended = 3.4e-4932;
MaxExtended = 1.1e+4932;
MinComp = -9.223372036854775807e+18;
MaxComp = 9.223372036854775807e+18;
{-----------------------------------------------------------------------
References:
1) P.J. Plauger, "The Standard C Library", Prentice-Hall, 1992, Ch. 7.
2) W.J. Cody, Jr., and W. Waite, "Software Manual For the Elementary
Functions", Prentice-Hall, 1980.
3) Namir Shammas, "C/C++ Mathematical Algorithms for Scientists and Engineers",
McGraw-Hill, 1995, Ch 8.
4) H.T. Lau, "A Numerical Library in C for Scientists and Engineers",
CRC Press, 1994, Ch. 6.
5) "Pentium(tm) Processor User's Manual, Volume 3: Architecture
and Programming Manual", Intel, 1994
All angle parameters and results of trig functions are in radians.
Most of the following trig and log routines map directly to Intel 80387 FPU
floating point machine instructions. Input domains, output ranges, and
error handling are determined largely by the FPU hardware.
Routines coded in assembler favor the Pentium FPU pipeline architecture.
-----------------------------------------------------------------------}
function EAbs( D: Double ): Double;
function EMax( const Values: array of Double ): Double;
function EMin( const Values: array of Double ): Double;
function iMax( const Values: array of Integer ): Integer;
function iMin( const Values: array of Integer ): Integer;
{ Trigonometric functions }
function ArcCos(X: Extended): Extended; { IN: |X| <= 1 OUT: [0..PI] radians }
function ArcSin(X: Extended): Extended; { IN: |X| <= 1 OUT: [-PI/2..PI/2] radians }
{ ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
IN: |Y| < 2^64, |X| < 2^64, X <> 0 OUT: [-PI..PI] radians }
function ArcTan2(Y, X: Extended): Extended;
{ SinCos is 2x faster than calling Sin and Cos separately for the same angle }
procedure SinCos(Theta: Extended; var Sin, Cos: Extended) register;
function Tan(X: Extended): Extended;
function Cotan(X: Extended): Extended; { 1 / tan(X), X <> 0 }
function Hypot(X, Y: Extended): Extended; { Sqrt(X**2 + Y**2) }
{ Angle unit conversion routines }
function DegToRad(Degrees: Extended): Extended; { Radians := Degrees * PI / 180}
function RadToDeg(Radians: Extended): Extended; { Degrees := Radians * 180 / PI }
function GradToRad(Grads: Extended): Extended; { Radians := Grads * PI / 200 }
function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }
function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
{ Hyperbolic functions and inverses }
function Cosh(X: Extended): Extended;
function Sinh(X: Extended): Extended;
function Tanh(X: Extended): Extended;
function ArcCosh(X: Extended): Extended; { IN: X >= 1 }
function ArcSinh(X: Extended): Extended;
function ArcTanh(X: Extended): Extended; { IN: |X| <= 1 }
{ Logorithmic functions }
function LnXP1(X: Extended): Extended; { Ln(X + 1), accurate for X near zero }
function Log10(X: Extended): Extended; { Log base 10 of X}
function Log2(X: Extended): Extended; { Log base 2 of X }
function LogN(Base, X: Extended): Extended; { Log base N of X }
{ Exponential functions }
{ IntPower: Raise base to an integral power. Fast. }
//function IntPower(Base: Extended; Exponent: Integer): Extended register;
// -- already defined in kol.pas
{ Power: Raise base to any power.
For fractional exponents, or |exponents| > MaxInt, base must be > 0. }
function Power(Base, Exponent: Extended): Extended;
{ Miscellaneous Routines }
{ Frexp: Separates the mantissa and exponent of X. }
procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer) register;
{ Ldexp: returns X*2**P }
function Ldexp(X: Extended; P: Integer): Extended register;
{ Ceil: Smallest integer >= X, |X| < MaxInt }
function Ceil(X: Extended):Integer;
{ Floor: Largest integer <= X, |X| < MaxInt }
function Floor(X: Extended): Integer;
{ Poly: Evaluates a uniform polynomial of one variable at value X.
The coefficients are ordered in increasing powers of X:
Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
function Poly(X: Extended; const Coefficients: array of Double): Extended;
{-----------------------------------------------------------------------
Statistical functions.
Common commercial spreadsheet macro names for these statistical and
financial functions are given in the comments preceding each function.
-----------------------------------------------------------------------}
{ Mean: Arithmetic average of values. (AVG): SUM / N }
function Mean(const Data: array of Double): Extended;
{ Sum: Sum of values. (SUM) }
function Sum(const Data: array of Double): Extended register;
function SumInt(const Data: array of Integer): Integer register;
function SumOfSquares(const Data: array of Double): Extended;
procedure SumsAndSquares(const Data: array of Double;
var Sum, SumOfSquares: Extended) register;
{ MinValue: Returns the smallest signed value in the data array (MIN) }
function MinValue(const Data: array of Double): Double;
function MinIntValue(const Data: array of Integer): Integer;
function Min(A,B: Integer): Integer;
{$IFDEF _D4orHigher}
overload;
function Min(A,B: I64): I64; overload;
function Min(A,B: Single): Single; overload;
function Min(A,B: Double): Double; overload;
function Min(A,B: Extended): Extended; overload;
{$ENDIF}
{ MaxValue: Returns the largest signed value in the data array (MAX) }
function MaxValue(const Data: array of Double): Double;
function MaxIntValue(const Data: array of Integer): Integer;
function Max(A,B: Integer): Integer;
{$IFDEF _D4orHigher}
overload;
function Max(A,B: I64): I64; overload;
function Max(A,B: Single): Single; overload;
function Max(A,B: Double): Double; overload;
function Max(A,B: Extended): Extended; overload;
{$ENDIF}
{ Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
function StdDev(const Data: array of Double): Extended;
{ MeanAndStdDev calculates Mean and StdDev in one call. }
procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
{ Population Standard Deviation (STDP): Sqrt(PopnVariance).
Used in some business and financial calculations. }
function PopnStdDev(const Data: array of Double): Extended;
{ Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
function Variance(const Data: array of Double): Extended;
{ Population Variance (VAR or VARP): TotalVariance/ N }
function PopnVariance(const Data: array of Double): Extended;
{ Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
function TotalVariance(const Data: array of Double): Extended;
{ Norm: The Euclidean L2-norm. Sqrt(SumOfSquares) }
function Norm(const Data: array of Double): Extended;
{ MomentSkewKurtosis: Calculates the core factors of statistical analysis:
the first four moments plus the coefficients of skewness and kurtosis.
M1 is the Mean. M2 is the Variance.
Skew reflects symmetry of distribution: M3 / (M2**(3/2))
Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
procedure MomentSkewKurtosis(const Data: array of Double;
var M1, M2, M3, M4, Skew, Kurtosis: Extended);
{ RandG produces random numbers with Gaussian distribution about the mean.
Useful for simulating data with sampling errors. }
function RandG(Mean, StdDev: Extended): Extended;
{-----------------------------------------------------------------------
Financial functions. Standard set from Quattro Pro.
Parameter conventions:
From the point of view of A, amounts received by A are positive and
amounts disbursed by A are negative (e.g. a borrower's loan repayments
are regarded by the borrower as negative).
Interest rates are per payment period. 11% annual percentage rate on a
loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667
-----------------------------------------------------------------------}
type
TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);
{ Double Declining Balance (DDB) }
function DoubleDecliningBalance(Cost, Salvage: Extended;
Life, Period: Integer): Extended;
{ Future Value (FVAL) }
function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
Extended; PaymentTime: TPaymentTime): Extended;
{ Interest Payment (IPAYMT) }
function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
{ Interest Rate (IRATE) }
function InterestRate(NPeriods: Integer;
Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
{ Internal Rate of Return. (IRR) Needs array of cash flows. }
function InternalRateOfReturn(Guess: Extended;
const CashFlows: array of Double): Extended;
{ Number of Periods (NPER) }
function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
PaymentTime: TPaymentTime): Extended;
{ Net Present Value. (NPV) Needs array of cash flows. }
function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
PaymentTime: TPaymentTime): Extended;
{ Payment (PAYMT) }
function Payment(Rate: Extended; NPeriods: Integer;
PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
{ Period Payment (PPAYMT) }
function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
{ Present Value (PVAL) }
function PresentValue(Rate: Extended; NPeriods: Integer;
Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;
{ Straight Line depreciation (SLN) }
function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;
{ Sum-of-Years-Digits depreciation (SYD) }
function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;
{type
EInvalidArgument = class(EMathError) end;}
implementation
{$IFNDEF _D2orD3}
uses SysConst;
{$ENDIF}
function EAbs( D: Double ): Double;
begin
Result := D;
if Result < 0.0 then
Result := -Result;
end;
function EMax( const Values: array of Double ): Double;
var I: Integer;
begin
Result := Values[ 0 ];
for I := 1 to High( Values ) do
if Result < Values[ I ] then Result := Values[ I ];
end;
function EMin( const Values: array of Double ): Double;
var I: Integer;
begin
Result := Values[ 0 ];
for I := 1 to High( Values ) do
if Result > Values[ I ] then Result := Values[ I ];
end;
function iMax( const Values: array of Integer ): Integer;
var I: Integer;
begin
Result := Values[ 0 ];
for I := 1 to High( Values ) do
if Result < Values[ I ] then Result := Values[ I ];
end;
function iMin( const Values: array of Integer ): Integer;
var I: Integer;
begin
Result := Values[ 0 ];
for I := 1 to High( Values ) do
if Result > Values[ I ] then Result := Values[ I ];
end;
function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
var CompoundRN: Extended): Extended; Forward;
function Compound(R: Extended; N: Integer): Extended; Forward;
function RelSmall(X, Y: Extended): Boolean; Forward;
type
TPoly = record
Neg, Pos, DNeg, DPos: Extended
end;
const
MaxIterations = 15;
procedure ArgError(const Msg: string);
begin
raise Exception.Create(e_Math_InvalidArgument, Msg);
end;
function DegToRad(Degrees: Extended): Extended; { Radians := Degrees * PI / 180 }
begin
Result := Degrees * (PI / 180);
end;
function RadToDeg(Radians: Extended): Extended; { Degrees := Radians * 180 / PI }
begin
Result := Radians * (180 / PI);
end;
function GradToRad(Grads: Extended): Extended; { Radians := Grads * PI / 200 }
begin
Result := Grads * (PI / 200);
end;
function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI}
begin
Result := Radians * (200 / PI);
end;
function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
begin
Result := Cycles * (2 * PI);
end;
function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }
begin
Result := Radians / (2 * PI);
end;
function LnXP1(X: Extended): Extended;
{ Return ln(1 + X). Accurate for X near 0. }
asm
FLDLN2
MOV AX,WORD PTR X+8 { exponent }
FLD X
CMP AX,$3FFD { .4225 }
JB @@1
FLD1
FADD
FYL2X
JMP @@2
@@1:
FYL2XP1
@@2:
FWAIT
end;
{ Invariant: Y >= 0 & Result*X**Y = X**I. Init Y = I and Result = 1. }
{function IntPower(X: Extended; I: Integer): Extended;
var
Y: Integer;
begin
Y := Abs(I);
Result := 1.0;
while Y > 0 do begin
while not Odd(Y) do
begin
Y := Y shr 1;
X := X * X
end;
Dec(Y);
Result := Result * X
end;
if I < 0 then Result := 1.0 / Result
end;
}
(* -- already defined in kol.pas
function IntPower(Base: Extended; Exponent: Integer): Extended;
asm
mov ecx, eax
cdq
fld1 { Result := 1 }
xor eax, edx
sub eax, edx { eax := Abs(Exponent) }
jz @@3
fld Base
jmp @@2
@@1: fmul ST, ST { X := Base * Base }
@@2: shr eax,1
jnc @@1
fmul ST(1),ST { Result := Result * X }
jnz @@1
fstp st { pop X from FPU stack }
cmp ecx, 0
jge @@3
fld1
fdivrp { Result := 1 / Result }
@@3:
fwait
end;
*)
function Compound(R: Extended; N: Integer): Extended;
{ Return (1 + R)**N. }
begin
Result := IntPower(1.0 + R, N)
end;
function Annuity2(R: Extended; N: Integer; PaymentTime: TPaymentTime;
var CompoundRN: Extended): Extended;
{ Set CompoundRN to Compound(R, N),
return (1+Rate*PaymentTime)*(Compound(R,N)-1)/R;
}
begin
if R = 0.0 then
begin
CompoundRN := 1.0;
Result := N;
end
else
begin
{ 6.1E-5 approx= 2**-14 }
if EAbs(R) < 6.1E-5 then
begin
CompoundRN := Exp(N * LnXP1(R));
Result := N*(1+(N-1)*R/2);
end
else
begin
CompoundRN := Compound(R, N);
Result := (CompoundRN-1) / R
end;
if PaymentTime = ptStartOfPeriod then
Result := Result * (1 + R);
end;
end; {Annuity2}
procedure PolyX(const A: array of Double; X: Extended; var Poly: TPoly);
{ Compute A[0] + A[1]*X + ... + A[N]*X**N and X * its derivative.
Accumulate positive and negative terms separately. }
var
I: Integer;
Neg, Pos, DNeg, DPos: Extended;
begin
Neg := 0.0;
Pos := 0.0;
DNeg := 0.0;
DPos := 0.0;
for I := High(A) downto Low(A) do
begin
DNeg := X * DNeg + Neg;
Neg := Neg * X;
DPos := X * DPos + Pos;
Pos := Pos * X;
if A[I] >= 0.0 then
Pos := Pos + A[I]
else
Neg := Neg + A[I]
end;
Poly.Neg := Neg;
Poly.Pos := Pos;
Poly.DNeg := DNeg * X;
Poly.DPos := DPos * X;
end; {PolyX}
function RelSmall(X, Y: Extended): Boolean;
{ Returns True if X is small relative to Y }
const
C1: Double = 1E-15;
C2: Double = 1E-12;
begin
Result := EAbs(X) < (C1 + C2 * EAbs(Y))
end;
{ Math functions. }
function ArcCos(X: Extended): Extended;
begin
Result := ArcTan2(Sqrt(1 - X*X), X);
end;
function ArcSin(X: Extended): Extended;
begin
Result := ArcTan2(X, Sqrt(1 - X*X))
end;
function ArcTan2(Y, X: Extended): Extended;
asm
FLD Y
FLD X
FPATAN
FWAIT
end;
function Tan(X: Extended): Extended;
{ Tan := Sin(X) / Cos(X) }
asm
FLD X
FPTAN
FSTP ST(0) { FPTAN pushes 1.0 after result }
FWAIT
end;
function CoTan(X: Extended): Extended;
{ CoTan := Cos(X) / Sin(X) = 1 / Tan(X) }
asm
FLD X
FPTAN
FDIVRP
FWAIT
end;
function Hypot(X, Y: Extended): Extended;
{ formula: Sqrt(X*X + Y*Y)
implemented as: |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision
var
Temp: Extended;
begin
X := Abs(X);
Y := Abs(Y);
if X > Y then
begin
Temp := X;
X := Y;
Y := Temp;
end;
if X = 0 then
Result := Y
else // Y > X, X <> 0, so Y > 0
Result := Y * Sqrt(1 + Sqr(X/Y));
end;
}
asm
FLD Y
FABS
FLD X
FABS
FCOM
FNSTSW AX
TEST AH,$45
JNZ @@1 // if ST > ST(1) then swap
FXCH ST(1) // put larger number in ST(1)
@@1: FLDZ
FCOMP
FNSTSW AX
TEST AH,$40 // if ST = 0, return ST(1)
JZ @@2
FSTP ST // eat ST(0)
JMP @@3
@@2: FDIV ST,ST(1) // ST := ST / ST(1)
FMUL ST,ST // ST := ST * ST
FLD1
FADD // ST := ST + 1
FSQRT // ST := Sqrt(ST)
FMUL // ST(1) := ST * ST(1); Pop ST
@@3: FWAIT
end;
procedure SinCos(Theta: Extended; var Sin, Cos: Extended);
asm
FLD Theta
FSINCOS
FSTP tbyte ptr [edx] // Cos
FSTP tbyte ptr [eax] // Sin
FWAIT
end;
{ Extract exponent and mantissa from X }
procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer);
{ Mantissa ptr in EAX, Exponent ptr in EDX }
asm
FLD X
PUSH EAX
MOV dword ptr [edx], 0 { if X = 0, return 0 }
FTST
FSTSW AX
FWAIT
SAHF
JZ @@Done
FXTRACT // ST(1) = exponent, (pushed) ST = fraction
FXCH
// The FXTRACT instruction normalizes the fraction 1 bit higher than
// wanted for the definition of frexp() so we need to tweak the result
// by scaling the fraction down and incrementing the exponent.
FISTP dword ptr [edx]
FLD1
FCHS
FXCH
FSCALE // scale fraction
INC dword ptr [edx] // exponent biased to match
FSTP ST(1) // discard -1, leave fraction as TOS
@@Done:
POP EAX
FSTP tbyte ptr [eax]
FWAIT
end;
function Ldexp(X: Extended; P: Integer): Extended;
{ Result := X * (2^P) }
asm
PUSH EAX
FILD dword ptr [ESP]
FLD X
FSCALE
POP EAX
FSTP ST(1)
FWAIT
end;
function Ceil(X: Extended): Integer;
begin
Result := Integer(Trunc(X));
if Frac(X) > 0 then
Inc(Result);
end;
function Floor(X: Extended): Integer;
begin
Result := Integer(Trunc(X));
if Frac(X) < 0 then
Dec(Result);
end;
{ Conversion of bases: Log.b(X) = Log.a(X) / Log.a(b) }
function Log10(X: Extended): Extended;
{ Log.10(X) := Log.2(X) * Log.10(2) }
asm
FLDLG2 { Log base ten of 2 }
FLD X
FYL2X
FWAIT
end;
function Log2(X: Extended): Extended;
asm
FLD1
FLD X
FYL2X
FWAIT
end;
function LogN(Base, X: Extended): Extended;
{ Log.N(X) := Log.2(X) / Log.2(N) }
asm
FLD1
FLD X
FYL2X
FLD1
FLD Base
FYL2X
FDIV
FWAIT
end;
function Poly(X: Extended; const Coefficients: array of Double): Extended;
{ Horner's method }
var
I: Integer;
begin
Result := Coefficients[High(Coefficients)];
for I := High(Coefficients)-1 downto Low(Coefficients) do
Result := Result * X + Coefficients[I];
end;
function Power(Base, Exponent: Extended): Extended;
begin
if Exponent = 0.0 then
Result := 1.0 { n**0 = 1 }
else if (Base = 0.0) and (Exponent > 0.0) then
Result := 0.0 { 0**n = 0, n > 0 }
else if (Frac(Exponent) = 0.0) and (EAbs(Exponent) <= MaxInt) then
Result := IntPower(Base, Integer(Trunc(Exponent)))
else
Result := Exp(Exponent * Ln(Base))
end;
{ Hyperbolic functions }
function CoshSinh(X: Extended; Factor: Double): Extended;
begin
Result := Exp(X) / 2;
Result := Result + Factor / Result;
end;
function Cosh(X: Extended): Extended;
begin
Result := CoshSinh(X, 0.25)
end;
function Sinh(X: Extended): Extended;
begin
Result := CoshSinh(X, -0.25)
end;
const
MaxTanhDomain = 5678.22249441322; // Ln(MaxExtended)/2
function Tanh(X: Extended): Extended;
begin
if X > MaxTanhDomain then
Result := 1.0
else if X < -MaxTanhDomain then
Result := -1.0
else
begin
Result := Exp(X);
Result := Result * Result;
Result := (Result - 1.0) / (Result + 1.0)
end;
end;
function ArcCosh(X: Extended): Extended;
begin
if X <= 1.0 then
Result := 0.0
else if X > 1.0e10 then
Result := Ln(2) + Ln(X)
else
Result := Ln(X + Sqrt((X - 1.0) * (X + 1.0)));
end;
function ArcSinh(X: Extended): Extended;
var
Neg: Boolean;
begin
if X = 0 then
Result := 0
else
begin
Neg := (X < 0);
X := EAbs(X);
if X > 1.0e10 then
Result := Ln(2) + Ln(X)
else
begin
Result := X*X;
Result := LnXP1(X + Result / (1 + Sqrt(1 + Result)));
end;
if Neg then Result := -Result;
end;
end;
function ArcTanh(X: Extended): Extended;
var
Neg: Boolean;
begin
if X = 0 then
Result := 0
else
begin
Neg := (X < 0);
X := EAbs(X);
if X >= 1 then
Result := MaxExtended
else
Result := 0.5 * LnXP1((2.0 * X) / (1.0 - X));
if Neg then Result := -Result;
end;
end;
{ Statistical functions }
function Mean(const Data: array of Double): Extended;
begin
Result := SUM(Data) / (High(Data) - Low(Data) + 1)
end;
function MinValue(const Data: array of Double): Double;
var
I: Integer;
begin
Result := Data[Low(Data)];
for I := Low(Data) + 1 to High(Data) do
if Result > Data[I] then
Result := Data[I];
end;
function MinIntValue(const Data: array of Integer): Integer;
var
I: Integer;
begin
Result := Data[Low(Data)];
for I := Low(Data) + 1 to High(Data) do
if Result > Data[I] then
Result := Data[I];
end;
function Min(A,B: Integer): Integer;
begin
if A < B then
Result := A
else
Result := B;
end;
{$IFDEF _D4orHigher}
function Min(A,B: I64): I64;
begin
if Cmp64( A, B ) < 0 then
Result := A
else
Result := B;
end;
function Min(A,B: Single): Single;
begin
if A < B then
Result := A
else
Result := B;
end;
function Min(A,B: Double): Double;
begin
if A < B then
Result := A
else
Result := B;
end;
function Min(A,B: Extended): Extended;
begin
if A < B then
Result := A
else
Result := B;
end;
{$ENDIF}
function MaxValue(const Data: array of Double): Double;
var
I: Integer;
begin
Result := Data[Low(Data)];
for I := Low(Data) + 1 to High(Data) do
if Result < Data[I] then
Result := Data[I];
end;
function MaxIntValue(const Data: array of Integer): Integer;
var
I: Integer;
begin
Result := Data[Low(Data)];
for I := Low(Data) + 1 to High(Data) do
if Result < Data[I] then
Result := Data[I];
end;
function Max(A,B: Integer): Integer;
begin
if A > B then
Result := A
else
Result := B;
end;
{$IFDEF _D4orHigher}
function Max(A,B: I64): I64;
begin
if Cmp64( A, B ) > 0 then
Result := A
else
Result := B;
end;
function Max(A,B: Single): Single;
begin
if A > B then
Result := A
else
Result := B;
end;
function Max(A,B: Double): Double;
begin
if A > B then
Result := A
else
Result := B;
end;
function Max(A,B: Extended): Extended;
begin
if A > B then
Result := A
else
Result := B;
end;
{$ENDIF}
procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);
var
S: Extended;
N,I: Integer;
begin
N := High(Data)- Low(Data) + 1;
if N = 1 then
begin
Mean := Data[0];
StdDev := Data[0];
Exit;
end;
Mean := Sum(Data) / N;
S := 0; // sum differences from the mean, for greater accuracy
for I := Low(Data) to High(Data) do
S := S + Sqr(Mean - Data[I]);
StdDev := Sqrt(S / (N - 1));
end;
procedure MomentSkewKurtosis(const Data: array of Double;
var M1, M2, M3, M4, Skew, Kurtosis: Extended);
var
Sum, SumSquares, SumCubes, SumQuads, OverN, Accum, M1Sqr, S2N, S3N: Extended;
I: Integer;
begin
OverN := 1 / (High(Data) - Low(Data) + 1);
Sum := 0;
SumSquares := 0;
SumCubes := 0;
SumQuads := 0;
for I := Low(Data) to High(Data) do
begin
Sum := Sum + Data[I];
Accum := Sqr(Data[I]);
SumSquares := SumSquares + Accum;
Accum := Accum*Data[I];
SumCubes := SumCubes + Accum;
SumQuads := SumQuads + Accum*Data[I];
end;
M1 := Sum * OverN;
M1Sqr := Sqr(M1);
S2N := SumSquares * OverN;
S3N := SumCubes * OverN;
M2 := S2N - M1Sqr;
M3 := S3N - (M1 * 3 * S2N) + 2*M1Sqr*M1;
M4 := (SumQuads * OverN) - (M1 * 4 * S3N) + (M1Sqr*6*S2N - 3*Sqr(M1Sqr));