math.pas 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887
  1. {
  2. This file is part of the Pas2JS run time library.
  3. Copyright (c) 2017 by Mattias Gaertner
  4. See the file COPYING.FPC, included in this distribution,
  5. for details about the copyright.
  6. This program is distributed in the hope that it will be useful,
  7. but WITHOUT ANY WARRANTY; without even the implied warranty of
  8. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  9. **********************************************************************}
  10. unit Math;
  11. {$mode objfpc}
  12. interface
  13. const
  14. MinInteger = -$fffffffffffff-1;
  15. MaxInteger = $fffffffffffff;
  16. MinDouble = 5.0e-324;
  17. MaxDouble = 1.7e+308;
  18. const
  19. NaN: Double; external name 'NaN';
  20. Infinity: Double; external name 'Infinity';
  21. NegInfinity: Double; external name '-Infinity';
  22. type
  23. float = double;
  24. //EInvalidArgument = class(EMathError);
  25. function Min(const a, b: Double): Double; varargs; external name 'Math.min'; overload;
  26. function Max(const a, b: Double): Double; varargs; external name 'Math.max'; overload;
  27. function Min(const a, b: NativeLargeUInt): NativeLargeUInt; varargs; external name 'Math.min'; overload;
  28. function Max(const a, b: NativeLargeUInt): NativeLargeUInt; varargs; external name 'Math.max'; overload;
  29. function Min(const a, b: NativeLargeInt): NativeLargeInt; varargs; external name 'Math.min'; overload;
  30. function Max(const a, b: NativeLargeInt): NativeLargeInt; varargs; external name 'Math.max'; overload;
  31. function Min(const a, b: Integer): Integer; varargs; external name 'Math.min'; overload;
  32. function Max(const a, b: Integer): Integer; varargs; external name 'Math.max'; overload;
  33. function InRange(const AValue, AMin, AMax: Integer): Boolean; assembler; overload;
  34. function InRange(const AValue, AMin, AMax: Double): Boolean; assembler; overload;
  35. function EnsureRange(const AValue, AMin, AMax: Integer): Integer; assembler; overload;
  36. function EnsureRange(const AValue, AMin, AMax: Double): Double; assembler; overload;
  37. type
  38. TRoundToRange = -37..37;
  39. function RoundTo(const AValue: Double; const Digits: TRoundToRange): Double;
  40. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange = -2): Double;
  41. function randg(mean,stddev : float) : float;
  42. function RandomRange(const aFrom, aTo: Integer): Integer;
  43. function RandomRange(const aFrom, aTo: NativeLargeInt): NativeLargeInt;
  44. Type
  45. TValueSign = -1..1;
  46. const
  47. NegativeValue = -1;
  48. ZeroValue = 0;
  49. PositiveValue = 1;
  50. function Sign(const AValue: NativeInt): TValueSign; external name 'Math.sign'; overload;
  51. function Sign(const AValue: Double): TValueSign; external name 'Math.sign'; overload;
  52. function IsZero(const d: Double; Epsilon: Double): Boolean; overload;
  53. function IsZero(const d: Double): Boolean; overload;
  54. function IsNaN(const v: JSValue): boolean; external name {$IFDEF ECMAScript5}'isNaN'{$ELSE}'Number.isNaN'{$ENDIF}; overload;
  55. function IsFinite(const d: JSValue): Boolean; external name 'isFinite'; overload;// false if NaN, positive or negative infinity
  56. function IsInfinite(const d: JSValue): Boolean; assembler; overload; // negative or positive infinity
  57. {$IFDEF ECMAScript6}
  58. function IsInteger(const d: JSValue): Boolean; external name 'Number.isInteger'; // any integer representable by a double
  59. function IsSafeInteger(const d: JSValue): Boolean; external name 'Number.isSafeInteger'; // an integer between MinInteger and MaxInteger, inclusive
  60. {$ENDIF}
  61. function SameValue(const A, B: Double; Epsilon: Double = 0.0): Boolean; overload;
  62. // function Abs is in System.pas
  63. function ArcCos(const A : Double): Double; external name 'Math.acos';
  64. function ArcCosH(const A: Double): Double; external name 'Math.acosh'; // not on IE
  65. function ArcSin(const A : Double): Double; external name 'Math.asin';
  66. function ArcSinH(const A : Double): Double; external name 'Math.asinh'; // not on IE
  67. function ArcTanH(const A: Double): Double; external name 'Math.atanh'; // not on IE
  68. function CosH(const A: Double): Double; external name 'Math.cosh'; // not on IE
  69. procedure SinCos(const A: Double; var B: Double; var C: Double);
  70. function ExpM1(const A: Double): Double; external name 'Math.expm1'; // not on IE
  71. function FRound(const A: Double): Double; overload; external name 'Math.fround'; // not on IE
  72. function FTrunc(const A: Double): double; overload; external name 'Math.trunc'; // not on IE
  73. function Hypot(const A : Double): Double; varargs; external name 'Math.hypot'; // not on IE
  74. function IMul(const A, B: Integer): Integer; external name 'Math.imul'; // not on IE
  75. function Log10(const A: Double): Double; external name 'Math.log10';
  76. function Log1p(const A: Double): Double; external name 'Math.log1p'; // not on IE
  77. function Log2(const A: Double): Double; external name 'Math.log2'; // not on IE
  78. function LogN(const Base, N: Double): Double;
  79. function Power(const Base, Exponent: Double): Double; external name 'Math.pow';
  80. // ln, round, sqrt, trunc, cos, sin, arctan, round, exp are in unit system
  81. function Ceil(const A: Double): Integer;
  82. function Floor(const A: Double): Integer;
  83. function Ceil64(const A: Double): NativeLargeInt;
  84. function Floor64(const A: Double): NativeLargeInt;
  85. function ldexp(x : double;const p : Integer) : double;
  86. procedure Frexp(X: double; out Mantissa: double; out Exponent: integer);
  87. function lnxp1(x : double) : double;
  88. function IntPower(base : float;const exponent : Integer) : double;
  89. procedure DivMod(Dividend: LongInt; Divisor: Word; out Result, Remainder: Word);
  90. procedure DivMod(Dividend: LongInt; Divisor: Word; out Result, Remainder: SmallInt);
  91. procedure DivMod(Dividend: DWord; Divisor: DWord; out Result, Remainder: DWord);
  92. procedure DivMod(Dividend: LongInt; Divisor: LongInt; out Result, Remainder: LongInt);
  93. { Angle conversion }
  94. function DegToRad(deg : double) : double;
  95. function RadToDeg(rad : double) : double;
  96. function GradToRad(grad : double) : double;
  97. function RadToGrad(rad : double) : double;
  98. function DegToGrad(deg : double) : double;
  99. function GradToDeg(grad : double) : double;
  100. { one cycle are 2*Pi rad }
  101. function CycleToRad(cycle : double) : double;
  102. function RadToCycle(rad : double) : double;
  103. Function DegNormalize(deg : double) : double;
  104. function Norm(const data : array of double) : double;
  105. // Statistical functions
  106. function Mean(const data : array of double) : double;
  107. function Sum(const data : array of double) : double;
  108. procedure SumsAndSquares(const data : Array of Double; out Sum,SumOfSquares : double);
  109. function StdDev(const data : array of Double) : float;
  110. procedure MeanAndStdDev(const data : array of Double; out Mean,StdDev : double);
  111. function Variance(const data : array of Double) : double;
  112. function TotalVariance(const data : array of Double) : double;
  113. function PopNStdDev(const data : array of Double) : double;
  114. function PopNVariance(const data : array of Double) : double;
  115. procedure MomentSkewKurtosis(const data : array of Double; out m1,m2,m3,m4,skew,kurtosis : double);
  116. // Financial functions
  117. Type
  118. TPaymentTime = (ptEndOfPeriod,ptStartOfPeriod);
  119. function FutureValue(ARate: double; NPeriods: Integer;
  120. APayment, APresentValue: double; APaymentTime: TPaymentTime): double;
  121. function InterestRate(NPeriods: Integer; APayment, APresentValue, AFutureValue: double;
  122. APaymentTime: TPaymentTime): double;
  123. function NumberOfPeriods(ARate, APayment, APresentValue, AFutureValue: double;
  124. APaymentTime: TPaymentTime): double;
  125. function Payment(ARate: double; NPeriods: Integer;
  126. APresentValue, AFutureValue: double; APaymentTime: TPaymentTime): double;
  127. function PresentValue(ARate: double; NPeriods: Integer;
  128. APayment, AFutureValue: double; APaymentTime: TPaymentTime): double;
  129. // Miscellaneous
  130. function IfThen(val:boolean;const ifTrue:integer; const ifFalse:integer= 0) :integer; overload;
  131. function IfThen(val:boolean;const ifTrue:double ; const ifFalse:double =0.0):double; overload;
  132. Type
  133. TValueRelationship = -1..1;
  134. const
  135. EqualsValue = 0;
  136. LessThanValue = Low(TValueRelationship);
  137. GreaterThanValue = High(TValueRelationship);
  138. function CompareValue ( const A, B : Integer) : TValueRelationship;
  139. function CompareValue ( const A, B : NativeLargeInt) : TValueRelationship;
  140. function CompareValue ( const A, B : NativeLargeUInt) : TValueRelationship;
  141. function CompareValue ( const A, B : Double; delta : Double = 0.0) : TValueRelationship;
  142. implementation
  143. function InRange(const AValue, AMin, AMax: Integer): Boolean; assembler;
  144. asm
  145. return (AValue >= AMin) && (AValue <= AMax);
  146. end;
  147. function InRange(const AValue, AMin, AMax: Double): Boolean; assembler;
  148. asm
  149. return (AValue >= AMin) && (AValue <= AMax);
  150. end;
  151. function EnsureRange(const AValue, AMin, AMax: Integer): Integer; assembler;
  152. asm
  153. if (AValue<AMin){ return AMin;
  154. } else if (AValue>AMax){ return AMax;
  155. } else return AValue;
  156. end;
  157. function EnsureRange(const AValue, AMin, AMax: Double): Double; assembler;
  158. asm
  159. if (AValue<AMin){ return AMin;
  160. } else if (AValue>AMax){ return AMax;
  161. } else return AValue;
  162. end;
  163. function RoundTo(const AValue: Double; const Digits: TRoundToRange): Double;
  164. var
  165. RV : Double;
  166. begin
  167. RV:=IntPower(10,Digits);
  168. Result:=Round(AValue/RV)*RV;
  169. end;
  170. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange): Double;
  171. var
  172. RV : Double;
  173. begin
  174. RV := IntPower(10, -Digits);
  175. if AValue < 0 then
  176. Result := Int((AValue*RV) - 0.5)/RV
  177. else
  178. Result := Int((AValue*RV) + 0.5)/RV;
  179. end;
  180. function randg(mean,stddev : float) : float;
  181. Var
  182. U1,S2 : Float;
  183. begin
  184. repeat
  185. u1:= 2*random-1;
  186. S2:=Sqr(U1)+sqr(2*random-1);
  187. until s2<1;
  188. Result:=Sqrt(-2*ln(S2)/S2)*u1*stddev+Mean;
  189. end;
  190. function RandomRange(const aFrom, aTo: Integer): Integer;
  191. begin
  192. Result:=Random(Abs(aFrom-aTo))+Min(aTo,AFrom);
  193. end;
  194. function RandomRange(const aFrom, aTo: NativeLargeInt): NativeLargeInt;
  195. Var
  196. m : NativeLargeInt;
  197. begin
  198. if aFrom<aTo then
  199. M:=aFrom
  200. else
  201. M:=aTo;
  202. Result:=Random(Abs(aFrom-aTo))+M;
  203. end;
  204. const
  205. DZeroResolution = 1E-12;
  206. function IsZero(const d: Double; Epsilon: Double): Boolean;
  207. begin
  208. if Epsilon=0 then
  209. Epsilon:=DZeroResolution;
  210. Result:=Abs(d)<=Epsilon;
  211. end;
  212. function IsZero(const d: Double): Boolean;
  213. begin
  214. Result:=Abs(d)<=DZeroResolution;
  215. end;
  216. function IsInfinite(const d: JSValue): Boolean; assembler;
  217. asm
  218. return (d==Infinity) || (d==-Infinity);
  219. end;
  220. function SameValue(const A, B: Double; Epsilon: Double): Boolean;
  221. begin
  222. if (Epsilon=0.0) then
  223. Epsilon:=Max(Min(Abs(A),Abs(B))*DZeroResolution,DZeroResolution);
  224. if (A>B) then
  225. Result:=((A-B)<=Epsilon)
  226. else
  227. Result:=((B-A)<=Epsilon);
  228. end;
  229. function JSFloor(const A: Double): Double; external name 'Math.floor';
  230. function JSCeil(const A: Double): Double; external name 'Math.ceil';
  231. function Ceil(const A: Double): Integer;
  232. begin
  233. // TODO: add Range check ?
  234. Result:=trunc(JSCeil(a));
  235. end;
  236. function Floor(const A: Double): Integer;
  237. begin
  238. // TODO: add Range check ?
  239. Result:=trunc(JSFloor(a));
  240. end;
  241. function Ceil64(const A: Double): NativeLargeInt;
  242. begin
  243. Result:=trunc(JSCeil(a));
  244. end;
  245. function Floor64(const A: Double): NativeLargeInt;
  246. begin
  247. Result:=trunc(JSCeil(a));
  248. end;
  249. procedure Frexp(X: double; out Mantissa: double; out Exponent: integer);
  250. begin
  251. Exponent:=0;
  252. if (X<>0) then
  253. if (abs(X)<0.5) then
  254. repeat
  255. X:=X*2;
  256. Dec(Exponent);
  257. until (abs(X)>=0.5)
  258. else
  259. while (abs(X)>=1) do
  260. begin
  261. X:=X/2;
  262. Inc(Exponent);
  263. end;
  264. Mantissa:=X;
  265. end;
  266. procedure SinCos(const A: Double; var B: Double; var C: Double);
  267. begin
  268. B:=Sin(A);
  269. C:=Cos(A);
  270. end;
  271. function LogN(const Base,N: Double): Double;
  272. begin
  273. Result:=Ln(N)/Ln(Base);
  274. end;
  275. function lnxp1(x: double): double;
  276. var
  277. y: float;
  278. begin
  279. if (x>=4.0) then
  280. result:=ln(1.0+x)
  281. else
  282. begin
  283. y:=1.0+x;
  284. if (y=1.0) then
  285. result:=x
  286. else
  287. begin
  288. result:=ln(y); { lnxp1(-1) = ln(0) = -Inf }
  289. if y>0.0 then
  290. result:=result+(x-(y-1.0))/y;
  291. end;
  292. end;
  293. end;
  294. function ldexp(x : double;const p : Integer) : double;
  295. begin
  296. result:=x*intpower(2.0,p);
  297. end;
  298. function IntPower(base: float; const exponent: Integer): double;
  299. var
  300. i : longint;
  301. begin
  302. if (base = 0.0) and (exponent = 0) then
  303. result:=1
  304. else
  305. begin
  306. i:=abs(exponent);
  307. Result:=1.0;
  308. while i>0 do
  309. begin
  310. while (i and 1)=0 do
  311. begin
  312. i:=i shr 1;
  313. base:=sqr(base);
  314. end;
  315. i:=i-1;
  316. Result:=Result*base;
  317. end;
  318. if exponent<0 then
  319. Result:=1.0/Result;
  320. end;
  321. end;
  322. procedure DivMod(Dividend: LongInt; Divisor: Word; out Result, Remainder: Word);
  323. begin
  324. if Dividend < 0 then
  325. begin
  326. Dividend:=-Dividend;
  327. Result:=-(Dividend Div Divisor);
  328. Remainder:=-(Dividend+(Result*Divisor));
  329. end
  330. else
  331. begin
  332. Result:=Dividend Div Divisor;
  333. Remainder:=Dividend-(Result*Divisor);
  334. end;
  335. end;
  336. procedure DivMod(Dividend: LongInt; Divisor: Word; out Result,
  337. Remainder: SmallInt);
  338. begin
  339. if Dividend < 0 then
  340. begin
  341. Dividend:=-Dividend;
  342. Result:=-(Dividend Div Divisor);
  343. Remainder:=-(Dividend+(Result*Divisor));
  344. end
  345. else
  346. begin
  347. Result:=Dividend Div Divisor;
  348. Remainder:=Dividend-(Result*Divisor);
  349. end;
  350. end;
  351. procedure DivMod(Dividend: DWord; Divisor: DWord; out Result, Remainder: DWord);
  352. begin
  353. Result:=Dividend Div Divisor;
  354. Remainder:=Dividend-(Result*Divisor);
  355. end;
  356. procedure DivMod(Dividend: LongInt; Divisor: LongInt; out Result,
  357. Remainder: LongInt);
  358. begin
  359. if Dividend < 0 then
  360. begin
  361. Dividend:=-Dividend;
  362. Result:=-(Dividend Div Divisor);
  363. Remainder:=-(Dividend+(Result*Divisor));
  364. end
  365. else
  366. begin
  367. Result:=Dividend Div Divisor;
  368. Remainder:=Dividend-(Result*Divisor);
  369. end;
  370. end;
  371. { ---------------------------------------------------------------------
  372. Angle conversion
  373. ---------------------------------------------------------------------}
  374. function DegToRad(deg: double): double;
  375. begin
  376. Result:=deg*(pi/180.0);
  377. end;
  378. function RadToDeg(rad: double): double;
  379. begin
  380. Result:=rad*(180.0/pi);
  381. end;
  382. function GradToRad(grad: double): double;
  383. begin
  384. Result:=grad*(pi/200.0);
  385. end;
  386. function RadToGrad(rad: double): double;
  387. begin
  388. Result:=rad*(200.0/pi);
  389. end;
  390. function DegToGrad(deg: double): double;
  391. begin
  392. Result:=deg*(200.0/180.0);
  393. end;
  394. function GradToDeg(grad: double): double;
  395. begin
  396. Result:=grad*(180.0/200.0);
  397. end;
  398. function CycleToRad(cycle: double): double;
  399. begin
  400. Result:=(2*pi)*cycle;
  401. end;
  402. function RadToCycle(rad: double): double;
  403. begin
  404. Result:=rad*(1/(2*pi));
  405. end;
  406. function DegNormalize(deg: double): double;
  407. begin
  408. Result:=Deg-Int(Deg/360)*360;
  409. If (Result<0) then Result:=Result+360;
  410. end;
  411. function sumofsquares(const data : array of double) : double;
  412. var
  413. i,N : longint;
  414. begin
  415. N:=Length(Data);
  416. Result:=0.0;
  417. for i:=0 to N-1 do
  418. Result:=Result+sqr(data[i]);
  419. end;
  420. function Norm(const data: array of double): double;
  421. begin
  422. Result:=sqrt(sumofsquares(data));
  423. end;
  424. { ---------------------------------------------------------------------
  425. Statistical functions
  426. ---------------------------------------------------------------------}
  427. function Sum(const data: array of double): double;
  428. var
  429. i,N : longint;
  430. begin
  431. N:=Length(Data);
  432. Result:=0.0;
  433. for i:=0 to N-1 do
  434. Result:=Result+data[i];
  435. end;
  436. function Mean(const data: array of double): double;
  437. Var
  438. N : integer;
  439. begin
  440. N:=Length(Data);
  441. if N=0 then
  442. Result:=0
  443. else
  444. Result:=Sum(Data)/N;
  445. end;
  446. procedure SumsAndSquares(const data: array of Double; out Sum,
  447. SumOfSquares: double);
  448. var
  449. i,n : Integer;
  450. t,s,ss: double;
  451. begin
  452. n:=length(Data);
  453. ss:=0.0; // Use local vars, var is very inefficient in js
  454. s:=0.0;
  455. for i:=0 to N-1 do
  456. begin
  457. t:=data[i];
  458. ss:=ss+sqr(t);
  459. s:=s+t;
  460. end;
  461. Sum:=s;
  462. SumOfSquares:=ss;
  463. end;
  464. function StdDev(const data: array of Double): float;
  465. begin
  466. Result:=Sqrt(Variance(Data));
  467. end;
  468. function Variance(const data: array of Double): double;
  469. var
  470. n : integer;
  471. begin
  472. N:=Length(Data);
  473. If N=1 then
  474. Result:=0
  475. else
  476. Result:=TotalVariance(Data)/(N-1);
  477. end;
  478. function TotalVariance(const data: array of Double): double;
  479. var
  480. S,SS : Float;
  481. N : integer;
  482. begin
  483. N:=Length(Data);
  484. If Length(Data)=1 then
  485. Result:=0
  486. else
  487. begin
  488. SumsAndSquares(Data,S,SS);
  489. Result := SS-Sqr(S)/N;
  490. end;
  491. end;
  492. procedure MeanAndStdDev(const data: array of Double; out Mean, StdDev: double);
  493. Var
  494. I,N : longint;
  495. M,S : Double;
  496. begin
  497. N:=Length(Data);
  498. M:=0;
  499. S:=0;
  500. For I:=0 to N-1 do
  501. begin
  502. M:=M+Data[i];
  503. S:=S+Sqr(Data[i]);
  504. end;
  505. M:=M/N;
  506. S:=(S-N*Sqr(M));
  507. If N>1 then
  508. S:=Sqrt(S/(N-1))
  509. else
  510. S:=0;
  511. Mean:=M;
  512. StdDev:=S;
  513. end;
  514. function PopNStdDev(const data : array of Double) : double;
  515. begin
  516. Result:=Sqrt(PopnVariance(Data));
  517. end;
  518. function PopNVariance(const data : array of Double) : double;
  519. Var
  520. N : integer;
  521. begin
  522. N:=Length(Data);
  523. if N=0 then
  524. Result:=0
  525. else
  526. Result:=TotalVariance(Data)/N;
  527. end;
  528. procedure MomentSkewKurtosis(const data: array of Double; out m1, m2, m3, m4, skew, kurtosis: double);
  529. var
  530. i,N: integer;
  531. deviation, deviation2: double;
  532. reciprocalN: float;
  533. // Use local vars for all calculations, var is very slow
  534. lm1, lm2, lm3, lm4, lskew, lkurtosis: double;
  535. begin
  536. N:=length(Data);
  537. lm1 := 0;
  538. reciprocalN := 1/N;
  539. for i := 0 to N-1 do
  540. lm1 := lm1 + data[i];
  541. lm1 := reciprocalN * lm1;
  542. lm2 := 0;
  543. lm3 := 0;
  544. lm4 := 0;
  545. for i := 0 to N-1 do
  546. begin
  547. deviation := (data[i]-lm1);
  548. deviation2 := deviation * deviation;
  549. lm2 := lm2 + deviation2;
  550. lm3 := lm3 + deviation2 * deviation;
  551. lm4 := lm4 + deviation2 * deviation2;
  552. end;
  553. lm2 := reciprocalN * lm2;
  554. lm3 := reciprocalN * lm3;
  555. lm4 := reciprocalN * lm4;
  556. lskew := lm3 / (sqrt(lm2)*lm2);
  557. lkurtosis := lm4 / (lm2 * lm2);
  558. m1:=lm1;
  559. m2:=lm2;
  560. m3:=lm3;
  561. m4:=lm4;
  562. skew:=lskew;
  563. kurtosis:=lkurtosis;
  564. end;
  565. { ---------------------------------------------------------------------
  566. Financial functions
  567. ---------------------------------------------------------------------}
  568. function FutureValue(ARate: double; NPeriods: Integer;
  569. APayment, APresentValue: double; APaymentTime: TPaymentTime): double;
  570. var
  571. q, qn, factor: double;
  572. begin
  573. if ARate = 0 then
  574. Result := -APresentValue - APayment * NPeriods
  575. else begin
  576. q := 1.0 + ARate;
  577. qn := power(q, NPeriods);
  578. factor := (qn - 1) / (q - 1);
  579. if APaymentTime = ptStartOfPeriod then
  580. factor := factor * q;
  581. Result := -(APresentValue * qn + APayment*factor);
  582. end;
  583. end;
  584. function InterestRate(NPeriods: Integer; APayment, APresentValue, AFutureValue: double;
  585. APaymentTime: TPaymentTime): double;
  586. { The interest rate cannot be calculated analytically. We solve the equation
  587. numerically by means of the Newton method:
  588. - guess value for the interest reate
  589. - calculate at which interest rate the tangent of the curve fv(rate)
  590. (straight line!) has the requested future vale.
  591. - use this rate for the next iteration. }
  592. const
  593. DELTA = 0.001;
  594. EPS = 1E-9; // required precision of interest rate (after typ. 6 iterations)
  595. MAXIT = 20; // max iteration count to protect agains non-convergence
  596. var
  597. r1, r2, dr: double;
  598. fv1, fv2: double;
  599. iteration: Integer;
  600. begin
  601. iteration := 0;
  602. r1 := 0.05; // inital guess
  603. repeat
  604. r2 := r1 + DELTA;
  605. fv1 := FutureValue(r1, NPeriods, APayment, APresentValue, APaymentTime);
  606. fv2 := FutureValue(r2, NPeriods, APayment, APresentValue, APaymentTime);
  607. dr := (AFutureValue - fv1) / (fv2 - fv1) * delta; // tangent at fv(r)
  608. r1 := r1 + dr; // next guess
  609. inc(iteration);
  610. until (abs(dr) < EPS) or (iteration >= MAXIT);
  611. Result := r1;
  612. end;
  613. function NumberOfPeriods(ARate, APayment, APresentValue, AFutureValue: double;
  614. APaymentTime: TPaymentTime): double;
  615. { Solve the cash flow equation (1) for q^n and take the logarithm }
  616. var
  617. q, x1, x2: double;
  618. begin
  619. if ARate = 0 then
  620. Result := -(APresentValue + AFutureValue) / APayment
  621. else begin
  622. q := 1.0 + ARate;
  623. if APaymentTime = ptStartOfPeriod then
  624. APayment := APayment * q;
  625. x1 := APayment - AFutureValue * ARate;
  626. x2 := APayment + APresentValue * ARate;
  627. if (x2 = 0) // we have to divide by x2
  628. or (sign(x1) * sign(x2) < 0) // the argument of the log is negative
  629. then
  630. Result := Infinity
  631. else begin
  632. Result := ln(x1/x2) / ln(q);
  633. end;
  634. end;
  635. end;
  636. function Payment(ARate: double; NPeriods: Integer;
  637. APresentValue, AFutureValue: double; APaymentTime: TPaymentTime): double;
  638. var
  639. q, qn, factor: double;
  640. begin
  641. if ARate = 0 then
  642. Result := -(AFutureValue + APresentValue) / NPeriods
  643. else begin
  644. q := 1.0 + ARate;
  645. qn := power(q, NPeriods);
  646. factor := (qn - 1) / (q - 1);
  647. if APaymentTime = ptStartOfPeriod then
  648. factor := factor * q;
  649. Result := -(AFutureValue + APresentValue * qn) / factor;
  650. end;
  651. end;
  652. function PresentValue(ARate: double; NPeriods: Integer;
  653. APayment, AFutureValue: double; APaymentTime: TPaymentTime): double;
  654. var
  655. q, qn, factor: double;
  656. begin
  657. if ARate = 0.0 then
  658. Result := -AFutureValue - APayment * NPeriods
  659. else begin
  660. q := 1.0 + ARate;
  661. qn := power(q, NPeriods);
  662. factor := (qn - 1) / (q - 1);
  663. if APaymentTime = ptStartOfPeriod then
  664. factor := factor * q;
  665. Result := -(AFutureValue + APayment*factor) / qn;
  666. end;
  667. end;
  668. { ---------------------------------------------------------------------
  669. Miscellaneous
  670. ---------------------------------------------------------------------}
  671. function IfThen(val: boolean; const ifTrue: integer; const ifFalse: integer): integer;
  672. begin
  673. if val then result:=iftrue else result:=iffalse;
  674. end;
  675. function IfThen(val: boolean; const ifTrue: double; const ifFalse: double): double;
  676. begin
  677. if val then result:=iftrue else result:=iffalse;
  678. end;
  679. function CompareValue(const A, B : Integer): TValueRelationship;
  680. begin
  681. result:=GreaterThanValue;
  682. if a=b then
  683. result:=EqualsValue
  684. else
  685. if a<b then
  686. result:=LessThanValue;
  687. end;
  688. function CompareValue(const A, B: NativeLargeInt): TValueRelationship;
  689. begin
  690. result:=GreaterThanValue;
  691. if a=b then
  692. result:=EqualsValue
  693. else
  694. if a<b then
  695. result:=LessThanValue;
  696. end;
  697. function CompareValue(const A, B: NativeLargeUInt): TValueRelationship;
  698. begin
  699. result:=GreaterThanValue;
  700. if a=b then
  701. result:=EqualsValue
  702. else
  703. if a<b then
  704. result:=LessThanValue;
  705. end;
  706. function CompareValue(const A, B: Double; delta: Double): TValueRelationship;
  707. begin
  708. result:=GreaterThanValue;
  709. if abs(a-b)<=delta then
  710. result:=EqualsValue
  711. else
  712. if a<b then
  713. result:=LessThanValue;
  714. end;
  715. end.