math.pas 22 KB

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