math.pp 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047
  1. {
  2. $Id$
  3. This file is part of the Free Pascal run time library.
  4. Copyright (c) 1999-2000 by Florian Klaempfl
  5. member of the Free Pascal development team
  6. See the file COPYING.FPC, included in this distribution,
  7. for details about the copyright.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  11. **********************************************************************}
  12. {
  13. This unit is an equivalent to the Delphi math unit
  14. (with some improvements)
  15. What's to do:
  16. o a lot of function :), search for !!!!
  17. o some statistical functions
  18. o all financial functions
  19. o optimizations
  20. }
  21. unit math;
  22. interface
  23. {$MODE objfpc}
  24. {$ifdef VER1_0}
  25. { we don't assume cross compiling from 1.0.x-m68k ... }
  26. {$define FPC_HAS_TYPE_EXTENDED}
  27. {$endif VER1_0}
  28. uses
  29. sysutils;
  30. { Ranges of the IEEE floating point types, including denormals }
  31. {$ifdef FPC_HAS_TYPE_SINGLE}
  32. const
  33. MinSingle = 1.5e-45;
  34. MaxSingle = 3.4e+38;
  35. {$endif FPC_HAS_TYPE_SINGLE}
  36. {$ifdef FPC_HAS_TYPE_DOUBLE}
  37. const
  38. MinDouble = 5.0e-324;
  39. MaxDouble = 1.7e+308;
  40. {$endif FPC_HAS_TYPE_DOUBLE}
  41. {$ifdef FPC_HAS_TYPE_EXTENDED}
  42. const
  43. MinExtended = 3.4e-4932;
  44. MaxExtended = 1.1e+4932;
  45. {$endif FPC_HAS_TYPE_EXTENDED}
  46. {$ifdef FPC_HAS_TYPE_COMP}
  47. const
  48. MinComp = -9.223372036854775807e+18;
  49. MaxComp = 9.223372036854775807e+18;
  50. {$endif FPC_HAS_TYPE_COMP}
  51. { the original delphi functions use extended as argument, }
  52. { but I would prefer double, because 8 bytes is a very }
  53. { natural size for the processor }
  54. { WARNING : changing float type will }
  55. { break all assembler code PM }
  56. {$ifdef FPC_HAS_TYPE_FLOAT128}
  57. type
  58. float = float128;
  59. const
  60. MinFloat = MinFloat128;
  61. MaxFloat = MaxFloat128;
  62. {$else FPC_HAS_TYPE_FLOAT128}
  63. {$ifdef FPC_HAS_TYPE_EXTENDED}
  64. type
  65. float = extended;
  66. const
  67. MinFloat = MinExtended;
  68. MaxFloat = MaxExtended;
  69. {$else FPC_HAS_TYPE_EXTENDED}
  70. {$ifdef FPC_HAS_TYPE_DOUBLE}
  71. type
  72. float = double;
  73. const
  74. MinFloat = MinDouble;
  75. MaxFloat = MaxDouble;
  76. {$else FPC_HAS_TYPE_DOUBLE}
  77. {$ifdef FPC_HAS_TYPE_SINGLE}
  78. type
  79. float = single;
  80. const
  81. MinFloat = MinSingle;
  82. MaxFloat = MaxSingle;
  83. {$else FPC_HAS_TYPE_SINGLE}
  84. {$fatal At least one floating point type must be supported}
  85. {$endif FPC_HAS_TYPE_SINGLE}
  86. {$endif FPC_HAS_TYPE_DOUBLE}
  87. {$endif FPC_HAS_TYPE_EXTENDED}
  88. {$endif FPC_HAS_TYPE_FLOAT128}
  89. type
  90. PFloat = ^Float;
  91. PInteger = ^Integer;
  92. tpaymenttime = (ptendofperiod,ptstartofperiod);
  93. einvalidargument = class(ematherror);
  94. TValueRelationship = -1..1;
  95. const
  96. EqualsValue = 0;
  97. LessThanValue = Low(TValueRelationship);
  98. GreaterThanValue = High(TValueRelationship);
  99. { Min/max determination }
  100. function MinIntValue(const Data: array of Integer): Integer;
  101. function MaxIntValue(const Data: array of Integer): Integer;
  102. { Extra, not present in Delphi, but used frequently }
  103. function Min(a, b: Integer): Integer;
  104. function Max(a, b: Integer): Integer;
  105. function Min(a, b: Cardinal): Cardinal;
  106. function Max(a, b: Cardinal): Cardinal;
  107. function Min(a, b: Int64): Int64;
  108. function Max(a, b: Int64): Int64;
  109. {$ifdef FPC_HAS_TYPE_SINGLE}
  110. function Min(a, b: Single): Single;
  111. function Max(a, b: Single): Single;
  112. {$endif FPC_HAS_TYPE_SINGLE}
  113. {$ifdef FPC_HAS_TYPE_DOUBLE}
  114. function Min(a, b: Double): Double;
  115. function Max(a, b: Double): Double;
  116. {$endif FPC_HAS_TYPE_DOUBLE}
  117. {$ifdef FPC_HAS_TYPE_EXTENDED}
  118. function Min(a, b: Extended): Extended;
  119. function Max(a, b: Extended): Extended;
  120. {$endif FPC_HAS_TYPE_EXTENDED}
  121. { angle conversion }
  122. function degtorad(deg : float) : float;
  123. function radtodeg(rad : float) : float;
  124. function gradtorad(grad : float) : float;
  125. function radtograd(rad : float) : float;
  126. function degtograd(deg : float) : float;
  127. function gradtodeg(grad : float) : float;
  128. { one cycle are 2*Pi rad }
  129. function cycletorad(cycle : float) : float;
  130. function radtocycle(rad : float) : float;
  131. { trigoniometric functions }
  132. function tan(x : float) : float;
  133. function cotan(x : float) : float;
  134. procedure sincos(theta : float;var sinus,cosinus : float);
  135. { inverse functions }
  136. function arccos(x : float) : float;
  137. function arcsin(x : float) : float;
  138. { calculates arctan(y/x) and returns an angle in the correct quadrant }
  139. function arctan2(y,x : float) : float;
  140. { hyperbolic functions }
  141. function cosh(x : float) : float;
  142. function sinh(x : float) : float;
  143. function tanh(x : float) : float;
  144. { area functions }
  145. { delphi names: }
  146. function arccosh(x : float) : float;
  147. function arcsinh(x : float) : float;
  148. function arctanh(x : float) : float;
  149. { IMHO the function should be called as follows (FK) }
  150. function arcosh(x : float) : float;
  151. function arsinh(x : float) : float;
  152. function artanh(x : float) : float;
  153. { triangle functions }
  154. { returns the length of the hypotenuse of a right triangle }
  155. { if x and y are the other sides }
  156. function hypot(x,y : float) : float;
  157. { logarithm functions }
  158. function log10(x : float) : float;
  159. function log2(x : float) : float;
  160. function logn(n,x : float) : float;
  161. { returns natural logarithm of x+1 }
  162. function lnxp1(x : float) : float;
  163. { exponential functions }
  164. function power(base,exponent : float) : float;
  165. { base^exponent }
  166. function intpower(base : float;const exponent : Integer) : float;
  167. { number converting }
  168. { rounds x towards positive infinity }
  169. function ceil(x : float) : Integer;
  170. { rounds x towards negative infinity }
  171. function floor(x : float) : Integer;
  172. { misc. functions }
  173. { splits x into mantissa and exponent (to base 2) }
  174. procedure Frexp(X: float; var Mantissa: float; var Exponent: integer);
  175. { returns x*(2^p) }
  176. function ldexp(x : float; const p : Integer) : float;
  177. { statistical functions }
  178. function mean(const data : array of float) : float;
  179. function sum(const data : array of float) : float;
  180. function mean(const data : PFloat; Const N : longint) : float;
  181. function sum(const data : PFloat; Const N : Longint) : float;
  182. function sumofsquares(const data : array of float) : float;
  183. function sumofsquares(const data : PFloat; Const N : Integer) : float;
  184. { calculates the sum and the sum of squares of data }
  185. procedure sumsandsquares(const data : array of float;
  186. var sum,sumofsquares : float);
  187. procedure sumsandsquares(const data : PFloat; Const N : Integer;
  188. var sum,sumofsquares : float);
  189. function minvalue(const data : array of float) : float;
  190. function minvalue(const data : array of integer) : Integer;
  191. function minvalue(const data : PFloat; Const N : Integer) : float;
  192. function MinValue(const Data : PInteger; Const N : Integer): Integer;
  193. function maxvalue(const data : array of float) : float;
  194. function maxvalue(const data : array of integer) : Integer;
  195. function maxvalue(const data : PFloat; Const N : Integer) : float;
  196. function maxvalue(const data : PInteger; Const N : Integer) : Integer;
  197. { calculates the standard deviation }
  198. function stddev(const data : array of float) : float;
  199. function stddev(const data : PFloat; Const N : Integer) : float;
  200. { calculates the mean and stddev }
  201. procedure meanandstddev(const data : array of float;
  202. var mean,stddev : float);
  203. procedure meanandstddev(const data : PFloat;
  204. Const N : Longint;var mean,stddev : float);
  205. function variance(const data : array of float) : float;
  206. function totalvariance(const data : array of float) : float;
  207. function variance(const data : PFloat; Const N : Integer) : float;
  208. function totalvariance(const data : PFloat; Const N : Integer) : float;
  209. { returns random values with gaussian distribution }
  210. function randg(mean,stddev : float) : float;
  211. { I don't know what the following functions do: }
  212. function popnstddev(const data : array of float) : float;
  213. function popnstddev(const data : PFloat; Const N : Integer) : float;
  214. function popnvariance(const data : PFloat; Const N : Integer) : float;
  215. function popnvariance(const data : array of float) : float;
  216. procedure momentskewkurtosis(const data : array of float;
  217. var m1,m2,m3,m4,skew,kurtosis : float);
  218. procedure momentskewkurtosis(const data : PFloat; Const N : Integer;
  219. var m1,m2,m3,m4,skew,kurtosis : float);
  220. { geometrical function }
  221. { returns the euclidean L2 norm }
  222. function norm(const data : array of float) : float;
  223. function norm(const data : PFloat; Const N : Integer) : float;
  224. { include cpu specific stuff }
  225. {$i mathuh.inc}
  226. implementation
  227. { include cpu specific stuff }
  228. {$i mathu.inc}
  229. ResourceString
  230. SMathError = 'Math Error : %s';
  231. SInvalidArgument = 'Invalid argument';
  232. Procedure DoMathError(Const S : String);
  233. begin
  234. Raise EMathError.CreateFmt(SMathError,[S]);
  235. end;
  236. Procedure InvalidArgument;
  237. begin
  238. Raise EInvalidArgument.Create(SInvalidArgument);
  239. end;
  240. function degtorad(deg : float) : float;
  241. begin
  242. degtorad:=deg*(pi/180.0);
  243. end;
  244. function radtodeg(rad : float) : float;
  245. begin
  246. radtodeg:=rad*(180.0/pi);
  247. end;
  248. function gradtorad(grad : float) : float;
  249. begin
  250. gradtorad:=grad*(pi/200.0);
  251. end;
  252. function radtograd(rad : float) : float;
  253. begin
  254. radtograd:=rad*(200.0/pi);
  255. end;
  256. function degtograd(deg : float) : float;
  257. begin
  258. degtograd:=deg*(200.0/180.0);
  259. end;
  260. function gradtodeg(grad : float) : float;
  261. begin
  262. gradtodeg:=grad*(180.0/200.0);
  263. end;
  264. function cycletorad(cycle : float) : float;
  265. begin
  266. cycletorad:=(2*pi)*cycle;
  267. end;
  268. function radtocycle(rad : float) : float;
  269. begin
  270. { avoid division }
  271. radtocycle:=rad*(1/(2*pi));
  272. end;
  273. function tan(x : float) : float;
  274. begin
  275. Tan:=Sin(x)/Cos(x)
  276. end;
  277. function cotan(x : float) : float;
  278. begin
  279. cotan:=Cos(X)/Sin(X);
  280. end;
  281. procedure sincos(theta : float;var sinus,cosinus : float);
  282. begin
  283. sinus:=sin(theta);
  284. cosinus:=cos(theta);
  285. end;
  286. { Sign, ArcSin and ArcCos from Arjan van Dijk ([email protected]) }
  287. function sign(x : float) : float;
  288. begin
  289. if x > 0 then sign := 1.0
  290. else if x < 0 then sign := -1.0
  291. else sign := 0.0;
  292. end;
  293. function arcsin(x : float) : float;
  294. begin
  295. if abs(x) > 1 then InvalidArgument
  296. else if abs(x) < 0.5 then
  297. arcsin := arctan(x/sqrt(1-sqr(x)))
  298. else
  299. arcsin := sign(x) * (pi*0.5 - arctan(sqrt(1 / sqr(x) - 1)));
  300. end;
  301. function Arccos(x : Float) : Float;
  302. begin
  303. arccos := pi*0.5 - arcsin(x);
  304. end;
  305. {$ifndef FPC_MATH_HAS_ARCTAN2}
  306. function arctan2(y,x : float) : float;
  307. begin
  308. if (x=0) then
  309. begin
  310. if y=0 then
  311. arctan2:=0.0
  312. else if y>0 then
  313. arctan2:=pi/2
  314. else if y<0 then
  315. arctan2:=-pi/2;
  316. end
  317. else
  318. ArcTan2:=ArcTan(y/x);
  319. end;
  320. {$endif FPC_MATH_HAS_ARCTAN2}
  321. function cosh(x : float) : float;
  322. var
  323. temp : float;
  324. begin
  325. temp:=exp(x);
  326. cosh:=0.5*(temp+1.0/temp);
  327. end;
  328. function sinh(x : float) : float;
  329. var
  330. temp : float;
  331. begin
  332. temp:=exp(x);
  333. sinh:=0.5*(temp-1.0/temp);
  334. end;
  335. Const MaxTanh = 5678.22249441322; // Ln(MaxExtended)/2
  336. function tanh(x : float) : float;
  337. var Temp : float;
  338. begin
  339. if x>MaxTanh then exit(1.0)
  340. else if x<-MaxTanh then exit (-1.0);
  341. temp:=exp(-2*x);
  342. tanh:=(1-temp)/(1+temp)
  343. end;
  344. function arccosh(x : float) : float;
  345. begin
  346. arccosh:=arcosh(x);
  347. end;
  348. function arcsinh(x : float) : float;
  349. begin
  350. arcsinh:=arsinh(x);
  351. end;
  352. function arctanh(x : float) : float;
  353. begin
  354. if x>1 then InvalidArgument;
  355. arctanh:=artanh(x);
  356. end;
  357. function arcosh(x : float) : float;
  358. begin
  359. if x<1 then InvalidArgument;
  360. arcosh:=Ln(x+Sqrt(x*x-1));
  361. end;
  362. function arsinh(x : float) : float;
  363. begin
  364. arsinh:=Ln(x+Sqrt(1+x*x));
  365. end;
  366. function artanh(x : float) : float;
  367. begin
  368. If abs(x)>1 then InvalidArgument;
  369. artanh:=(Ln((1+x)/(1-x)))*0.5;
  370. end;
  371. function hypot(x,y : float) : float;
  372. begin
  373. hypot:=Sqrt(x*x+y*y)
  374. end;
  375. function log10(x : float) : float;
  376. begin
  377. log10:=ln(x)/ln(10);
  378. end;
  379. function log2(x : float) : float;
  380. begin
  381. log2:=ln(x)/ln(2)
  382. end;
  383. function logn(n,x : float) : float;
  384. begin
  385. if n<0 then InvalidArgument;
  386. logn:=ln(x)/ln(n);
  387. end;
  388. function lnxp1(x : float) : float;
  389. begin
  390. if x<-1 then
  391. InvalidArgument;
  392. lnxp1:=ln(1+x);
  393. end;
  394. function power(base,exponent : float) : float;
  395. begin
  396. If Exponent=0.0 then
  397. Result:=1.0
  398. else
  399. If base>0.0 then
  400. Power:=exp(exponent * ln (base))
  401. else if base=0.0 then
  402. Result:=0.0
  403. else
  404. InvalidArgument
  405. end;
  406. function intpower(base : float;const exponent : Integer) : float;
  407. var
  408. i : longint;
  409. begin
  410. i:=abs(exponent);
  411. intpower:=1.0;
  412. while i>0 do
  413. begin
  414. while (i and 1)=0 do
  415. begin
  416. i:=i shr 1;
  417. base:=sqr(base);
  418. end;
  419. i:=i-1;
  420. intpower:=intpower*base;
  421. end;
  422. if exponent<0 then
  423. intpower:=1.0/intpower;
  424. end;
  425. function ceil(x : float) : integer;
  426. begin
  427. Ceil:=Trunc(x);
  428. If Frac(x)>0 then
  429. Ceil:=Ceil+1;
  430. end;
  431. function floor(x : float) : integer;
  432. begin
  433. Floor:=Trunc(x);
  434. If Frac(x)<0 then
  435. Floor := Floor-1;
  436. end;
  437. procedure Frexp(X: float; var Mantissa: float; var Exponent: integer);
  438. begin
  439. Exponent :=0;
  440. if (abs(x)<0.5) then
  441. While (abs(x)<0.5) do
  442. begin
  443. x := x*2;
  444. Dec(Exponent);
  445. end
  446. else
  447. While (abs(x)>1) do
  448. begin
  449. x := x/2;
  450. Inc(Exponent);
  451. end;
  452. mantissa := x;
  453. end;
  454. function ldexp(x : float;const p : Integer) : float;
  455. begin
  456. ldexp:=x*intpower(2.0,p);
  457. end;
  458. function mean(const data : array of float) : float;
  459. begin
  460. Result:=Mean(@data[0],High(Data)+1);
  461. end;
  462. function mean(const data : PFloat; Const N : longint) : float;
  463. begin
  464. mean:=sum(Data,N);
  465. mean:=mean/N;
  466. end;
  467. function sum(const data : array of float) : float;
  468. begin
  469. Result:=Sum(@Data[0],High(Data)+1);
  470. end;
  471. function sum(const data : PFloat;Const N : longint) : float;
  472. var
  473. i : longint;
  474. begin
  475. sum:=0.0;
  476. for i:=0 to N-1 do
  477. sum:=sum+data[i];
  478. end;
  479. function sumofsquares(const data : array of float) : float;
  480. begin
  481. Result:=sumofsquares(@data[0],High(Data)+1);
  482. end;
  483. function sumofsquares(const data : PFloat; Const N : Integer) : float;
  484. var
  485. i : longint;
  486. begin
  487. sumofsquares:=0.0;
  488. for i:=0 to N-1 do
  489. sumofsquares:=sumofsquares+sqr(data[i]);
  490. end;
  491. procedure sumsandsquares(const data : array of float;
  492. var sum,sumofsquares : float);
  493. begin
  494. sumsandsquares (@Data[0],High(Data)+1,Sum,sumofsquares);
  495. end;
  496. procedure sumsandsquares(const data : PFloat; Const N : Integer;
  497. var sum,sumofsquares : float);
  498. var
  499. i : Integer;
  500. temp : float;
  501. begin
  502. sumofsquares:=0.0;
  503. sum:=0.0;
  504. for i:=0 to N-1 do
  505. begin
  506. temp:=data[i];
  507. sumofsquares:=sumofsquares+sqr(temp);
  508. sum:=sum+temp;
  509. end;
  510. end;
  511. function stddev(const data : array of float) : float;
  512. begin
  513. Result:=Stddev(@Data[0],High(Data)+1)
  514. end;
  515. function stddev(const data : PFloat; Const N : Integer) : float;
  516. begin
  517. StdDev:=Sqrt(Variance(Data,N));
  518. end;
  519. procedure meanandstddev(const data : array of float;
  520. var mean,stddev : float);
  521. begin
  522. Meanandstddev(@Data[0],High(Data)+1,Mean,stddev);
  523. end;
  524. procedure meanandstddev(const data : PFloat;
  525. Const N : Longint;var mean,stddev : float);
  526. Var I : longint;
  527. begin
  528. Mean:=0;
  529. StdDev:=0;
  530. For I:=0 to N-1 do
  531. begin
  532. Mean:=Mean+Data[i];
  533. StdDev:=StdDev+Sqr(Data[i]);
  534. end;
  535. Mean:=Mean/N;
  536. StdDev:=(StdDev-N*Sqr(Mean));
  537. If N>1 then
  538. StdDev:=Sqrt(Stddev/(N-1))
  539. else
  540. StdDev:=0;
  541. end;
  542. function variance(const data : array of float) : float;
  543. begin
  544. Variance:=Variance(@Data[0],High(Data)+1);
  545. end;
  546. function variance(const data : PFloat; Const N : Integer) : float;
  547. begin
  548. If N=1 then
  549. Result:=0
  550. else
  551. Result:=TotalVariance(Data,N)/(N-1);
  552. end;
  553. function totalvariance(const data : array of float) : float;
  554. begin
  555. Result:=TotalVariance(@Data[0],High(Data)+1);
  556. end;
  557. function totalvariance(const data : Pfloat;Const N : Integer) : float;
  558. var S,SS : Float;
  559. begin
  560. If N=1 then
  561. Result:=0
  562. else
  563. begin
  564. SumsAndSquares(Data,N,S,SS);
  565. Result := SS-Sqr(S)/N;
  566. end;
  567. end;
  568. function randg(mean,stddev : float) : float;
  569. Var U1,S2 : Float;
  570. begin
  571. repeat
  572. u1:= 2*random-1;
  573. S2:=Sqr(U1)+sqr(2*random-1);
  574. until s2<1;
  575. randg:=Sqrt(-2*ln(S2)/S2)*u1*stddev+Mean;
  576. end;
  577. function popnstddev(const data : array of float) : float;
  578. begin
  579. PopnStdDev:=Sqrt(PopnVariance(@Data[0],High(Data)+1));
  580. end;
  581. function popnstddev(const data : PFloat; Const N : Integer) : float;
  582. begin
  583. PopnStdDev:=Sqrt(PopnVariance(Data,N));
  584. end;
  585. function popnvariance(const data : array of float) : float;
  586. begin
  587. popnvariance:=popnvariance(@data[0],high(Data)+1);
  588. end;
  589. function popnvariance(const data : PFloat; Const N : Integer) : float;
  590. begin
  591. PopnVariance:=TotalVariance(Data,N)/N;
  592. end;
  593. procedure momentskewkurtosis(const data : array of float;
  594. var m1,m2,m3,m4,skew,kurtosis : float);
  595. begin
  596. momentskewkurtosis(@Data[0],High(Data)+1,m1,m2,m3,m4,skew,kurtosis);
  597. end;
  598. procedure momentskewkurtosis(const data : PFloat; Const N : Integer;
  599. var m1,m2,m3,m4,skew,kurtosis : float);
  600. Var S,SS,SC,SQ,invN,Acc,M1S,S2N,S3N,temp : Float;
  601. I : Longint;
  602. begin
  603. invN:=1.0/N;
  604. s:=0;
  605. ss:=0;
  606. sq:=0;
  607. sc:=0;
  608. for i:=0 to N-1 do
  609. begin
  610. temp:=Data[i]; { faster }
  611. S:=S+temp;
  612. acc:=temp*temp;
  613. ss:=ss+acc;
  614. Acc:=acc*temp;
  615. Sc:=sc+acc;
  616. acc:=acc*temp;
  617. sq:=sq+acc;
  618. end;
  619. M1:=s*invN;
  620. M1S:=M1*M1;
  621. S2N:=SS*invN;
  622. S3N:=SC*invN;
  623. M2:=S2N-M1S;
  624. M3:=S3N-(M1*3*S2N) + 2*M1S*M1;
  625. M4:=SQ*invN - (M1 * 4 * S3N) + (M1S*6*S2N-3*Sqr(M1S));
  626. Skew:=M3*power(M2,-3/2);
  627. Kurtosis:=M4 / Sqr(M2);
  628. end;
  629. function norm(const data : array of float) : float;
  630. begin
  631. norm:=Norm(@data[0],High(Data)+1);
  632. end;
  633. function norm(const data : PFloat; Const N : Integer) : float;
  634. begin
  635. norm:=sqrt(sumofsquares(data,N));
  636. end;
  637. function MinIntValue(const Data: array of Integer): Integer;
  638. var
  639. I: Integer;
  640. begin
  641. Result := Data[Low(Data)];
  642. For I := Succ(Low(Data)) To High(Data) Do
  643. If Data[I] < Result Then Result := Data[I];
  644. end;
  645. function MinValue(const Data: array of Integer): Integer;
  646. begin
  647. Result:=MinValue(Pinteger(@Data[0]),High(Data)+1)
  648. end;
  649. function MinValue(const Data: PInteger; Const N : Integer): Integer;
  650. var
  651. I: Integer;
  652. begin
  653. Result := Data[0];
  654. For I := 1 To N-1 do
  655. If Data[I] < Result Then Result := Data[I];
  656. end;
  657. function minvalue(const data : array of float) : float;
  658. begin
  659. Result:=minvalue(PFloat(@data[0]),High(Data)+1);
  660. end;
  661. function minvalue(const data : PFloat; Const N : Integer) : float;
  662. var
  663. i : longint;
  664. begin
  665. { get an initial value }
  666. minvalue:=data[0];
  667. for i:=1 to N-1 do
  668. if data[i]<minvalue then
  669. minvalue:=data[i];
  670. end;
  671. function MaxIntValue(const Data: array of Integer): Integer;
  672. var
  673. I: Integer;
  674. begin
  675. Result := Data[Low(Data)];
  676. For I := Succ(Low(Data)) To High(Data) Do
  677. If Data[I] > Result Then Result := Data[I];
  678. end;
  679. function maxvalue(const data : array of float) : float;
  680. begin
  681. Result:=maxvalue(PFloat(@data[0]),High(Data)+1);
  682. end;
  683. function maxvalue(const data : PFloat; Const N : Integer) : float;
  684. var
  685. i : longint;
  686. begin
  687. { get an initial value }
  688. maxvalue:=data[0];
  689. for i:=1 to N-1 do
  690. if data[i]>maxvalue then
  691. maxvalue:=data[i];
  692. end;
  693. function MaxValue(const Data: array of Integer): Integer;
  694. begin
  695. Result:=MaxValue(PInteger(@Data[0]),High(Data)+1)
  696. end;
  697. function maxvalue(const data : PInteger; Const N : Integer) : Integer;
  698. var
  699. i : longint;
  700. begin
  701. { get an initial value }
  702. maxvalue:=data[0];
  703. for i:=1 to N-1 do
  704. if data[i]>maxvalue then
  705. maxvalue:=data[i];
  706. end;
  707. function Min(a, b: Integer): Integer;
  708. begin
  709. if a < b then
  710. Result := a
  711. else
  712. Result := b;
  713. end;
  714. function Max(a, b: Integer): Integer;
  715. begin
  716. if a > b then
  717. Result := a
  718. else
  719. Result := b;
  720. end;
  721. function Min(a, b: Cardinal): Cardinal;
  722. begin
  723. if a < b then
  724. Result := a
  725. else
  726. Result := b;
  727. end;
  728. function Max(a, b: Cardinal): Cardinal;
  729. begin
  730. if a > b then
  731. Result := a
  732. else
  733. Result := b;
  734. end;
  735. function Min(a, b: Int64): Int64;
  736. begin
  737. if a < b then
  738. Result := a
  739. else
  740. Result := b;
  741. end;
  742. function Max(a, b: Int64): Int64;
  743. begin
  744. if a > b then
  745. Result := a
  746. else
  747. Result := b;
  748. end;
  749. {$ifdef FPC_HAS_TYPE_SINGLE}
  750. function Min(a, b: Single): Single;
  751. begin
  752. if a < b then
  753. Result := a
  754. else
  755. Result := b;
  756. end;
  757. function Max(a, b: Single): Single;
  758. begin
  759. if a > b then
  760. Result := a
  761. else
  762. Result := b;
  763. end;
  764. {$endif FPC_HAS_TYPE_SINGLE}
  765. {$ifdef FPC_HAS_TYPE_DOUBLE}
  766. function Min(a, b: Double): Double;
  767. begin
  768. if a < b then
  769. Result := a
  770. else
  771. Result := b;
  772. end;
  773. function Max(a, b: Double): Double;
  774. begin
  775. if a > b then
  776. Result := a
  777. else
  778. Result := b;
  779. end;
  780. {$endif FPC_HAS_TYPE_DOUBLE}
  781. {$ifdef FPC_HAS_TYPE_EXTENDED}
  782. function Min(a, b: Extended): Extended;
  783. begin
  784. if a < b then
  785. Result := a
  786. else
  787. Result := b;
  788. end;
  789. function Max(a, b: Extended): Extended;
  790. begin
  791. if a > b then
  792. Result := a
  793. else
  794. Result := b;
  795. end;
  796. {$endif FPC_HAS_TYPE_EXTENDED}
  797. end.
  798. {
  799. $Log$
  800. Revision 1.14 2003-10-29 19:10:07 jonas
  801. * fixed arctan2
  802. Revision 1.13 2003/10/26 15:58:05 florian
  803. * fixed arctan2 to handle x=0 correctly as well
  804. Revision 1.12 2003/09/01 20:46:59 peter
  805. * small fixes for sparc
  806. Revision 1.11 2003/04/24 09:38:12 florian
  807. * min/max must check the compiler capabilities
  808. Revision 1.10 2003/04/24 09:21:59 florian
  809. + moved cpu dependend code to mathuh.inc and mathu.inc
  810. Revision 1.9 2003/01/03 20:34:02 peter
  811. * i386 fpu controlword functions added
  812. Revision 1.8 2002/09/07 21:06:12 carl
  813. * cleanup of parameters
  814. - remove assembler code
  815. Revision 1.7 2002/09/07 16:01:22 peter
  816. * old logs removed and tabs fixed
  817. }