math.pp 96 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397339833993400340134023403340434053406340734083409341034113412341334143415341634173418341934203421342234233424342534263427342834293430343134323433343434353436343734383439344034413442344334443445344634473448344934503451345234533454345534563457345834593460346134623463346434653466346734683469347034713472347334743475347634773478347934803481348234833484348534863487348834893490349134923493349434953496349734983499350035013502350335043505350635073508350935103511351235133514351535163517351835193520352135223523352435253526352735283529353035313532353335343535353635373538353935403541354235433544354535463547354835493550355135523553355435553556355735583559356035613562356335643565356635673568356935703571357235733574357535763577357835793580358135823583358435853586358735883589359035913592359335943595359635973598359936003601360236033604360536063607360836093610361136123613361436153616361736183619362036213622362336243625362636273628362936303631363236333634363536363637363836393640364136423643364436453646364736483649365036513652365336543655365636573658365936603661366236633664366536663667366836693670367136723673367436753676367736783679368036813682368336843685368636873688368936903691369236933694369536963697369836993700370137023703370437053706370737083709371037113712371337143715371637173718371937203721372237233724372537263727372837293730373137323733373437353736373737383739374037413742374337443745374637473748374937503751375237533754375537563757375837593760376137623763376437653766376737683769377037713772377337743775377637773778377937803781378237833784378537863787378837893790379137923793379437953796379737983799380038013802380338043805380638073808380938103811381238133814381538163817381838193820382138223823382438253826382738283829383038313832383338343835
  1. {
  2. This file is part of the Free Pascal run time library.
  3. Copyright (c) 1999-2005 by Florian Klaempfl
  4. member of the Free Pascal development team
  5. See the file COPYING.FPC, included in this distribution,
  6. for details about the copyright.
  7. This program is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY; without even the implied warranty of
  9. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  10. **********************************************************************}
  11. {-------------------------------------------------------------------------
  12. Using functions from AMath/DAMath libraries, which are covered by the
  13. following license:
  14. (C) Copyright 2009-2013 Wolfgang Ehrhardt
  15. This software is provided 'as-is', without any express or implied warranty.
  16. In no event will the authors be held liable for any damages arising from
  17. the use of this software.
  18. Permission is granted to anyone to use this software for any purpose,
  19. including commercial applications, and to alter it and redistribute it
  20. freely, subject to the following restrictions:
  21. 1. The origin of this software must not be misrepresented; you must not
  22. claim that you wrote the original software. If you use this software in
  23. a product, an acknowledgment in the product documentation would be
  24. appreciated but is not required.
  25. 2. Altered source versions must be plainly marked as such, and must not be
  26. misrepresented as being the original software.
  27. 3. This notice may not be removed or altered from any source distribution.
  28. ----------------------------------------------------------------------------}
  29. {
  30. This unit is an equivalent to the Delphi Math unit
  31. (with some improvements)
  32. What's to do:
  33. o some statistical functions
  34. o optimizations
  35. }
  36. {$MODE objfpc}
  37. {$inline on }
  38. {$GOTO on}
  39. {$IFNDEF FPC_DOTTEDUNITS}
  40. unit Math;
  41. {$ENDIF FPC_DOTTEDUNITS}
  42. interface
  43. {$ifndef FPUNONE}
  44. {$IFDEF FPC_DOTTEDUNITS}
  45. uses
  46. System.SysUtils, System.Types;
  47. {$ELSE FPC_DOTTEDUNITS}
  48. uses
  49. sysutils, types;
  50. {$ENDIF FPC_DOTTEDUNITS}
  51. {$IFDEF FPDOC_MATH}
  52. Type
  53. Float = MaxFloatType;
  54. Const
  55. MinFloat = 0;
  56. MaxFloat = 0;
  57. {$ENDIF}
  58. { Ranges of the IEEE floating point types, including denormals }
  59. {$ifdef FPC_HAS_TYPE_SINGLE}
  60. const
  61. { values according to
  62. https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Single-precision_examples
  63. }
  64. MinSingle = 1.1754943508e-38;
  65. MaxSingle = 3.4028234664e+38;
  66. {$endif FPC_HAS_TYPE_SINGLE}
  67. {$ifdef FPC_HAS_TYPE_DOUBLE}
  68. const
  69. { values according to
  70. https://en.wikipedia.org/wiki/Double-precision_floating-point_format#Double-precision_examples
  71. }
  72. MinDouble = 2.2250738585072014e-308;
  73. MaxDouble = 1.7976931348623157e+308;
  74. {$endif FPC_HAS_TYPE_DOUBLE}
  75. {$ifdef FPC_HAS_TYPE_EXTENDED}
  76. const
  77. MinExtended = 3.4e-4932;
  78. MaxExtended = 1.1e+4932;
  79. {$endif FPC_HAS_TYPE_EXTENDED}
  80. {$ifdef FPC_HAS_TYPE_COMP}
  81. const
  82. MinComp = -9.223372036854775807e+18;
  83. MaxComp = 9.223372036854775807e+18;
  84. {$endif FPC_HAS_TYPE_COMP}
  85. { the original delphi functions use extended as argument, }
  86. { but I would prefer double, because 8 bytes is a very }
  87. { natural size for the processor }
  88. { WARNING : changing float type will }
  89. { break all assembler code PM }
  90. {$if defined(FPC_HAS_TYPE_FLOAT128)}
  91. type
  92. Float = Float128;
  93. const
  94. MinFloat = MinFloat128;
  95. MaxFloat = MaxFloat128;
  96. {$elseif defined(FPC_HAS_TYPE_EXTENDED)}
  97. type
  98. Float = extended;
  99. const
  100. MinFloat = MinExtended;
  101. MaxFloat = MaxExtended;
  102. {$elseif defined(FPC_HAS_TYPE_DOUBLE)}
  103. type
  104. Float = double;
  105. const
  106. MinFloat = MinDouble;
  107. MaxFloat = MaxDouble;
  108. {$elseif defined(FPC_HAS_TYPE_SINGLE)}
  109. type
  110. Float = single;
  111. const
  112. MinFloat = MinSingle;
  113. MaxFloat = MaxSingle;
  114. {$else}
  115. {$fatal At least one floating point type must be supported}
  116. {$endif}
  117. type
  118. PFloat = ^Float;
  119. PInteger = ObjPas.PInteger;
  120. TPaymentTime = (ptEndOfPeriod,ptStartOfPeriod);
  121. EInvalidArgument = class(ematherror);
  122. {$IFDEF FPC_DOTTEDUNITS}
  123. TValueRelationship = System.Types.TValueRelationship;
  124. {$ELSE FPC_DOTTEDUNITS}
  125. TValueRelationship = types.TValueRelationship;
  126. {$ENDIF FPC_DOTTEDUNITS}
  127. const
  128. {$IFDEF FPC_DOTTEDUNITS}
  129. EqualsValue = System.Types.EqualsValue;
  130. LessThanValue = System.Types.LessThanValue;
  131. GreaterThanValue = System.Types.GreaterThanValue;
  132. {$ELSE FPC_DOTTEDUNITS}
  133. EqualsValue = types.EqualsValue;
  134. LessThanValue = types.LessThanValue;
  135. GreaterThanValue = types.GreaterThanValue;
  136. {$ENDIF FPC_DOTTEDUNITS}
  137. {$push}
  138. {$R-}
  139. {$Q-}
  140. NaN = 0.0/0.0;
  141. Infinity = 1.0/0.0;
  142. NegInfinity = -1.0/0.0;
  143. {$pop}
  144. {$IFDEF FPDOC_MATH}
  145. // This must be after the above defines.
  146. {$DEFINE FPC_HAS_TYPE_SINGLE}
  147. {$DEFINE FPC_HAS_TYPE_DOUBLE}
  148. {$DEFINE FPC_HAS_TYPE_EXTENDED}
  149. {$DEFINE FPC_HAS_TYPE_COMP}
  150. {$ENDIF}
  151. { Min/max determination }
  152. function MinIntValue(const Data: array of Integer): Integer;
  153. function MaxIntValue(const Data: array of Integer): Integer;
  154. { Extra, not present in Delphi, but used frequently }
  155. function Min(a, b: Integer): Integer;inline; overload;
  156. function Max(a, b: Integer): Integer;inline; overload;
  157. { this causes more trouble than it solves
  158. function Min(a, b: Cardinal): Cardinal; overload;
  159. function Max(a, b: Cardinal): Cardinal; overload;
  160. }
  161. function Min(a, b: Int64): Int64;inline; overload;
  162. function Max(a, b: Int64): Int64;inline; overload;
  163. function Min(a, b: QWord): QWord;inline; overload;
  164. function Max(a, b: QWord): QWord;inline; overload;
  165. {$ifdef FPC_HAS_TYPE_SINGLE}
  166. function Min(a, b: Single): Single;inline; overload;
  167. function Max(a, b: Single): Single;inline; overload;
  168. {$endif FPC_HAS_TYPE_SINGLE}
  169. {$ifdef FPC_HAS_TYPE_DOUBLE}
  170. function Min(a, b: Double): Double;inline; overload;
  171. function Max(a, b: Double): Double;inline; overload;
  172. {$endif FPC_HAS_TYPE_DOUBLE}
  173. {$ifdef FPC_HAS_TYPE_EXTENDED}
  174. function Min(a, b: Extended): Extended;inline; overload;
  175. function Max(a, b: Extended): Extended;inline; overload;
  176. {$endif FPC_HAS_TYPE_EXTENDED}
  177. function InRange(const AValue, AMin, AMax: Integer): Boolean;inline; overload;
  178. function InRange(const AValue, AMin, AMax: Int64): Boolean;inline; overload;
  179. {$ifdef FPC_HAS_TYPE_DOUBLE}
  180. function InRange(const AValue, AMin, AMax: Double): Boolean;inline; overload;
  181. {$endif FPC_HAS_TYPE_DOUBLE}
  182. function EnsureRange(const AValue, AMin, AMax: Integer): Integer;inline; overload;
  183. function EnsureRange(const AValue, AMin, AMax: Int64): Int64;inline; overload;
  184. {$ifdef FPC_HAS_TYPE_DOUBLE}
  185. function EnsureRange(const AValue, AMin, AMax: Double): Double;inline; overload;
  186. {$endif FPC_HAS_TYPE_DOUBLE}
  187. procedure DivMod(Dividend: LongInt; Divisor: Word; var Result, Remainder: Word);
  188. procedure DivMod(Dividend: LongInt; Divisor: Word; var Result, Remainder: SmallInt);
  189. procedure DivMod(Dividend: DWord; Divisor: DWord; var Result, Remainder: DWord);
  190. procedure DivMod(Dividend: LongInt; Divisor: LongInt; var Result, Remainder: LongInt);
  191. { Floating point modulo}
  192. {$ifdef FPC_HAS_TYPE_SINGLE}
  193. function FMod(const a, b: Single): Single;inline;overload;
  194. {$endif FPC_HAS_TYPE_SINGLE}
  195. {$ifdef FPC_HAS_TYPE_DOUBLE}
  196. function FMod(const a, b: Double): Double;inline;overload;
  197. {$endif FPC_HAS_TYPE_DOUBLE}
  198. {$ifdef FPC_HAS_TYPE_EXTENDED}
  199. function FMod(const a, b: Extended): Extended;inline;overload;
  200. {$endif FPC_HAS_TYPE_EXTENDED}
  201. operator mod(const a,b:float) c:float;inline;
  202. // Sign functions
  203. Type
  204. TValueSign = -1..1;
  205. const
  206. NegativeValue = Low(TValueSign);
  207. ZeroValue = 0;
  208. PositiveValue = High(TValueSign);
  209. function Sign(const AValue: Integer): TValueSign;inline; overload;
  210. function Sign(const AValue: Int64): TValueSign;inline; overload;
  211. {$ifdef FPC_HAS_TYPE_SINGLE}
  212. function Sign(const AValue: Single): TValueSign;inline; overload;
  213. {$endif}
  214. function Sign(const AValue: Double): TValueSign;inline; overload;
  215. {$ifdef FPC_HAS_TYPE_EXTENDED}
  216. function Sign(const AValue: Extended): TValueSign;inline; overload;
  217. {$endif}
  218. function IsZero(const A: Single; Epsilon: Single): Boolean; overload;
  219. function IsZero(const A: Single): Boolean;inline; overload;
  220. {$ifdef FPC_HAS_TYPE_DOUBLE}
  221. function IsZero(const A: Double; Epsilon: Double): Boolean; overload;
  222. function IsZero(const A: Double): Boolean;inline; overload;
  223. {$endif FPC_HAS_TYPE_DOUBLE}
  224. {$ifdef FPC_HAS_TYPE_EXTENDED}
  225. function IsZero(const A: Extended; Epsilon: Extended): Boolean; overload;
  226. function IsZero(const A: Extended): Boolean;inline; overload;
  227. {$endif FPC_HAS_TYPE_EXTENDED}
  228. function IsNan(const d : Single): Boolean; overload;
  229. {$ifdef FPC_HAS_TYPE_DOUBLE}
  230. function IsNan(const d : Double): Boolean; overload;
  231. {$endif FPC_HAS_TYPE_DOUBLE}
  232. {$ifdef FPC_HAS_TYPE_EXTENDED}
  233. function IsNan(const d : Extended): Boolean; overload;
  234. {$endif FPC_HAS_TYPE_EXTENDED}
  235. function IsInfinite(const d : Single): Boolean; overload;
  236. {$ifdef FPC_HAS_TYPE_DOUBLE}
  237. function IsInfinite(const d : Double): Boolean; overload;
  238. {$endif FPC_HAS_TYPE_DOUBLE}
  239. {$ifdef FPC_HAS_TYPE_EXTENDED}
  240. function IsInfinite(const d : Extended): Boolean; overload;
  241. {$endif FPC_HAS_TYPE_EXTENDED}
  242. {$ifdef FPC_HAS_TYPE_EXTENDED}
  243. function SameValue(const A, B: Extended): Boolean;inline; overload;
  244. {$endif}
  245. {$ifdef FPC_HAS_TYPE_DOUBLE}
  246. function SameValue(const A, B: Double): Boolean;inline; overload;
  247. {$endif}
  248. function SameValue(const A, B: Single): Boolean;inline; overload;
  249. {$ifdef FPC_HAS_TYPE_EXTENDED}
  250. function SameValue(const A, B: Extended; Epsilon: Extended): Boolean; overload;
  251. {$endif}
  252. {$ifdef FPC_HAS_TYPE_DOUBLE}
  253. function SameValue(const A, B: Double; Epsilon: Double): Boolean; overload;
  254. {$endif}
  255. function SameValue(const A, B: Single; Epsilon: Single): Boolean; overload;
  256. type
  257. TRoundToRange = -37..37;
  258. {$ifdef FPC_HAS_TYPE_DOUBLE}
  259. function RoundTo(const AValue: Double; const Digits: TRoundToRange): Double;
  260. {$endif}
  261. {$ifdef FPC_HAS_TYPE_EXTENDED}
  262. function RoundTo(const AVAlue: Extended; const Digits: TRoundToRange): Extended;
  263. {$endif}
  264. {$ifdef FPC_HAS_TYPE_SINGLE}
  265. function RoundTo(const AValue: Single; const Digits: TRoundToRange): Single;
  266. {$endif}
  267. {$ifdef FPC_HAS_TYPE_SINGLE}
  268. function SimpleRoundTo(const AValue: Single; const Digits: TRoundToRange = -2): Single;
  269. {$endif}
  270. {$ifdef FPC_HAS_TYPE_DOUBLE}
  271. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange = -2): Double;
  272. {$endif}
  273. {$ifdef FPC_HAS_TYPE_EXTENDED}
  274. function SimpleRoundTo(const AValue: Extended; const Digits: TRoundToRange = -2): Extended;
  275. {$endif}
  276. { angle conversion }
  277. function DegToRad(deg : float) : float;inline;
  278. function RadToDeg(rad : float) : float;inline;
  279. function GradToRad(grad : float) : float;inline;
  280. function RadToGrad(rad : float) : float;inline;
  281. function DegToGrad(deg : float) : float;inline;
  282. function GradToDeg(grad : float) : float;inline;
  283. {$ifdef FPC_HAS_TYPE_SINGLE}
  284. function CycleToDeg(const Cycles: Single): Single;
  285. {$ENDIF}
  286. {$ifdef FPC_HAS_TYPE_DOUBLE}
  287. function CycleToDeg(const Cycles: Double): Double;
  288. {$ENDIF}
  289. {$ifdef FPC_HAS_TYPE_EXTENDED}
  290. function CycleToDeg(const Cycles: Extended): Extended;
  291. {$ENDIF}
  292. {$ifdef FPC_HAS_TYPE_SINGLE}
  293. function DegToCycle(const Degrees: Single): Single;
  294. {$ENDIF}
  295. {$ifdef FPC_HAS_TYPE_DOUBLE}
  296. function DegToCycle(const Degrees: Double): Double;
  297. {$ENDIF}
  298. {$ifdef FPC_HAS_TYPE_EXTENDED}
  299. function DegToCycle(const Degrees: Extended): Extended;
  300. {$ENDIF}
  301. {$ifdef FPC_HAS_TYPE_SINGLE}
  302. function CycleToGrad(const Cycles: Single): Single;
  303. {$ENDIF}
  304. {$ifdef FPC_HAS_TYPE_DOUBLE}
  305. function CycleToGrad(const Cycles: Double): Double;
  306. {$ENDIF}
  307. {$ifdef FPC_HAS_TYPE_EXTENDED}
  308. function CycleToGrad(const Cycles: Extended): Extended;
  309. {$ENDIF}
  310. {$ifdef FPC_HAS_TYPE_SINGLE}
  311. function GradToCycle(const Grads: Single): Single;
  312. {$ENDIF}
  313. {$ifdef FPC_HAS_TYPE_DOUBLE}
  314. function GradToCycle(const Grads: Double): Double;
  315. {$ENDIF}
  316. {$ifdef FPC_HAS_TYPE_EXTENDED}
  317. function GradToCycle(const Grads: Extended): Extended;
  318. {$ENDIF}
  319. {$ifdef FPC_HAS_TYPE_SINGLE}
  320. function CycleToRad(const Cycles: Single): Single;
  321. {$ENDIF}
  322. {$ifdef FPC_HAS_TYPE_DOUBLE}
  323. function CycleToRad(const Cycles: Double): Double;
  324. {$ENDIF}
  325. {$ifdef FPC_HAS_TYPE_EXTENDED}
  326. function CycleToRad(const Cycles: Extended): Extended;
  327. {$ENDIF}
  328. {$ifdef FPC_HAS_TYPE_SINGLE}
  329. function RadToCycle(const Rads: Single): Single;
  330. {$ENDIF}
  331. {$ifdef FPC_HAS_TYPE_DOUBLE}
  332. function RadToCycle(const Rads: Double): Double;
  333. {$ENDIF}
  334. {$ifdef FPC_HAS_TYPE_EXTENDED}
  335. function RadToCycle(const Rads: Extended): Extended;
  336. {$ENDIF}
  337. {$ifdef FPC_HAS_TYPE_SINGLE}
  338. Function DegNormalize(deg : single) : single; inline;
  339. {$ENDIF}
  340. {$ifdef FPC_HAS_TYPE_DOUBLE}
  341. Function DegNormalize(deg : double) : double; inline;
  342. {$ENDIF}
  343. {$ifdef FPC_HAS_TYPE_EXTENDED}
  344. Function DegNormalize(deg : extended) : extended; inline;
  345. {$ENDIF}
  346. { trigoniometric functions }
  347. function Tan(x : float) : float;
  348. function Cotan(x : float) : float;
  349. function Cot(x : float) : float; inline;
  350. {$ifdef FPC_HAS_TYPE_SINGLE}
  351. procedure SinCos(theta : single;out sinus,cosinus : single);
  352. {$endif}
  353. {$ifdef FPC_HAS_TYPE_DOUBLE}
  354. procedure SinCos(theta : double;out sinus,cosinus : double);
  355. {$endif}
  356. {$ifdef FPC_HAS_TYPE_EXTENDED}
  357. procedure SinCos(theta : extended;out sinus,cosinus : extended);
  358. {$endif}
  359. function Secant(x : float) : float; inline;
  360. function Cosecant(x : float) : float; inline;
  361. function Sec(x : float) : float; inline;
  362. function Csc(x : float) : float; inline;
  363. { inverse functions }
  364. {$ifdef FPC_HAS_TYPE_SINGLE}
  365. function ArcCos(x : Single) : Single;
  366. {$ENDIF}
  367. {$ifdef FPC_HAS_TYPE_DOUBLE}
  368. function ArcCos(x : Double) : Double;
  369. {$ENDIF}
  370. {$ifdef FPC_HAS_TYPE_EXTENDED}
  371. function ArcCos(x : Extended) : Extended;
  372. {$ENDIF}
  373. {$ifdef FPC_HAS_TYPE_SINGLE}
  374. function ArcSin(x : Single) : Single;
  375. {$ENDIF}
  376. {$ifdef FPC_HAS_TYPE_DOUBLE}
  377. function ArcSin(x : Double) : Double;
  378. {$ENDIF}
  379. {$ifdef FPC_HAS_TYPE_EXTENDED}
  380. function ArcSin(x : Extended) : Extended;
  381. {$ENDIF}
  382. { calculates arctan(y/x) and returns an angle in the correct quadrant }
  383. function ArcTan2(y,x : float) : float;
  384. { hyperbolic functions }
  385. {$ifdef FPC_HAS_TYPE_SINGLE}
  386. function cosh(x : Single) : Single;
  387. {$ENDIF}
  388. {$ifdef FPC_HAS_TYPE_DOUBLE}
  389. function cosh(x : Double) : Double;
  390. {$ENDIF}
  391. {$ifdef FPC_HAS_TYPE_EXTENDED}
  392. function cosh(x : Extended) : Extended;
  393. {$ENDIF}
  394. {$ifdef FPC_HAS_TYPE_SINGLE}
  395. function sinh(x : Single) : Single;
  396. {$ENDIF}
  397. {$ifdef FPC_HAS_TYPE_DOUBLE}
  398. function sinh(x : Double) : Double;
  399. {$ENDIF}
  400. {$ifdef FPC_HAS_TYPE_EXTENDED}
  401. function sinh(x : Extended) : Extended;
  402. {$ENDIF}
  403. {$ifdef FPC_HAS_TYPE_SINGLE}
  404. function tanh(x : Single) : Single;
  405. {$ENDIF}
  406. {$ifdef FPC_HAS_TYPE_DOUBLE}
  407. function tanh(x : Double) : Double;
  408. {$ENDIF}
  409. {$ifdef FPC_HAS_TYPE_EXTENDED}
  410. function tanh(x : Extended) : Extended;
  411. {$ENDIF}
  412. {$ifdef FPC_HAS_TYPE_SINGLE}
  413. function SecH(const X: Single): Single;
  414. {$ENDIF}
  415. {$ifdef FPC_HAS_TYPE_DOUBLE}
  416. function SecH(const X: Double): Double;
  417. {$ENDIF}
  418. {$ifdef FPC_HAS_TYPE_EXTENDED}
  419. function SecH(const X: Extended): Extended;
  420. {$ENDIF}
  421. {$ifdef FPC_HAS_TYPE_SINGLE}
  422. function CscH(const X: Single): Single;
  423. {$ENDIF}
  424. {$ifdef FPC_HAS_TYPE_DOUBLE}
  425. function CscH(const X: Double): Double;
  426. {$ENDIF}
  427. {$ifdef FPC_HAS_TYPE_EXTENDED}
  428. function CscH(const X: Extended): Extended;
  429. {$ENDIF}
  430. {$ifdef FPC_HAS_TYPE_SINGLE}
  431. function CotH(const X: Single): Single;
  432. {$ENDIF}
  433. {$ifdef FPC_HAS_TYPE_DOUBLE}
  434. function CotH(const X: Double): Double;
  435. {$ENDIF}
  436. {$ifdef FPC_HAS_TYPE_EXTENDED}
  437. function CotH(const X: Extended): Extended;
  438. {$ENDIF}
  439. { area functions }
  440. { delphi names: }
  441. function ArcCosH(x : float) : float;inline;
  442. function ArcSinH(x : float) : float;inline;
  443. function ArcTanH(x : float) : float;inline;
  444. { IMHO the function should be called as follows (FK) }
  445. function ArCosH(x : float) : float;
  446. function ArSinH(x : float) : float;
  447. function ArTanH(x : float) : float;
  448. {$ifdef FPC_HAS_TYPE_SINGLE}
  449. function ArcSec(X: Single): Single;
  450. {$ENDIF}
  451. {$ifdef FPC_HAS_TYPE_DOUBLE}
  452. function ArcSec(X: Double): Double;
  453. {$ENDIF}
  454. {$ifdef FPC_HAS_TYPE_EXTENDED}
  455. function ArcSec(X: Extended): Extended;
  456. {$ENDIF}
  457. {$ifdef FPC_HAS_TYPE_SINGLE}
  458. function ArcCsc(X: Single): Single;
  459. {$ENDIF}
  460. {$ifdef FPC_HAS_TYPE_DOUBLE}
  461. function ArcCsc(X: Double): Double;
  462. {$ENDIF}
  463. {$ifdef FPC_HAS_TYPE_EXTENDED}
  464. function ArcCsc(X: Extended): Extended;
  465. {$ENDIF}
  466. {$ifdef FPC_HAS_TYPE_SINGLE}
  467. function ArcCot(X: Single): Single;
  468. {$ENDIF}
  469. {$ifdef FPC_HAS_TYPE_DOUBLE}
  470. function ArcCot(X: Double): Double;
  471. {$ENDIF}
  472. {$ifdef FPC_HAS_TYPE_EXTENDED}
  473. function ArcCot(X: Extended): Extended;
  474. {$ENDIF}
  475. {$ifdef FPC_HAS_TYPE_SINGLE}
  476. function ArcSecH(X : Single): Single;
  477. {$ENDIF}
  478. {$ifdef FPC_HAS_TYPE_DOUBLE}
  479. function ArcSecH(X : Double): Double;
  480. {$ENDIF}
  481. {$ifdef FPC_HAS_TYPE_EXTENDED}
  482. function ArcSecH(X : Extended): Extended;
  483. {$ENDIF}
  484. {$ifdef FPC_HAS_TYPE_SINGLE}
  485. function ArcCscH(X: Single): Single;
  486. {$ENDIF}
  487. {$ifdef FPC_HAS_TYPE_DOUBLE}
  488. function ArcCscH(X: Double): Double;
  489. {$ENDIF}
  490. {$ifdef FPC_HAS_TYPE_EXTENDED}
  491. function ArcCscH(X: Extended): Extended;
  492. {$ENDIF}
  493. {$ifdef FPC_HAS_TYPE_SINGLE}
  494. function ArcCotH(X: Single): Single;
  495. {$ENDIF}
  496. {$ifdef FPC_HAS_TYPE_DOUBLE}
  497. function ArcCotH(X: Double): Double;
  498. {$ENDIF}
  499. {$ifdef FPC_HAS_TYPE_EXTENDED}
  500. function ArcCotH(X: Extended): Extended;
  501. {$ENDIF}
  502. { triangle functions }
  503. { returns the length of the hypotenuse of a right triangle }
  504. { if x and y are the other sides }
  505. function Hypot(x,y : float) : float;
  506. { logarithm functions }
  507. function Log10(x : float) : float;
  508. function Log2(x : float) : float;
  509. function LogN(n,x : float) : float;
  510. { returns natural logarithm of x+1, accurate for x values near zero }
  511. function LnXP1(x : float) : float;
  512. { exponential functions }
  513. function Power(base,exponent : float) : float;
  514. { base^exponent }
  515. function IntPower(base : float;exponent : longint) : float;
  516. operator ** (base,exponent : float) e: float; inline;
  517. operator ** (base,exponent : int64) res: int64;
  518. { number converting }
  519. { rounds x towards positive infinity }
  520. function Ceil(x : float) : Integer;
  521. function Ceil64(x: float): Int64;
  522. { rounds x towards negative infinity }
  523. function Floor(x : float) : Integer;
  524. function Floor64(x: float): Int64;
  525. { misc. functions }
  526. {$ifdef FPC_HAS_TYPE_SINGLE}
  527. { splits x into mantissa and exponent (to base 2) }
  528. procedure Frexp(X: single; out Mantissa: single; out Exponent: integer);
  529. { returns x*(2^p) }
  530. function Ldexp(X: single; p: Integer) : single;
  531. {$endif}
  532. {$ifdef FPC_HAS_TYPE_DOUBLE}
  533. procedure Frexp(X: double; out Mantissa: double; out Exponent: integer);
  534. function Ldexp(X: double; p: Integer) : double;
  535. {$endif}
  536. {$ifdef FPC_HAS_TYPE_EXTENDED}
  537. procedure Frexp(X: extended; out Mantissa: extended; out Exponent: integer);
  538. function Ldexp(X: extended; p: Integer) : extended;
  539. {$endif}
  540. { statistical functions }
  541. {$ifdef FPC_HAS_TYPE_SINGLE}
  542. function Mean(const data : array of Single) : float;
  543. function Sum(const data : array of Single) : float;inline;
  544. function Mean(const data : PSingle; Const N : longint) : float;
  545. function Sum(const data : PSingle; Const N : Longint) : float;
  546. {$endif FPC_HAS_TYPE_SINGLE}
  547. {$ifdef FPC_HAS_TYPE_DOUBLE}
  548. function Mean(const data : array of double) : float;inline;
  549. function Sum(const data : array of double) : float;inline;
  550. function Mean(const data : PDouble; Const N : longint) : float;
  551. function Sum(const data : PDouble; Const N : Longint) : float;
  552. {$endif FPC_HAS_TYPE_DOUBLE}
  553. {$ifdef FPC_HAS_TYPE_EXTENDED}
  554. function Mean(const data : array of Extended) : float;
  555. function Sum(const data : array of Extended) : float;inline;
  556. function Mean(const data : PExtended; Const N : longint) : float;
  557. function Sum(const data : PExtended; Const N : Longint) : float;
  558. {$endif FPC_HAS_TYPE_EXTENDED}
  559. function SumInt(const data : PInt64;Const N : longint) : Int64;
  560. function SumInt(const data : array of Int64) : Int64;inline;
  561. function Mean(const data : PInt64; const N : Longint):Float;
  562. function Mean(const data: array of Int64):Float;
  563. function SumInt(const data : PInteger; Const N : longint) : Int64;
  564. function SumInt(const data : array of Integer) : Int64;inline;
  565. function Mean(const data : PInteger; const N : Longint):Float;
  566. function Mean(const data: array of Integer):Float;
  567. {$ifdef FPC_HAS_TYPE_SINGLE}
  568. function SumOfSquares(const data : array of Single) : float;inline;
  569. function SumOfSquares(const data : PSingle; Const N : Integer) : float;
  570. { calculates the sum and the sum of squares of data }
  571. procedure SumsAndSquares(const data : array of Single;
  572. var sum,sumofsquares : float);inline;
  573. procedure SumsAndSquares(const data : PSingle; Const N : Integer;
  574. var sum,sumofsquares : float);
  575. {$endif FPC_HAS_TYPE_SINGLE}
  576. {$ifdef FPC_HAS_TYPE_DOUBLE}
  577. function SumOfSquares(const data : array of double) : float;inline;
  578. function SumOfSquares(const data : PDouble; Const N : Integer) : float;
  579. { calculates the sum and the sum of squares of data }
  580. procedure SumsAndSquares(const data : array of Double;
  581. var sum,sumofsquares : float);inline;
  582. procedure SumsAndSquares(const data : PDouble; Const N : Integer;
  583. var sum,sumofsquares : float);
  584. {$endif FPC_HAS_TYPE_DOUBLE}
  585. {$ifdef FPC_HAS_TYPE_EXTENDED}
  586. function SumOfSquares(const data : array of Extended) : float;inline;
  587. function SumOfSquares(const data : PExtended; Const N : Integer) : float;
  588. { calculates the sum and the sum of squares of data }
  589. procedure SumsAndSquares(const data : array of Extended;
  590. var sum,sumofsquares : float);inline;
  591. procedure SumsAndSquares(const data : PExtended; Const N : Integer;
  592. var sum,sumofsquares : float);
  593. {$endif FPC_HAS_TYPE_EXTENDED}
  594. {$ifdef FPC_HAS_TYPE_SINGLE}
  595. function MinValue(const data : array of Single) : Single;inline;
  596. function MinValue(const data : PSingle; Const N : Integer) : Single;
  597. function MaxValue(const data : array of Single) : Single;inline;
  598. function MaxValue(const data : PSingle; Const N : Integer) : Single;
  599. {$endif FPC_HAS_TYPE_SINGLE}
  600. {$ifdef FPC_HAS_TYPE_DOUBLE}
  601. function MinValue(const data : array of Double) : Double;inline;
  602. function MinValue(const data : PDouble; Const N : Integer) : Double;
  603. function MaxValue(const data : array of Double) : Double;inline;
  604. function MaxValue(const data : PDouble; Const N : Integer) : Double;
  605. {$endif FPC_HAS_TYPE_DOUBLE}
  606. {$ifdef FPC_HAS_TYPE_EXTENDED}
  607. function MinValue(const data : array of Extended) : Extended;inline;
  608. function MinValue(const data : PExtended; Const N : Integer) : Extended;
  609. function MaxValue(const data : array of Extended) : Extended;inline;
  610. function MaxValue(const data : PExtended; Const N : Integer) : Extended;
  611. {$endif FPC_HAS_TYPE_EXTENDED}
  612. function MinValue(const data : array of integer) : Integer;inline;
  613. function MinValue(const Data : PInteger; Const N : Integer): Integer;
  614. function MaxValue(const data : array of integer) : Integer;inline;
  615. function MaxValue(const data : PInteger; Const N : Integer) : Integer;
  616. { returns random values with gaussian distribution }
  617. function RandG(mean,stddev : float) : float;
  618. function RandomRange(const aFrom, aTo: Integer): Integer;
  619. function RandomRange(const aFrom, aTo: Int64): Int64;
  620. {$ifdef FPC_HAS_TYPE_SINGLE}
  621. { calculates the standard deviation }
  622. function StdDev(const data : array of Single) : float;inline;
  623. function StdDev(const data : PSingle; Const N : Integer) : float;
  624. { calculates the mean and stddev }
  625. procedure MeanAndStdDev(const data : array of Single;
  626. var mean,stddev : float);inline;
  627. procedure MeanAndStdDev(const data : PSingle;
  628. Const N : Longint;var mean,stddev : float);
  629. function Variance(const data : array of Single) : float;inline;
  630. function TotalVariance(const data : array of Single) : float;inline;
  631. function Variance(const data : PSingle; Const N : Integer) : float;
  632. function TotalVariance(const data : PSingle; Const N : Integer) : float;
  633. { Population (aka uncorrected) variance and standard deviation }
  634. function PopnStdDev(const data : array of Single) : float;inline;
  635. function PopnStdDev(const data : PSingle; Const N : Integer) : float;
  636. function PopnVariance(const data : PSingle; Const N : Integer) : float;
  637. function PopnVariance(const data : array of Single) : float;inline;
  638. procedure MomentSkewKurtosis(const data : array of Single;
  639. out m1,m2,m3,m4,skew,kurtosis : float);inline;
  640. procedure MomentSkewKurtosis(const data : PSingle; Const N : Integer;
  641. out m1,m2,m3,m4,skew,kurtosis : float);
  642. { geometrical function }
  643. { returns the euclidean L2 norm }
  644. function Norm(const data : array of Single) : float;inline;
  645. function Norm(const data : PSingle; Const N : Integer) : float;
  646. {$endif FPC_HAS_TYPE_SINGLE}
  647. {$ifdef FPC_HAS_TYPE_DOUBLE}
  648. { calculates the standard deviation }
  649. function StdDev(const data : array of Double) : float;inline;
  650. function StdDev(const data : PDouble; Const N : Integer) : float;
  651. { calculates the mean and stddev }
  652. procedure MeanAndStdDev(const data : array of Double;
  653. var mean,stddev : float);inline;
  654. procedure MeanAndStdDev(const data : PDouble;
  655. Const N : Longint;var mean,stddev : float);
  656. function Variance(const data : array of Double) : float;inline;
  657. function TotalVariance(const data : array of Double) : float;inline;
  658. function Variance(const data : PDouble; Const N : Integer) : float;
  659. function TotalVariance(const data : PDouble; Const N : Integer) : float;
  660. { Population (aka uncorrected) variance and standard deviation }
  661. function PopnStdDev(const data : array of Double) : float;inline;
  662. function PopnStdDev(const data : PDouble; Const N : Integer) : float;
  663. function PopnVariance(const data : PDouble; Const N : Integer) : float;
  664. function PopnVariance(const data : array of Double) : float;inline;
  665. procedure MomentSkewKurtosis(const data : array of Double;
  666. out m1,m2,m3,m4,skew,kurtosis : float);inline;
  667. procedure MomentSkewKurtosis(const data : PDouble; Const N : Integer;
  668. out m1,m2,m3,m4,skew,kurtosis : float);
  669. { geometrical function }
  670. { returns the euclidean L2 norm }
  671. function Norm(const data : array of double) : float;inline;
  672. function Norm(const data : PDouble; Const N : Integer) : float;
  673. {$endif FPC_HAS_TYPE_DOUBLE}
  674. {$ifdef FPC_HAS_TYPE_EXTENDED}
  675. { calculates the standard deviation }
  676. function StdDev(const data : array of Extended) : float;inline;
  677. function StdDev(const data : PExtended; Const N : Integer) : float;
  678. { calculates the mean and stddev }
  679. procedure MeanAndStdDev(const data : array of Extended;
  680. var mean,stddev : float);inline;
  681. procedure MeanAndStdDev(const data : PExtended;
  682. Const N : Longint;var mean,stddev : float);
  683. function Variance(const data : array of Extended) : float;inline;
  684. function TotalVariance(const data : array of Extended) : float;inline;
  685. function Variance(const data : PExtended; Const N : Integer) : float;
  686. function TotalVariance(const data : PExtended; Const N : Integer) : float;
  687. { Population (aka uncorrected) variance and standard deviation }
  688. function PopnStdDev(const data : array of Extended) : float;inline;
  689. function PopnStdDev(const data : PExtended; Const N : Integer) : float;
  690. function PopnVariance(const data : PExtended; Const N : Integer) : float;
  691. function PopnVariance(const data : array of Extended) : float;inline;
  692. procedure MomentSkewKurtosis(const data : array of Extended;
  693. out m1,m2,m3,m4,skew,kurtosis : float);inline;
  694. procedure MomentSkewKurtosis(const data : PExtended; Const N : Integer;
  695. out m1,m2,m3,m4,skew,kurtosis : float);
  696. { geometrical function }
  697. { returns the euclidean L2 norm }
  698. function Norm(const data : array of Extended) : float;inline;
  699. function Norm(const data : PExtended; Const N : Integer) : float;
  700. {$endif FPC_HAS_TYPE_EXTENDED}
  701. { Financial functions }
  702. function FutureValue(ARate: Float; NPeriods: Integer;
  703. APayment, APresentValue: Float; APaymentTime: TPaymentTime): Float;
  704. function InterestRate(NPeriods: Integer; APayment, APresentValue, AFutureValue: Float;
  705. APaymentTime: TPaymentTime): Float;
  706. function NumberOfPeriods(ARate, APayment, APresentValue, AFutureValue: Float;
  707. APaymentTime: TPaymentTime): Float;
  708. function Payment(ARate: Float; NPeriods: Integer;
  709. APresentValue, AFutureValue: Float; APaymentTime: TPaymentTime): Float;
  710. function PresentValue(ARate: Float; NPeriods: Integer;
  711. APayment, AFutureValue: Float; APaymentTime: TPaymentTime): Float;
  712. { Misc functions }
  713. function IfThen(val:boolean;const iftrue:integer; const iffalse:integer= 0) :integer; inline; overload;
  714. function IfThen(val:boolean;const iftrue:int64 ; const iffalse:int64 = 0) :int64; inline; overload;
  715. function IfThen(val:boolean;const iftrue:double ; const iffalse:double =0.0):double; inline; overload;
  716. function CompareValue ( const A, B : Integer) : TValueRelationship; inline;
  717. function CompareValue ( const A, B : Int64) : TValueRelationship; inline;
  718. function CompareValue ( const A, B : QWord) : TValueRelationship; inline;
  719. {$ifdef FPC_HAS_TYPE_SINGLE}
  720. function CompareValue ( const A, B : Single; delta : Single = 0.0 ) : TValueRelationship; inline;
  721. {$endif}
  722. {$ifdef FPC_HAS_TYPE_DOUBLE}
  723. function CompareValue ( const A, B : Double; delta : Double = 0.0) : TValueRelationship; inline;
  724. {$endif}
  725. {$ifdef FPC_HAS_TYPE_EXTENDED}
  726. function CompareValue ( const A, B : Extended; delta : Extended = 0.0 ) : TValueRelationship; inline;
  727. {$endif}
  728. function RandomFrom(const AValues: array of Double): Double; overload;
  729. function RandomFrom(const AValues: array of Integer): Integer; overload;
  730. function RandomFrom(const AValues: array of Int64): Int64; overload;
  731. {$if FPC_FULLVERSION >=30101}
  732. generic function RandomFrom<T>(const AValues:array of T):T;
  733. {$endif}
  734. { cpu specific stuff }
  735. type
  736. TFPURoundingMode = system.TFPURoundingMode;
  737. TFPUPrecisionMode = system.TFPUPrecisionMode;
  738. TFPUException = system.TFPUException;
  739. TFPUExceptionMask = system.TFPUExceptionMask;
  740. function GetRoundMode: TFPURoundingMode;
  741. function SetRoundMode(const RoundMode: TFPURoundingMode): TFPURoundingMode;
  742. function GetPrecisionMode: TFPUPrecisionMode;
  743. function SetPrecisionMode(const Precision: TFPUPrecisionMode): TFPUPrecisionMode;
  744. function GetExceptionMask: TFPUExceptionMask;
  745. function SetExceptionMask(const Mask: TFPUExceptionMask): TFPUExceptionMask;
  746. procedure ClearExceptions(RaisePending: Boolean =true);
  747. implementation
  748. function copysign(x,y: float): float; forward; { returns abs(x)*sign(y) }
  749. { include cpu specific stuff }
  750. {$i mathu.inc}
  751. ResourceString
  752. SMathError = 'Math Error : %s';
  753. SInvalidArgument = 'Invalid argument';
  754. Procedure DoMathError(Const S : String);
  755. begin
  756. Raise EMathError.CreateFmt(SMathError,[S]);
  757. end;
  758. Procedure InvalidArgument;
  759. begin
  760. Raise EInvalidArgument.Create(SInvalidArgument);
  761. end;
  762. function Sign(const AValue: Integer): TValueSign;inline;
  763. begin
  764. result:=TValueSign(
  765. SarLongint(AValue,sizeof(AValue)*8-1) or { gives -1 for negative values, 0 otherwise }
  766. (longint(-AValue) shr (sizeof(AValue)*8-1)) { gives 1 for positive values, 0 otherwise }
  767. );
  768. end;
  769. function Sign(const AValue: Int64): TValueSign;inline;
  770. begin
  771. {$ifdef cpu64}
  772. result:=TValueSign(
  773. SarInt64(AValue,sizeof(AValue)*8-1) or
  774. (-AValue shr (sizeof(AValue)*8-1))
  775. );
  776. {$else cpu64}
  777. If Avalue<0 then
  778. Result:=NegativeValue
  779. else If Avalue>0 then
  780. Result:=PositiveValue
  781. else
  782. Result:=ZeroValue;
  783. {$endif}
  784. end;
  785. {$ifdef FPC_HAS_TYPE_SINGLE}
  786. function Sign(const AValue: Single): TValueSign;inline;
  787. begin
  788. Result:=ord(AValue>0.0)-ord(AValue<0.0);
  789. end;
  790. {$endif}
  791. function Sign(const AValue: Double): TValueSign;inline;
  792. begin
  793. Result:=ord(AValue>0.0)-ord(AValue<0.0);
  794. end;
  795. {$ifdef FPC_HAS_TYPE_EXTENDED}
  796. function Sign(const AValue: Extended): TValueSign;inline;
  797. begin
  798. Result:=ord(AValue>0.0)-ord(AValue<0.0);
  799. end;
  800. {$endif}
  801. function degtorad(deg : float) : float;inline;
  802. begin
  803. degtorad:=deg*(pi/180.0);
  804. end;
  805. function radtodeg(rad : float) : float;inline;
  806. begin
  807. radtodeg:=rad*(180.0/pi);
  808. end;
  809. function gradtorad(grad : float) : float;inline;
  810. begin
  811. gradtorad:=grad*(pi/200.0);
  812. end;
  813. function radtograd(rad : float) : float;inline;
  814. begin
  815. radtograd:=rad*(200.0/pi);
  816. end;
  817. function degtograd(deg : float) : float;inline;
  818. begin
  819. degtograd:=deg*(200.0/180.0);
  820. end;
  821. function gradtodeg(grad : float) : float;inline;
  822. begin
  823. gradtodeg:=grad*(180.0/200.0);
  824. end;
  825. {$ifdef FPC_HAS_TYPE_SINGLE}
  826. function CycleToDeg(const Cycles: Single): Single;
  827. begin
  828. CycleToDeg:=Cycles*360.0;
  829. end;
  830. {$ENDIF}
  831. {$ifdef FPC_HAS_TYPE_DOUBLE}
  832. function CycleToDeg(const Cycles: Double): Double;
  833. begin
  834. CycleToDeg:=Cycles*360.0;
  835. end;
  836. {$ENDIF}
  837. {$ifdef FPC_HAS_TYPE_EXTENDED}
  838. function CycleToDeg(const Cycles: Extended): Extended;
  839. begin
  840. CycleToDeg:=Cycles*360.0;
  841. end;
  842. {$ENDIF}
  843. {$ifdef FPC_HAS_TYPE_SINGLE}
  844. function DegToCycle(const Degrees: Single): Single;
  845. begin
  846. DegToCycle:=Degrees*(1/360.0);
  847. end;
  848. {$ENDIF}
  849. {$ifdef FPC_HAS_TYPE_DOUBLE}
  850. function DegToCycle(const Degrees: Double): Double;
  851. begin
  852. DegToCycle:=Degrees*(1/360.0);
  853. end;
  854. {$ENDIF}
  855. {$ifdef FPC_HAS_TYPE_EXTENDED}
  856. function DegToCycle(const Degrees: Extended): Extended;
  857. begin
  858. DegToCycle:=Degrees*(1/360.0);
  859. end;
  860. {$ENDIF}
  861. {$ifdef FPC_HAS_TYPE_SINGLE}
  862. function CycleToGrad(const Cycles: Single): Single;
  863. begin
  864. CycleToGrad:=Cycles*400.0;
  865. end;
  866. {$ENDIF}
  867. {$ifdef FPC_HAS_TYPE_DOUBLE}
  868. function CycleToGrad(const Cycles: Double): Double;
  869. begin
  870. CycleToGrad:=Cycles*400.0;
  871. end;
  872. {$ENDIF}
  873. {$ifdef FPC_HAS_TYPE_EXTENDED}
  874. function CycleToGrad(const Cycles: Extended): Extended;
  875. begin
  876. CycleToGrad:=Cycles*400.0;
  877. end;
  878. {$ENDIF}
  879. {$ifdef FPC_HAS_TYPE_SINGLE}
  880. function GradToCycle(const Grads: Single): Single;
  881. begin
  882. GradToCycle:=Grads*(1/400.0);
  883. end;
  884. {$ENDIF}
  885. {$ifdef FPC_HAS_TYPE_DOUBLE}
  886. function GradToCycle(const Grads: Double): Double;
  887. begin
  888. GradToCycle:=Grads*(1/400.0);
  889. end;
  890. {$ENDIF}
  891. {$ifdef FPC_HAS_TYPE_EXTENDED}
  892. function GradToCycle(const Grads: Extended): Extended;
  893. begin
  894. GradToCycle:=Grads*(1/400.0);
  895. end;
  896. {$ENDIF}
  897. {$ifdef FPC_HAS_TYPE_SINGLE}
  898. function CycleToRad(const Cycles: Single): Single;
  899. begin
  900. CycleToRad:=Cycles*2*pi;
  901. end;
  902. {$ENDIF}
  903. {$ifdef FPC_HAS_TYPE_DOUBLE}
  904. function CycleToRad(const Cycles: Double): Double;
  905. begin
  906. CycleToRad:=Cycles*2*pi;
  907. end;
  908. {$ENDIF}
  909. {$ifdef FPC_HAS_TYPE_EXTENDED}
  910. function CycleToRad(const Cycles: Extended): Extended;
  911. begin
  912. CycleToRad:=Cycles*2*pi;
  913. end;
  914. {$ENDIF}
  915. {$ifdef FPC_HAS_TYPE_SINGLE}
  916. function RadToCycle(const Rads: Single): Single;
  917. begin
  918. RadToCycle:=Rads*(1/(2*pi));
  919. end;
  920. {$ENDIF}
  921. {$ifdef FPC_HAS_TYPE_DOUBLE}
  922. function RadToCycle(const Rads: Double): Double;
  923. begin
  924. RadToCycle:=Rads*(1/(2*pi));
  925. end;
  926. {$ENDIF}
  927. {$ifdef FPC_HAS_TYPE_EXTENDED}
  928. function RadToCycle(const Rads: Extended): Extended;
  929. begin
  930. RadToCycle:=Rads*(1/(2*pi));
  931. end;
  932. {$ENDIF}
  933. {$ifdef FPC_HAS_TYPE_SINGLE}
  934. Function DegNormalize(deg : single) : single;
  935. begin
  936. Result:=Deg-Int(Deg/360)*360;
  937. If Result<0 then Result:=Result+360;
  938. end;
  939. {$ENDIF}
  940. {$ifdef FPC_HAS_TYPE_DOUBLE}
  941. Function DegNormalize(deg : double) : double; inline;
  942. begin
  943. Result:=Deg-Int(Deg/360)*360;
  944. If (Result<0) then Result:=Result+360;
  945. end;
  946. {$ENDIF}
  947. {$ifdef FPC_HAS_TYPE_EXTENDED}
  948. Function DegNormalize(deg : extended) : extended; inline;
  949. begin
  950. Result:=Deg-Int(Deg/360)*360;
  951. If Result<0 then Result:=Result+360;
  952. end;
  953. {$ENDIF}
  954. {$ifndef FPC_MATH_HAS_TAN}
  955. function tan(x : float) : float;
  956. var
  957. _sin,_cos : float;
  958. begin
  959. sincos(x,_sin,_cos);
  960. tan:=_sin/_cos;
  961. end;
  962. {$endif FPC_MATH_HAS_TAN}
  963. {$ifndef FPC_MATH_HAS_COTAN}
  964. function cotan(x : float) : float;
  965. var
  966. _sin,_cos : float;
  967. begin
  968. sincos(x,_sin,_cos);
  969. cotan:=_cos/_sin;
  970. end;
  971. {$endif FPC_MATH_HAS_COTAN}
  972. function cot(x : float) : float; inline;
  973. begin
  974. cot := cotan(x);
  975. end;
  976. {$ifndef FPC_MATH_HAS_SINCOS}
  977. {$ifdef FPC_HAS_TYPE_SINGLE}
  978. procedure sincos(theta : single;out sinus,cosinus : single);
  979. begin
  980. sinus:=sin(theta);
  981. cosinus:=cos(theta);
  982. end;
  983. {$endif}
  984. {$ifdef FPC_HAS_TYPE_DOUBLE}
  985. procedure sincos(theta : double;out sinus,cosinus : double);
  986. begin
  987. sinus:=sin(theta);
  988. cosinus:=cos(theta);
  989. end;
  990. {$endif}
  991. {$ifdef FPC_HAS_TYPE_EXTENDED}
  992. procedure sincos(theta : extended;out sinus,cosinus : extended);
  993. begin
  994. sinus:=sin(theta);
  995. cosinus:=cos(theta);
  996. end;
  997. {$endif}
  998. {$endif FPC_MATH_HAS_SINCOS}
  999. function secant(x : float) : float; inline;
  1000. begin
  1001. secant := 1 / cos(x);
  1002. end;
  1003. function cosecant(x : float) : float; inline;
  1004. begin
  1005. cosecant := 1 / sin(x);
  1006. end;
  1007. function sec(x : float) : float; inline;
  1008. begin
  1009. sec := secant(x);
  1010. end;
  1011. function csc(x : float) : float; inline;
  1012. begin
  1013. csc := cosecant(x);
  1014. end;
  1015. { arcsin and arccos functions from AMath library (C) Copyright 2009-2013 Wolfgang Ehrhardt }
  1016. {$ifdef FPC_HAS_TYPE_SINGLE}
  1017. function arcsin(x : Single) : Single;
  1018. begin
  1019. arcsin:=arctan2(x,sqrt((1.0-x)*(1.0+x)));
  1020. end;
  1021. {$ENDIF}
  1022. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1023. function arcsin(x : Double) : Double;
  1024. begin
  1025. arcsin:=arctan2(x,sqrt((1.0-x)*(1.0+x)));
  1026. end;
  1027. {$ENDIF}
  1028. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1029. function arcsin(x : Extended) : Extended;
  1030. begin
  1031. arcsin:=arctan2(x,sqrt((1.0-x)*(1.0+x)));
  1032. end;
  1033. {$ENDIF}
  1034. {$ifdef FPC_HAS_TYPE_SINGLE}
  1035. function Arccos(x : Single) : Single;
  1036. begin
  1037. arccos:=arctan2(sqrt((1.0-x)*(1.0+x)),x);
  1038. end;
  1039. {$ENDIF}
  1040. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1041. function Arccos(x : Double) : Double;
  1042. begin
  1043. arccos:=arctan2(sqrt((1.0-x)*(1.0+x)),x);
  1044. end;
  1045. {$ENDIF}
  1046. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1047. function Arccos(x : Extended) : Extended;
  1048. begin
  1049. arccos:=arctan2(sqrt((1.0-x)*(1.0+x)),x);
  1050. end;
  1051. {$ENDIF}
  1052. {$ifndef FPC_MATH_HAS_ARCTAN2}
  1053. function arctan2(y,x : float) : float;
  1054. begin
  1055. if x=0 then
  1056. begin
  1057. if y=0 then
  1058. result:=0.0
  1059. else if y>0 then
  1060. result:=pi/2
  1061. else
  1062. result:=-pi/2;
  1063. end
  1064. else
  1065. begin
  1066. result:=ArcTan(y/x);
  1067. if x<0 then
  1068. if y<0 then
  1069. result:=result-pi
  1070. else
  1071. result:=result+pi;
  1072. end;
  1073. end;
  1074. {$endif FPC_MATH_HAS_ARCTAN2}
  1075. {$ifdef FPC_HAS_TYPE_SINGLE}
  1076. function cosh(x : Single) : Single;
  1077. var
  1078. temp : ValReal;
  1079. begin
  1080. temp:=exp(x);
  1081. cosh:=0.5*(temp+1.0/temp);
  1082. end;
  1083. {$ENDIF}
  1084. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1085. function cosh(x : Double) : Double;
  1086. var
  1087. temp : ValReal;
  1088. begin
  1089. temp:=exp(x);
  1090. cosh:=0.5*(temp+1.0/temp);
  1091. end;
  1092. {$ENDIF}
  1093. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1094. function cosh(x : Extended) : Extended;
  1095. var
  1096. temp : Extended;
  1097. begin
  1098. temp:=exp(x);
  1099. cosh:=0.5*(temp+1.0/temp);
  1100. end;
  1101. {$ENDIF}
  1102. {$ifdef FPC_HAS_TYPE_SINGLE}
  1103. function sinh(x : Single) : Single;
  1104. var
  1105. temp : ValReal;
  1106. begin
  1107. temp:=exp(x);
  1108. { gives better behavior around zero, and in particular ensures that sinh(-0.0)=-0.0 }
  1109. if temp=1 then
  1110. exit(x);
  1111. sinh:=0.5*(temp-1.0/temp);
  1112. end;
  1113. {$ENDIF}
  1114. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1115. function sinh(x : Double) : Double;
  1116. var
  1117. temp : ValReal;
  1118. begin
  1119. temp:=exp(x);
  1120. if temp=1 then
  1121. exit(x);
  1122. sinh:=0.5*(temp-1.0/temp);
  1123. end;
  1124. {$ENDIF}
  1125. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1126. function sinh(x : Extended) : Extended;
  1127. var
  1128. temp : Extended;
  1129. begin
  1130. temp:=exp(x);
  1131. if temp=1 then
  1132. exit(x);
  1133. sinh:=0.5*(temp-1.0/temp);
  1134. end;
  1135. {$ENDIF}
  1136. {$ifdef FPC_HAS_TYPE_SINGLE}
  1137. function tanh(x : Single) : Single;
  1138. var
  1139. tmp:ValReal;
  1140. begin
  1141. if x < 0 then begin
  1142. tmp:=exp(2*x);
  1143. if tmp=1 then
  1144. exit(x);
  1145. result:=(tmp-1)/(1+tmp)
  1146. end
  1147. else begin
  1148. tmp:=exp(-2*x);
  1149. if tmp=1 then
  1150. exit(x);
  1151. result:=(1-tmp)/(1+tmp)
  1152. end;
  1153. end;
  1154. {$ENDIF}
  1155. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1156. function tanh(x : Double) : Double;
  1157. var
  1158. tmp:ValReal;
  1159. begin
  1160. if x < 0 then begin
  1161. tmp:=exp(2*x);
  1162. if tmp=1 then
  1163. exit(x);
  1164. result:=(tmp-1)/(1+tmp)
  1165. end
  1166. else begin
  1167. tmp:=exp(-2*x);
  1168. if tmp=1 then
  1169. exit(x);
  1170. result:=(1-tmp)/(1+tmp)
  1171. end;
  1172. end;
  1173. {$ENDIF}
  1174. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1175. function tanh(x : Extended) : Extended;
  1176. var
  1177. tmp:Extended;
  1178. begin
  1179. if x < 0 then begin
  1180. tmp:=exp(2*x);
  1181. if tmp=1 then
  1182. exit(x);
  1183. result:=(tmp-1)/(1+tmp)
  1184. end
  1185. else begin
  1186. tmp:=exp(-2*x);
  1187. if tmp=1 then
  1188. exit(x);
  1189. result:=(1-tmp)/(1+tmp)
  1190. end;
  1191. end;
  1192. {$ENDIF}
  1193. {$ifdef FPC_HAS_TYPE_SINGLE}
  1194. function SecH(const X: Single): Single;
  1195. var
  1196. Ex: ValReal;
  1197. begin
  1198. //https://en.wikipedia.org/wiki/Hyperbolic_functions#Definitions
  1199. //SecH = 2 / (e^X + e^-X)
  1200. Ex:=Exp(X);
  1201. SecH:=2/(Ex+1/Ex);
  1202. end;
  1203. {$ENDIF}
  1204. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1205. function SecH(const X: Double): Double;
  1206. var
  1207. Ex: ValReal;
  1208. begin
  1209. Ex:=Exp(X);
  1210. SecH:=2/(Ex+1/Ex);
  1211. end;
  1212. {$ENDIF}
  1213. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1214. function SecH(const X: Extended): Extended;
  1215. var
  1216. Ex: Extended;
  1217. begin
  1218. Ex:=Exp(X);
  1219. SecH:=2/(Ex+1/Ex);
  1220. end;
  1221. {$ENDIF}
  1222. {$ifdef FPC_HAS_TYPE_SINGLE}
  1223. function CscH(const X: Single): Single;
  1224. var
  1225. Ex: ValReal;
  1226. begin
  1227. //CscH = 2 / (e^X - e^-X)
  1228. Ex:=Exp(X);
  1229. CscH:=2/(Ex-1/Ex);
  1230. end;
  1231. {$ENDIF}
  1232. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1233. function CscH(const X: Double): Double;
  1234. var
  1235. Ex: ValReal;
  1236. begin
  1237. Ex:=Exp(X);
  1238. CscH:=2/(Ex-1/Ex);
  1239. end;
  1240. {$ENDIF}
  1241. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1242. function CscH(const X: Extended): Extended;
  1243. var
  1244. Ex: Extended;
  1245. begin
  1246. Ex:=Exp(X);
  1247. CscH:=2/(Ex-1/Ex);
  1248. end;
  1249. {$ENDIF}
  1250. {$ifdef FPC_HAS_TYPE_SINGLE}
  1251. function CotH(const X: Single): Single;
  1252. var
  1253. e2: ValReal;
  1254. begin
  1255. if x < 0 then begin
  1256. e2:=exp(2*x);
  1257. if e2=1 then
  1258. exit(1/x);
  1259. result:=(1+e2)/(e2-1)
  1260. end
  1261. else begin
  1262. e2:=exp(-2*x);
  1263. if e2=1 then
  1264. exit(1/x);
  1265. result:=(1+e2)/(1-e2)
  1266. end;
  1267. end;
  1268. {$ENDIF}
  1269. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1270. function CotH(const X: Double): Double;
  1271. var
  1272. e2: ValReal;
  1273. begin
  1274. if x < 0 then begin
  1275. e2:=exp(2*x);
  1276. if e2=1 then
  1277. exit(1/x);
  1278. result:=(1+e2)/(e2-1)
  1279. end
  1280. else begin
  1281. e2:=exp(-2*x);
  1282. if e2=1 then
  1283. exit(1/x);
  1284. result:=(1+e2)/(1-e2)
  1285. end;
  1286. end;
  1287. {$ENDIF}
  1288. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1289. function CotH(const X: Extended): Extended;
  1290. var
  1291. e2: Extended;
  1292. begin
  1293. if x < 0 then begin
  1294. e2:=exp(2*x);
  1295. if e2=1 then
  1296. exit(1/x);
  1297. result:=(1+e2)/(e2-1)
  1298. end
  1299. else begin
  1300. e2:=exp(-2*x);
  1301. if e2=1 then
  1302. exit(1/x);
  1303. result:=(1+e2)/(1-e2)
  1304. end;
  1305. end;
  1306. {$ENDIF}
  1307. function arccosh(x : float) : float; inline;
  1308. begin
  1309. arccosh:=arcosh(x);
  1310. end;
  1311. function arcsinh(x : float) : float;inline;
  1312. begin
  1313. arcsinh:=arsinh(x);
  1314. end;
  1315. function arctanh(x : float) : float;inline;
  1316. begin
  1317. arctanh:=artanh(x);
  1318. end;
  1319. function arcosh(x : float) : float;
  1320. begin
  1321. { Provides accuracy about 4*eps near 1.0 }
  1322. arcosh:=Ln(x+Sqrt((x-1.0)*(x+1.0)));
  1323. end;
  1324. function arsinh(x : float) : float;
  1325. var
  1326. z: float;
  1327. begin
  1328. z:=abs(x);
  1329. z:=Ln(z+Sqrt(1+z*z));
  1330. { copysign ensures that arsinh(-Inf)=-Inf and arsinh(-0.0)=-0.0 }
  1331. arsinh:=copysign(z,x);
  1332. end;
  1333. function artanh(x : float) : float;
  1334. begin
  1335. artanh:=(lnxp1(x)-lnxp1(-x))*0.5;
  1336. end;
  1337. {$ifdef FPC_HAS_TYPE_SINGLE}
  1338. function ArcSec(X: Single): Single;
  1339. begin
  1340. ArcSec:=ArcCos(1/X);
  1341. end;
  1342. {$ENDIF}
  1343. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1344. function ArcSec(X: Double): Double;
  1345. begin
  1346. ArcSec:=ArcCos(1/X);
  1347. end;
  1348. {$ENDIF}
  1349. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1350. function ArcSec(X: Extended): Extended;
  1351. begin
  1352. ArcSec:=ArcCos(1/X);
  1353. end;
  1354. {$ENDIF}
  1355. {$ifdef FPC_HAS_TYPE_SINGLE}
  1356. function ArcCsc(X: Single): Single;
  1357. begin
  1358. ArcCsc:=ArcSin(1/X);
  1359. end;
  1360. {$ENDIF}
  1361. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1362. function ArcCsc(X: Double): Double;
  1363. begin
  1364. ArcCsc:=ArcSin(1/X);
  1365. end;
  1366. {$ENDIF}
  1367. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1368. function ArcCsc(X: Extended): Extended;
  1369. begin
  1370. ArcCsc:=ArcSin(1/X);
  1371. end;
  1372. {$ENDIF}
  1373. {$ifdef FPC_HAS_TYPE_SINGLE}
  1374. function ArcCot(X: Single): Single;
  1375. begin
  1376. if x=0 then
  1377. ArcCot:=0.5*pi
  1378. else
  1379. ArcCot:=ArcTan(1/X);
  1380. end;
  1381. {$ENDIF}
  1382. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1383. function ArcCot(X: Double): Double;
  1384. begin
  1385. begin
  1386. if x=0 then
  1387. ArcCot:=0.5*pi
  1388. else
  1389. ArcCot:=ArcTan(1/X);
  1390. end;
  1391. end;
  1392. {$ENDIF}
  1393. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1394. function ArcCot(X: Extended): Extended;
  1395. begin
  1396. begin
  1397. if x=0 then
  1398. ArcCot:=0.5*pi
  1399. else
  1400. ArcCot:=ArcTan(1/X);
  1401. end;
  1402. end;
  1403. {$ENDIF}
  1404. {$ifdef FPC_HAS_TYPE_SINGLE}
  1405. function ArcSecH(X : Single): Single;
  1406. begin
  1407. ArcSecH:=ln((1+(sqrt(1.0-sqr(X))))/X); //replacing division inside ln() by subtracting 2 ln()'s seems to be slower
  1408. end;
  1409. {$ENDIF}
  1410. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1411. function ArcSecH(X : Double): Double;
  1412. begin
  1413. ArcSecH:=ln((1+(sqrt(1.0-sqr(X))))/X);
  1414. end;
  1415. {$ENDIF}
  1416. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1417. function ArcSecH(X : Extended): Extended;
  1418. begin
  1419. ArcSecH:=ln((1+(sqrt(1.0-sqr(X))))/X);
  1420. end;
  1421. {$ENDIF}
  1422. {$ifdef FPC_HAS_TYPE_SINGLE}
  1423. function ArcCscH(X: Single): Single;
  1424. begin
  1425. ArcCscH:=ln((1.0/X)+sqrt(1.0/(sqr(x))+1.0));
  1426. end;
  1427. {$ENDIF}
  1428. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1429. function ArcCscH(X: Double): Double;
  1430. begin
  1431. ArcCscH:=ln((1.0/X)+sqrt(1.0/(sqr(x))+1.0));
  1432. end;
  1433. {$ENDIF}
  1434. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1435. function ArcCscH(X: Extended): Extended;
  1436. begin
  1437. ArcCscH:=ln((1.0/X)+sqrt(1.0/(sqr(x))+1.0));
  1438. end;
  1439. {$ENDIF}
  1440. {$ifdef FPC_HAS_TYPE_SINGLE}
  1441. function ArcCotH(X: Single): Single;
  1442. begin
  1443. ArcCotH:=0.5*ln((x + 1.0)/(x - 1.0));
  1444. end;
  1445. {$ENDIF}
  1446. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1447. function ArcCotH(X: Double): Double;
  1448. begin
  1449. ArcCotH:=0.5*ln((x + 1.0)/(x - 1.0));
  1450. end;
  1451. {$ENDIF}
  1452. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1453. function ArcCotH(X: Extended): Extended;
  1454. begin
  1455. ArcCotH:=0.5*ln((x + 1.0)/(x - 1.0));
  1456. end;
  1457. {$ENDIF}
  1458. { hypot function from AMath library (C) Copyright 2009-2013 Wolfgang Ehrhardt }
  1459. function hypot(x,y : float) : float;
  1460. begin
  1461. x:=abs(x);
  1462. y:=abs(y);
  1463. if (x>y) then
  1464. hypot:=x*sqrt(1.0+sqr(y/x))
  1465. else if (x>0.0) then
  1466. hypot:=y*sqrt(1.0+sqr(x/y))
  1467. else
  1468. hypot:=y;
  1469. end;
  1470. function log10(x : float) : float;
  1471. begin
  1472. log10:=ln(x)*0.43429448190325182765; { 1/ln(10) }
  1473. end;
  1474. {$ifndef FPC_MATH_HAS_LOG2}
  1475. function log2(x : float) : float;
  1476. begin
  1477. log2:=ln(x)*1.4426950408889634079; { 1/ln(2) }
  1478. end;
  1479. {$endif FPC_MATH_HAS_LOG2}
  1480. function logn(n,x : float) : float;
  1481. begin
  1482. logn:=ln(x)/ln(n);
  1483. end;
  1484. { lnxp1 function from AMath library (C) Copyright 2009-2013 Wolfgang Ehrhardt }
  1485. function lnxp1(x : float) : float;
  1486. var
  1487. y: float;
  1488. begin
  1489. if (x>=4.0) then
  1490. lnxp1:=ln(1.0+x)
  1491. else
  1492. begin
  1493. y:=1.0+x;
  1494. if (y=1.0) then
  1495. lnxp1:=x
  1496. else
  1497. begin
  1498. lnxp1:=ln(y); { lnxp1(-1) = ln(0) = -Inf }
  1499. if y>0.0 then
  1500. lnxp1:=lnxp1+(x-(y-1.0))/y;
  1501. end;
  1502. end;
  1503. end;
  1504. function power(base,exponent : float) : float;
  1505. begin
  1506. if Exponent=0.0 then
  1507. result:=1.0
  1508. else if (base=0.0) and (exponent>0.0) then
  1509. result:=0.0
  1510. else if (frac(exponent)=0.0) and (abs(exponent)<=maxint) then
  1511. result:=intpower(base,trunc(exponent))
  1512. else
  1513. result:=exp(exponent * ln (base));
  1514. end;
  1515. function intpower(base : float;exponent : longint) : float;
  1516. begin
  1517. if exponent<0 then
  1518. begin
  1519. base:=1.0/base;
  1520. exponent:=-exponent;
  1521. end;
  1522. intpower:=1.0;
  1523. while exponent<>0 do
  1524. begin
  1525. if exponent and 1<>0 then
  1526. intpower:=intpower*base;
  1527. exponent:=exponent shr 1;
  1528. base:=sqr(base);
  1529. end;
  1530. end;
  1531. operator ** (base,exponent : float) e: float; inline;
  1532. begin
  1533. e:=power(base,exponent);
  1534. end;
  1535. operator ** (base,exponent : int64) res: int64;
  1536. begin
  1537. if exponent<0 then
  1538. begin
  1539. if base<=0 then
  1540. raise EInvalidArgument.Create('Non-positive base with negative exponent in **');
  1541. if base=1 then
  1542. res:=1
  1543. else
  1544. res:=0;
  1545. exit;
  1546. end;
  1547. res:=1;
  1548. while exponent<>0 do
  1549. begin
  1550. if exponent and 1<>0 then
  1551. res:=res*base;
  1552. exponent:=exponent shr 1;
  1553. base:=base*base;
  1554. end;
  1555. end;
  1556. function ceil(x : float) : integer;
  1557. begin
  1558. Result:=Trunc(x)+ord(Frac(x)>0);
  1559. end;
  1560. function ceil64(x: float): Int64;
  1561. begin
  1562. Result:=Trunc(x)+ord(Frac(x)>0);
  1563. end;
  1564. function floor(x : float) : integer;
  1565. begin
  1566. Result:=Trunc(x)-ord(Frac(x)<0);
  1567. end;
  1568. function floor64(x: float): Int64;
  1569. begin
  1570. Result:=Trunc(x)-ord(Frac(x)<0);
  1571. end;
  1572. // Correction for "rounding to nearest, ties to even".
  1573. // RoundToNearestTieToEven(QWE.RTYUIOP) = QWE + TieToEven(ER, TYUIOP <> 0).
  1574. function TieToEven(AB: cardinal; somethingAfter: boolean): cardinal;
  1575. begin
  1576. result := AB and 1;
  1577. if (result <> 0) and not somethingAfter then
  1578. result := AB shr 1;
  1579. end;
  1580. {$ifdef FPC_HAS_TYPE_SINGLE}
  1581. procedure Frexp(X: single; out Mantissa: single; out Exponent: integer);
  1582. var
  1583. M: uint32;
  1584. E, ExtraE: int32;
  1585. begin
  1586. Mantissa := X;
  1587. E := TSingleRec(X).Exp;
  1588. if (E > 0) and (E < 2 * TSingleRec.Bias + 1) then
  1589. begin
  1590. // Normal.
  1591. TSingleRec(Mantissa).Exp := TSingleRec.Bias - 1;
  1592. Exponent := E - (TSingleRec.Bias - 1);
  1593. exit;
  1594. end;
  1595. if E = 0 then
  1596. begin
  1597. M := TSingleRec(X).Frac;
  1598. if M <> 0 then
  1599. begin
  1600. // Subnormal.
  1601. ExtraE := 23 - BsrDWord(M);
  1602. TSingleRec(Mantissa).Frac := M shl ExtraE; // "and (1 shl 23 - 1)" required to remove starting 1, but .SetFrac already does it.
  1603. TSingleRec(Mantissa).Exp := TSingleRec.Bias - 1;
  1604. Exponent := -TSingleRec.Bias + 2 - ExtraE;
  1605. exit;
  1606. end;
  1607. end;
  1608. // ±0, ±Inf, NaN.
  1609. Exponent := 0;
  1610. end;
  1611. function Ldexp(X: single; p: integer): single;
  1612. var
  1613. M, E: uint32;
  1614. xp, sh: integer;
  1615. begin
  1616. E := TSingleRec(X).Exp;
  1617. if (E = 0) and (TSingleRec(X).Frac = 0) or (E = 2 * TSingleRec.Bias + 1) then
  1618. // ±0, ±Inf, NaN.
  1619. exit(X);
  1620. Frexp(X, result, xp);
  1621. inc(xp, p);
  1622. if (xp >= -TSingleRec.Bias + 2) and (xp <= TSingleRec.Bias + 1) then
  1623. // Normalized.
  1624. TSingleRec(result).Exp := xp + (TSingleRec.Bias - 1)
  1625. else if xp > TSingleRec.Bias + 1 then
  1626. begin
  1627. // Overflow.
  1628. TSingleRec(result).Exp := 2 * TSingleRec.Bias + 1;
  1629. TSingleRec(result).Frac := 0;
  1630. end else
  1631. begin
  1632. TSingleRec(result).Exp := 0;
  1633. if xp >= -TSingleRec.Bias + 2 - 23 then
  1634. begin
  1635. // Denormalized.
  1636. M := TSingleRec(result).Frac or uint32(1) shl 23;
  1637. sh := -TSingleRec.Bias + 1 - xp;
  1638. TSingleRec(result).Frac := M shr (sh + 1) + TieToEven(M shr sh and 3, M and (uint32(1) shl sh - 1) <> 0);
  1639. end else
  1640. // Underflow.
  1641. TSingleRec(result).Frac := 0;
  1642. end;
  1643. end;
  1644. {$endif}
  1645. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1646. procedure Frexp(X: double; out Mantissa: double; out Exponent: integer);
  1647. var
  1648. M: uint64;
  1649. E, ExtraE: int32;
  1650. begin
  1651. Mantissa := X;
  1652. E := TDoubleRec(X).Exp;
  1653. if (E > 0) and (E < 2 * TDoubleRec.Bias + 1) then
  1654. begin
  1655. // Normal.
  1656. TDoubleRec(Mantissa).Exp := TDoubleRec.Bias - 1;
  1657. Exponent := E - (TDoubleRec.Bias - 1);
  1658. exit;
  1659. end;
  1660. if E = 0 then
  1661. begin
  1662. M := TDoubleRec(X).Frac;
  1663. if M <> 0 then
  1664. begin
  1665. // Subnormal.
  1666. ExtraE := 52 - BsrQWord(M);
  1667. TDoubleRec(Mantissa).Frac := M shl ExtraE; // "and (1 shl 52 - 1)" required to remove starting 1, but .SetFrac already does it.
  1668. TDoubleRec(Mantissa).Exp := TDoubleRec.Bias - 1;
  1669. Exponent := -TDoubleRec.Bias + 2 - ExtraE;
  1670. exit;
  1671. end;
  1672. end;
  1673. // ±0, ±Inf, NaN.
  1674. Exponent := 0;
  1675. end;
  1676. function Ldexp(X: double; p: integer): double;
  1677. var
  1678. M: uint64;
  1679. E: uint32;
  1680. xp, sh: integer;
  1681. begin
  1682. E := TDoubleRec(X).Exp;
  1683. if (E = 0) and (TDoubleRec(X).Frac = 0) or (E = 2 * TDoubleRec.Bias + 1) then
  1684. // ±0, ±Inf, NaN.
  1685. exit(X);
  1686. Frexp(X, result, xp);
  1687. inc(xp, p);
  1688. if (xp >= -TDoubleRec.Bias + 2) and (xp <= TDoubleRec.Bias + 1) then
  1689. // Normalized.
  1690. TDoubleRec(result).Exp := xp + (TDoubleRec.Bias - 1)
  1691. else if xp > TDoubleRec.Bias + 1 then
  1692. begin
  1693. // Overflow.
  1694. TDoubleRec(result).Exp := 2 * TDoubleRec.Bias + 1;
  1695. TDoubleRec(result).Frac := 0;
  1696. end else
  1697. begin
  1698. TDoubleRec(result).Exp := 0;
  1699. if xp >= -TDoubleRec.Bias + 2 - 52 then
  1700. begin
  1701. // Denormalized.
  1702. M := TDoubleRec(result).Frac or uint64(1) shl 52;
  1703. sh := -TSingleRec.Bias + 1 - xp;
  1704. TDoubleRec(result).Frac := M shr (sh + 1) + TieToEven(M shr sh and 3, M and (uint64(1) shl sh - 1) <> 0);
  1705. end else
  1706. // Underflow.
  1707. TDoubleRec(result).Frac := 0;
  1708. end;
  1709. end;
  1710. {$endif}
  1711. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1712. procedure Frexp(X: extended; out Mantissa: extended; out Exponent: integer);
  1713. var
  1714. M: uint64;
  1715. E, ExtraE: int32;
  1716. begin
  1717. Mantissa := X;
  1718. E := TExtended80Rec(X).Exp;
  1719. if (E > 0) and (E < 2 * TExtended80Rec.Bias + 1) then
  1720. begin
  1721. // Normal.
  1722. TExtended80Rec(Mantissa).Exp := TExtended80Rec.Bias - 1;
  1723. Exponent := E - (TExtended80Rec.Bias - 1);
  1724. exit;
  1725. end;
  1726. if E = 0 then
  1727. begin
  1728. M := TExtended80Rec(X).Frac;
  1729. if M <> 0 then
  1730. begin
  1731. // Subnormal. Extended has explicit starting 1.
  1732. ExtraE := 63 - BsrQWord(M);
  1733. TExtended80Rec(Mantissa).Frac := M shl ExtraE;
  1734. TExtended80Rec(Mantissa).Exp := TExtended80Rec.Bias - 1;
  1735. Exponent := -TExtended80Rec.Bias + 2 - ExtraE;
  1736. exit;
  1737. end;
  1738. end;
  1739. // ±0, ±Inf, NaN.
  1740. Exponent := 0;
  1741. end;
  1742. function Ldexp(X: extended; p: integer): extended;
  1743. var
  1744. M: uint64;
  1745. E: uint32;
  1746. xp, sh: integer;
  1747. begin
  1748. E := TExtended80Rec(X).Exp;
  1749. if (E = 0) and (TExtended80Rec(X).Frac = 0) or (E = 2 * TExtended80Rec.Bias + 1) then
  1750. // ±0, ±Inf, NaN.
  1751. exit(X);
  1752. Frexp(X, result, xp);
  1753. inc(xp, p);
  1754. if (xp >= -TExtended80Rec.Bias + 2) and (xp <= TExtended80Rec.Bias + 1) then
  1755. // Normalized.
  1756. TExtended80Rec(result).Exp := xp + (TExtended80Rec.Bias - 1)
  1757. else if xp > TExtended80Rec.Bias + 1 then
  1758. begin
  1759. // Overflow.
  1760. TExtended80Rec(result).Exp := 2 * TExtended80Rec.Bias + 1;
  1761. TExtended80Rec(result).Frac := uint64(1) shl 63;
  1762. end
  1763. else if xp >= -TExtended80Rec.Bias + 2 - 63 then
  1764. begin
  1765. // Denormalized... usually.
  1766. // Mantissa of subnormal 'extended' (Exp = 0) must always start with 0.
  1767. // If the calculated mantissa starts with 1, extended instead becomes normalized with Exp = 1.
  1768. M := TExtended80Rec(result).Frac;
  1769. sh := -TExtended80Rec.Bias + 1 - xp;
  1770. M := M shr (sh + 1) + TieToEven(M shr sh and 3, M and (uint64(1) shl sh - 1) <> 0);
  1771. TExtended80Rec(result).Exp := M shr 63;
  1772. TExtended80Rec(result).Frac := M;
  1773. end else
  1774. begin
  1775. // Underflow.
  1776. TExtended80Rec(result).Exp := 0;
  1777. TExtended80Rec(result).Frac := 0;
  1778. end;
  1779. end;
  1780. {$endif}
  1781. const
  1782. { Cutoff for https://en.wikipedia.org/wiki/Pairwise_summation; sums of at least this many elements are split in two halves. }
  1783. RecursiveSumThreshold=12;
  1784. {$ifdef FPC_HAS_TYPE_SINGLE}
  1785. function mean(const data : array of Single) : float;
  1786. begin
  1787. Result:=Mean(PSingle(@data[0]),High(Data)+1);
  1788. end;
  1789. function mean(const data : PSingle; Const N : longint) : float;
  1790. begin
  1791. mean:=sum(Data,N);
  1792. mean:=mean/N;
  1793. end;
  1794. function sum(const data : array of Single) : float;inline;
  1795. begin
  1796. Result:=Sum(PSingle(@Data[0]),High(Data)+1);
  1797. end;
  1798. function sum(const data : PSingle;Const N : longint) : float;
  1799. var
  1800. i : SizeInt;
  1801. begin
  1802. if N>=RecursiveSumThreshold then
  1803. result:=sum(data,longword(N) div 2)+sum(data+longword(N) div 2,N-longword(N) div 2)
  1804. else
  1805. begin
  1806. result:=0;
  1807. for i:=0 to N-1 do
  1808. result:=result+data[i];
  1809. end;
  1810. end;
  1811. {$endif FPC_HAS_TYPE_SINGLE}
  1812. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1813. function mean(const data : array of Double) : float; inline;
  1814. begin
  1815. Result:=Mean(PDouble(@data[0]),High(Data)+1);
  1816. end;
  1817. function mean(const data : PDouble; Const N : longint) : float;
  1818. begin
  1819. mean:=sum(Data,N);
  1820. mean:=mean/N;
  1821. end;
  1822. function sum(const data : array of Double) : float; inline;
  1823. begin
  1824. Result:=Sum(PDouble(@Data[0]),High(Data)+1);
  1825. end;
  1826. function sum(const data : PDouble;Const N : longint) : float;
  1827. var
  1828. i : SizeInt;
  1829. begin
  1830. if N>=RecursiveSumThreshold then
  1831. result:=sum(data,longword(N) div 2)+sum(data+longword(N) div 2,N-longword(N) div 2)
  1832. else
  1833. begin
  1834. result:=0;
  1835. for i:=0 to N-1 do
  1836. result:=result+data[i];
  1837. end;
  1838. end;
  1839. {$endif FPC_HAS_TYPE_DOUBLE}
  1840. {$ifdef FPC_HAS_TYPE_EXTENDED}
  1841. function mean(const data : array of Extended) : float;
  1842. begin
  1843. Result:=Mean(PExtended(@data[0]),High(Data)+1);
  1844. end;
  1845. function mean(const data : PExtended; Const N : longint) : float;
  1846. begin
  1847. mean:=sum(Data,N);
  1848. mean:=mean/N;
  1849. end;
  1850. function sum(const data : array of Extended) : float; inline;
  1851. begin
  1852. Result:=Sum(PExtended(@Data[0]),High(Data)+1);
  1853. end;
  1854. function sum(const data : PExtended;Const N : longint) : float;
  1855. var
  1856. i : SizeInt;
  1857. begin
  1858. if N>=RecursiveSumThreshold then
  1859. result:=sum(data,longword(N) div 2)+sum(data+longword(N) div 2,N-longword(N) div 2)
  1860. else
  1861. begin
  1862. result:=0;
  1863. for i:=0 to N-1 do
  1864. result:=result+data[i];
  1865. end;
  1866. end;
  1867. {$endif FPC_HAS_TYPE_EXTENDED}
  1868. function sumInt(const data : PInt64;Const N : longint) : Int64;
  1869. var
  1870. i : SizeInt;
  1871. begin
  1872. sumInt:=0;
  1873. for i:=0 to N-1 do
  1874. sumInt:=sumInt+data[i];
  1875. end;
  1876. function sumInt(const data : array of Int64) : Int64; inline;
  1877. begin
  1878. Result:=SumInt(PInt64(@Data[0]),High(Data)+1);
  1879. end;
  1880. function mean(const data : PInt64; const N : Longint):Float;
  1881. begin
  1882. mean:=sumInt(Data,N);
  1883. mean:=mean/N;
  1884. end;
  1885. function mean(const data: array of Int64):Float;
  1886. begin
  1887. mean:=mean(PInt64(@data[0]),High(Data)+1);
  1888. end;
  1889. function sumInt(const data : PInteger; Const N : longint) : Int64;
  1890. var
  1891. i : SizeInt;
  1892. begin
  1893. sumInt:=0;
  1894. for i:=0 to N-1 do
  1895. sumInt:=sumInt+data[i];
  1896. end;
  1897. function sumInt(const data : array of Integer) : Int64;inline;
  1898. begin
  1899. Result:=sumInt(PInteger(@Data[0]),High(Data)+1);
  1900. end;
  1901. function mean(const data : PInteger; const N : Longint):Float;
  1902. begin
  1903. mean:=sumInt(Data,N);
  1904. mean:=mean/N;
  1905. end;
  1906. function mean(const data: array of Integer):Float;
  1907. begin
  1908. mean:=mean(PInteger(@data[0]),High(Data)+1);
  1909. end;
  1910. {$ifdef FPC_HAS_TYPE_SINGLE}
  1911. function sumofsquares(const data : array of Single) : float; inline;
  1912. begin
  1913. Result:=sumofsquares(PSingle(@data[0]),High(Data)+1);
  1914. end;
  1915. function sumofsquares(const data : PSingle; Const N : Integer) : float;
  1916. var
  1917. i : SizeInt;
  1918. begin
  1919. if N>=RecursiveSumThreshold then
  1920. result:=sumofsquares(data,cardinal(N) div 2)+sumofsquares(data+cardinal(N) div 2,N-cardinal(N) div 2)
  1921. else
  1922. begin
  1923. result:=0;
  1924. for i:=0 to N-1 do
  1925. result:=result+sqr(data[i]);
  1926. end;
  1927. end;
  1928. procedure sumsandsquares(const data : array of Single;
  1929. var sum,sumofsquares : float); inline;
  1930. begin
  1931. sumsandsquares (PSingle(@Data[0]),High(Data)+1,Sum,sumofsquares);
  1932. end;
  1933. procedure sumsandsquares(const data : PSingle; Const N : Integer;
  1934. var sum,sumofsquares : float);
  1935. var
  1936. i : SizeInt;
  1937. temp,tsum,tsumofsquares,sum0,sumofsquares0,sum1,sumofsquares1 : float;
  1938. begin
  1939. if N>=RecursiveSumThreshold then
  1940. begin
  1941. sumsandsquares(data,cardinal(N) div 2,sum0,sumofsquares0);
  1942. sumsandsquares(data+cardinal(N) div 2,N-cardinal(N) div 2,sum1,sumofsquares1);
  1943. sum:=sum0+sum1;
  1944. sumofsquares:=sumofsquares0+sumofsquares1;
  1945. end
  1946. else
  1947. begin
  1948. tsum:=0;
  1949. tsumofsquares:=0;
  1950. for i:=0 to N-1 do
  1951. begin
  1952. temp:=data[i];
  1953. tsum:=tsum+temp;
  1954. tsumofsquares:=tsumofsquares+sqr(temp);
  1955. end;
  1956. sum:=tsum;
  1957. sumofsquares:=tsumofsquares;
  1958. end;
  1959. end;
  1960. {$endif FPC_HAS_TYPE_SINGLE}
  1961. {$ifdef FPC_HAS_TYPE_DOUBLE}
  1962. function sumofsquares(const data : array of Double) : float; inline;
  1963. begin
  1964. Result:=sumofsquares(PDouble(@data[0]),High(Data)+1);
  1965. end;
  1966. function sumofsquares(const data : PDouble; Const N : Integer) : float;
  1967. var
  1968. i : SizeInt;
  1969. begin
  1970. if N>=RecursiveSumThreshold then
  1971. result:=sumofsquares(data,cardinal(N) div 2)+sumofsquares(data+cardinal(N) div 2,N-cardinal(N) div 2)
  1972. else
  1973. begin
  1974. result:=0;
  1975. for i:=0 to N-1 do
  1976. result:=result+sqr(data[i]);
  1977. end;
  1978. end;
  1979. procedure sumsandsquares(const data : array of Double;
  1980. var sum,sumofsquares : float); inline;
  1981. begin
  1982. sumsandsquares (PDouble(@Data[0]),High(Data)+1,Sum,sumofsquares);
  1983. end;
  1984. procedure sumsandsquares(const data : PDouble; Const N : Integer;
  1985. var sum,sumofsquares : float);
  1986. var
  1987. i : SizeInt;
  1988. temp,tsum,tsumofsquares,sum0,sumofsquares0,sum1,sumofsquares1 : float;
  1989. begin
  1990. if N>=RecursiveSumThreshold then
  1991. begin
  1992. sumsandsquares(data,cardinal(N) div 2,sum0,sumofsquares0);
  1993. sumsandsquares(data+cardinal(N) div 2,N-cardinal(N) div 2,sum1,sumofsquares1);
  1994. sum:=sum0+sum1;
  1995. sumofsquares:=sumofsquares0+sumofsquares1;
  1996. end
  1997. else
  1998. begin
  1999. tsum:=0;
  2000. tsumofsquares:=0;
  2001. for i:=0 to N-1 do
  2002. begin
  2003. temp:=data[i];
  2004. tsum:=tsum+temp;
  2005. tsumofsquares:=tsumofsquares+sqr(temp);
  2006. end;
  2007. sum:=tsum;
  2008. sumofsquares:=tsumofsquares;
  2009. end;
  2010. end;
  2011. {$endif FPC_HAS_TYPE_DOUBLE}
  2012. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2013. function sumofsquares(const data : array of Extended) : float; inline;
  2014. begin
  2015. Result:=sumofsquares(PExtended(@data[0]),High(Data)+1);
  2016. end;
  2017. function sumofsquares(const data : PExtended; Const N : Integer) : float;
  2018. var
  2019. i : SizeInt;
  2020. begin
  2021. if N>=RecursiveSumThreshold then
  2022. result:=sumofsquares(data,cardinal(N) div 2)+sumofsquares(data+cardinal(N) div 2,N-cardinal(N) div 2)
  2023. else
  2024. begin
  2025. result:=0;
  2026. for i:=0 to N-1 do
  2027. result:=result+sqr(data[i]);
  2028. end;
  2029. end;
  2030. procedure sumsandsquares(const data : array of Extended;
  2031. var sum,sumofsquares : float); inline;
  2032. begin
  2033. sumsandsquares (PExtended(@Data[0]),High(Data)+1,Sum,sumofsquares);
  2034. end;
  2035. procedure sumsandsquares(const data : PExtended; Const N : Integer;
  2036. var sum,sumofsquares : float);
  2037. var
  2038. i : SizeInt;
  2039. temp,tsum,tsumofsquares,sum0,sumofsquares0,sum1,sumofsquares1 : float;
  2040. begin
  2041. if N>=RecursiveSumThreshold then
  2042. begin
  2043. sumsandsquares(data,cardinal(N) div 2,sum0,sumofsquares0);
  2044. sumsandsquares(data+cardinal(N) div 2,N-cardinal(N) div 2,sum1,sumofsquares1);
  2045. sum:=sum0+sum1;
  2046. sumofsquares:=sumofsquares0+sumofsquares1;
  2047. end
  2048. else
  2049. begin
  2050. tsum:=0;
  2051. tsumofsquares:=0;
  2052. for i:=0 to N-1 do
  2053. begin
  2054. temp:=data[i];
  2055. tsum:=tsum+temp;
  2056. tsumofsquares:=tsumofsquares+sqr(temp);
  2057. end;
  2058. sum:=tsum;
  2059. sumofsquares:=tsumofsquares;
  2060. end;
  2061. end;
  2062. {$endif FPC_HAS_TYPE_EXTENDED}
  2063. function randg(mean,stddev : float) : float;
  2064. Var U1,S2 : Float;
  2065. begin
  2066. repeat
  2067. u1:= 2*random-1;
  2068. S2:=Sqr(U1)+sqr(2*random-1);
  2069. until s2<1;
  2070. randg:=Sqrt(-2*ln(S2)/S2)*u1*stddev+Mean;
  2071. end;
  2072. function RandomRange(const aFrom, aTo: Integer): Integer;
  2073. begin
  2074. Result:=Random(Abs(aFrom-aTo))+Min(aTo,AFrom);
  2075. end;
  2076. function RandomRange(const aFrom, aTo: Int64): Int64;
  2077. begin
  2078. Result:=Random(Abs(aFrom-aTo))+Min(aTo,AFrom);
  2079. end;
  2080. {$ifdef FPC_HAS_TYPE_SINGLE}
  2081. procedure MeanAndTotalVariance
  2082. (const data: PSingle; N: LongInt; var mu, variance: float);
  2083. function CalcVariance(data: PSingle; N: SizeInt; mu: float): float;
  2084. var
  2085. i: SizeInt;
  2086. begin
  2087. if N>=RecursiveSumThreshold then
  2088. result:=CalcVariance(data,SizeUint(N) div 2,mu)+CalcVariance(data+SizeUint(N) div 2,N-SizeUint(N) div 2,mu)
  2089. else
  2090. begin
  2091. result:=0;
  2092. for i:=0 to N-1 do
  2093. result:=result+Sqr(data[i]-mu);
  2094. end;
  2095. end;
  2096. begin
  2097. mu := Mean( data, N );
  2098. variance := CalcVariance( data, N, mu );
  2099. end;
  2100. function stddev(const data : array of Single) : float; inline;
  2101. begin
  2102. Result:=Stddev(PSingle(@Data[0]),High(Data)+1);
  2103. end;
  2104. function stddev(const data : PSingle; Const N : Integer) : float;
  2105. begin
  2106. StdDev:=Sqrt(Variance(Data,N));
  2107. end;
  2108. procedure meanandstddev(const data : array of Single;
  2109. var mean,stddev : float); inline;
  2110. begin
  2111. Meanandstddev(PSingle(@Data[0]),High(Data)+1,Mean,stddev);
  2112. end;
  2113. procedure meanandstddev
  2114. ( const data: PSingle;
  2115. const N: Longint;
  2116. var mean,
  2117. stdDev: Float
  2118. );
  2119. var totalVariance: float;
  2120. begin
  2121. MeanAndTotalVariance( data, N, mean, totalVariance );
  2122. if N < 2 then stdDev := 0
  2123. else stdDev := Sqrt( totalVariance / ( N - 1 ) );
  2124. end;
  2125. function variance(const data : array of Single) : float; inline;
  2126. begin
  2127. Variance:=Variance(PSingle(@Data[0]),High(Data)+1);
  2128. end;
  2129. function variance(const data : PSingle; Const N : Integer) : float;
  2130. begin
  2131. If N=1 then
  2132. Result:=0
  2133. else
  2134. Result:=TotalVariance(Data,N)/(N-1);
  2135. end;
  2136. function totalvariance(const data : array of Single) : float; inline;
  2137. begin
  2138. Result:=TotalVariance(PSingle(@Data[0]),High(Data)+1);
  2139. end;
  2140. function totalvariance(const data : PSingle; const N : Integer) : float;
  2141. var mu: float;
  2142. begin
  2143. MeanAndTotalVariance( data, N, mu, result );
  2144. end;
  2145. function popnstddev(const data : array of Single) : float;
  2146. begin
  2147. PopnStdDev:=Sqrt(PopnVariance(PSingle(@Data[0]),High(Data)+1));
  2148. end;
  2149. function popnstddev(const data : PSingle; Const N : Integer) : float;
  2150. begin
  2151. PopnStdDev:=Sqrt(PopnVariance(Data,N));
  2152. end;
  2153. function popnvariance(const data : array of Single) : float; inline;
  2154. begin
  2155. popnvariance:=popnvariance(PSingle(@data[0]),high(Data)+1);
  2156. end;
  2157. function popnvariance(const data : PSingle; Const N : Integer) : float;
  2158. begin
  2159. PopnVariance:=TotalVariance(Data,N)/N;
  2160. end;
  2161. procedure momentskewkurtosis(const data : array of single;
  2162. out m1,m2,m3,m4,skew,kurtosis : float); inline;
  2163. begin
  2164. momentskewkurtosis(PSingle(@Data[0]),High(Data)+1,m1,m2,m3,m4,skew,kurtosis);
  2165. end;
  2166. type
  2167. TMoments2to4 = array[2 .. 4] of float;
  2168. procedure momentskewkurtosis(
  2169. const data: pSingle;
  2170. Const N: integer;
  2171. out m1: float;
  2172. out m2: float;
  2173. out m3: float;
  2174. out m4: float;
  2175. out skew: float;
  2176. out kurtosis: float
  2177. );
  2178. procedure CalcDevSums2to4(data: PSingle; N: SizeInt; m1: float; out m2to4: TMoments2to4);
  2179. var
  2180. tm2, tm3, tm4, dev, dev2: float;
  2181. i: SizeInt;
  2182. m2to4Part0, m2to4Part1: TMoments2to4;
  2183. begin
  2184. if N >= RecursiveSumThreshold then
  2185. begin
  2186. CalcDevSums2to4(data, SizeUint(N) div 2, m1, m2to4Part0);
  2187. CalcDevSums2to4(data + SizeUint(N) div 2, N - SizeUint(N) div 2, m1, m2to4Part1);
  2188. for i := Low(TMoments2to4) to High(TMoments2to4) do
  2189. m2to4[i] := m2to4Part0[i] + m2to4Part1[i];
  2190. end
  2191. else
  2192. begin
  2193. tm2 := 0;
  2194. tm3 := 0;
  2195. tm4 := 0;
  2196. for i := 0 to N - 1 do
  2197. begin
  2198. dev := data[i] - m1;
  2199. dev2 := sqr(dev);
  2200. tm2 := tm2 + dev2;
  2201. tm3 := tm3 + dev2 * dev;
  2202. tm4 := tm4 + sqr(dev2);
  2203. end;
  2204. m2to4[2] := tm2;
  2205. m2to4[3] := tm3;
  2206. m2to4[4] := tm4;
  2207. end;
  2208. end;
  2209. var
  2210. reciprocalN: float;
  2211. m2to4: TMoments2to4;
  2212. begin
  2213. m1 := 0;
  2214. reciprocalN := 1/N;
  2215. m1 := reciprocalN * sum(data, N);
  2216. CalcDevSums2to4(data, N, m1, m2to4);
  2217. m2 := reciprocalN * m2to4[2];
  2218. m3 := reciprocalN * m2to4[3];
  2219. m4 := reciprocalN * m2to4[4];
  2220. skew := m3 / (sqrt(m2)*m2);
  2221. kurtosis := m4 / (m2 * m2);
  2222. end;
  2223. function norm(const data : array of Single) : float; inline;
  2224. begin
  2225. norm:=Norm(PSingle(@data[0]),High(Data)+1);
  2226. end;
  2227. function norm(const data : PSingle; Const N : Integer) : float;
  2228. begin
  2229. norm:=sqrt(sumofsquares(data,N));
  2230. end;
  2231. {$endif FPC_HAS_TYPE_SINGLE}
  2232. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2233. procedure MeanAndTotalVariance
  2234. (const data: PDouble; N: LongInt; var mu, variance: float);
  2235. function CalcVariance(data: PDouble; N: SizeInt; mu: float): float;
  2236. var
  2237. i: SizeInt;
  2238. begin
  2239. if N>=RecursiveSumThreshold then
  2240. result:=CalcVariance(data,SizeUint(N) div 2,mu)+CalcVariance(data+SizeUint(N) div 2,N-SizeUint(N) div 2,mu)
  2241. else
  2242. begin
  2243. result:=0;
  2244. for i:=0 to N-1 do
  2245. result:=result+Sqr(data[i]-mu);
  2246. end;
  2247. end;
  2248. begin
  2249. mu := Mean( data, N );
  2250. variance := CalcVariance( data, N, mu );
  2251. end;
  2252. function stddev(const data : array of Double) : float; inline;
  2253. begin
  2254. Result:=Stddev(PDouble(@Data[0]),High(Data)+1)
  2255. end;
  2256. function stddev(const data : PDouble; Const N : Integer) : float;
  2257. begin
  2258. StdDev:=Sqrt(Variance(Data,N));
  2259. end;
  2260. procedure meanandstddev(const data : array of Double;
  2261. var mean,stddev : float);
  2262. begin
  2263. Meanandstddev(PDouble(@Data[0]),High(Data)+1,Mean,stddev);
  2264. end;
  2265. procedure meanandstddev
  2266. ( const data: PDouble;
  2267. const N: Longint;
  2268. var mean,
  2269. stdDev: Float
  2270. );
  2271. var totalVariance: float;
  2272. begin
  2273. MeanAndTotalVariance( data, N, mean, totalVariance );
  2274. if N < 2 then stdDev := 0
  2275. else stdDev := Sqrt( totalVariance / ( N - 1 ) );
  2276. end;
  2277. function variance(const data : array of Double) : float; inline;
  2278. begin
  2279. Variance:=Variance(PDouble(@Data[0]),High(Data)+1);
  2280. end;
  2281. function variance(const data : PDouble; Const N : Integer) : float;
  2282. begin
  2283. If N=1 then
  2284. Result:=0
  2285. else
  2286. Result:=TotalVariance(Data,N)/(N-1);
  2287. end;
  2288. function totalvariance(const data : array of Double) : float; inline;
  2289. begin
  2290. Result:=TotalVariance(PDouble(@Data[0]),High(Data)+1);
  2291. end;
  2292. function totalvariance(const data : PDouble; const N : Integer) : float;
  2293. var mu: float;
  2294. begin
  2295. MeanAndTotalVariance( data, N, mu, result );
  2296. end;
  2297. function popnstddev(const data : array of Double) : float;
  2298. begin
  2299. PopnStdDev:=Sqrt(PopnVariance(PDouble(@Data[0]),High(Data)+1));
  2300. end;
  2301. function popnstddev(const data : PDouble; Const N : Integer) : float;
  2302. begin
  2303. PopnStdDev:=Sqrt(PopnVariance(Data,N));
  2304. end;
  2305. function popnvariance(const data : array of Double) : float; inline;
  2306. begin
  2307. popnvariance:=popnvariance(PDouble(@data[0]),high(Data)+1);
  2308. end;
  2309. function popnvariance(const data : PDouble; Const N : Integer) : float;
  2310. begin
  2311. PopnVariance:=TotalVariance(Data,N)/N;
  2312. end;
  2313. procedure momentskewkurtosis(const data : array of Double;
  2314. out m1,m2,m3,m4,skew,kurtosis : float);
  2315. begin
  2316. momentskewkurtosis(PDouble(@Data[0]),High(Data)+1,m1,m2,m3,m4,skew,kurtosis);
  2317. end;
  2318. procedure momentskewkurtosis(
  2319. const data: pdouble;
  2320. Const N: integer;
  2321. out m1: float;
  2322. out m2: float;
  2323. out m3: float;
  2324. out m4: float;
  2325. out skew: float;
  2326. out kurtosis: float
  2327. );
  2328. procedure CalcDevSums2to4(data: PDouble; N: SizeInt; m1: float; out m2to4: TMoments2to4);
  2329. var
  2330. tm2, tm3, tm4, dev, dev2: float;
  2331. i: SizeInt;
  2332. m2to4Part0, m2to4Part1: TMoments2to4;
  2333. begin
  2334. if N >= RecursiveSumThreshold then
  2335. begin
  2336. CalcDevSums2to4(data, SizeUint(N) div 2, m1, m2to4Part0);
  2337. CalcDevSums2to4(data + SizeUint(N) div 2, N - SizeUint(N) div 2, m1, m2to4Part1);
  2338. for i := Low(TMoments2to4) to High(TMoments2to4) do
  2339. m2to4[i] := m2to4Part0[i] + m2to4Part1[i];
  2340. end
  2341. else
  2342. begin
  2343. tm2 := 0;
  2344. tm3 := 0;
  2345. tm4 := 0;
  2346. for i := 0 to N - 1 do
  2347. begin
  2348. dev := data[i] - m1;
  2349. dev2 := sqr(dev);
  2350. tm2 := tm2 + dev2;
  2351. tm3 := tm3 + dev2 * dev;
  2352. tm4 := tm4 + sqr(dev2);
  2353. end;
  2354. m2to4[2] := tm2;
  2355. m2to4[3] := tm3;
  2356. m2to4[4] := tm4;
  2357. end;
  2358. end;
  2359. var
  2360. reciprocalN: float;
  2361. m2to4: TMoments2to4;
  2362. begin
  2363. m1 := 0;
  2364. reciprocalN := 1/N;
  2365. m1 := reciprocalN * sum(data, N);
  2366. CalcDevSums2to4(data, N, m1, m2to4);
  2367. m2 := reciprocalN * m2to4[2];
  2368. m3 := reciprocalN * m2to4[3];
  2369. m4 := reciprocalN * m2to4[4];
  2370. skew := m3 / (sqrt(m2)*m2);
  2371. kurtosis := m4 / (m2 * m2);
  2372. end;
  2373. function norm(const data : array of Double) : float; inline;
  2374. begin
  2375. norm:=Norm(PDouble(@data[0]),High(Data)+1);
  2376. end;
  2377. function norm(const data : PDouble; Const N : Integer) : float;
  2378. begin
  2379. norm:=sqrt(sumofsquares(data,N));
  2380. end;
  2381. {$endif FPC_HAS_TYPE_DOUBLE}
  2382. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2383. procedure MeanAndTotalVariance
  2384. (const data: PExtended; N: LongInt; var mu, variance: float);
  2385. function CalcVariance(data: PExtended; N: SizeInt; mu: float): float;
  2386. var
  2387. i: SizeInt;
  2388. begin
  2389. if N>=RecursiveSumThreshold then
  2390. result:=CalcVariance(data,SizeUint(N) div 2,mu)+CalcVariance(data+SizeUint(N) div 2,N-SizeUint(N) div 2,mu)
  2391. else
  2392. begin
  2393. result:=0;
  2394. for i:=0 to N-1 do
  2395. result:=result+Sqr(data[i]-mu);
  2396. end;
  2397. end;
  2398. begin
  2399. mu := Mean( data, N );
  2400. variance := CalcVariance( data, N, mu );
  2401. end;
  2402. function stddev(const data : array of Extended) : float; inline;
  2403. begin
  2404. Result:=Stddev(PExtended(@Data[0]),High(Data)+1)
  2405. end;
  2406. function stddev(const data : PExtended; Const N : Integer) : float;
  2407. begin
  2408. StdDev:=Sqrt(Variance(Data,N));
  2409. end;
  2410. procedure meanandstddev(const data : array of Extended;
  2411. var mean,stddev : float); inline;
  2412. begin
  2413. Meanandstddev(PExtended(@Data[0]),High(Data)+1,Mean,stddev);
  2414. end;
  2415. procedure meanandstddev
  2416. ( const data: PExtended;
  2417. const N: Longint;
  2418. var mean,
  2419. stdDev: Float
  2420. );
  2421. var totalVariance: float;
  2422. begin
  2423. MeanAndTotalVariance( data, N, mean, totalVariance );
  2424. if N < 2 then stdDev := 0
  2425. else stdDev := Sqrt( totalVariance / ( N - 1 ) );
  2426. end;
  2427. function variance(const data : array of Extended) : float; inline;
  2428. begin
  2429. Variance:=Variance(PExtended(@Data[0]),High(Data)+1);
  2430. end;
  2431. function variance(const data : PExtended; Const N : Integer) : float;
  2432. begin
  2433. If N=1 then
  2434. Result:=0
  2435. else
  2436. Result:=TotalVariance(Data,N)/(N-1);
  2437. end;
  2438. function totalvariance(const data : array of Extended) : float; inline;
  2439. begin
  2440. Result:=TotalVariance(PExtended(@Data[0]),High(Data)+1);
  2441. end;
  2442. function totalvariance(const data : PExtended;Const N : Integer) : float;
  2443. var mu: float;
  2444. begin
  2445. MeanAndTotalVariance( data, N, mu, result );
  2446. end;
  2447. function popnstddev(const data : array of Extended) : float;
  2448. begin
  2449. PopnStdDev:=Sqrt(PopnVariance(PExtended(@Data[0]),High(Data)+1));
  2450. end;
  2451. function popnstddev(const data : PExtended; Const N : Integer) : float;
  2452. begin
  2453. PopnStdDev:=Sqrt(PopnVariance(Data,N));
  2454. end;
  2455. function popnvariance(const data : array of Extended) : float; inline;
  2456. begin
  2457. popnvariance:=popnvariance(PExtended(@data[0]),high(Data)+1);
  2458. end;
  2459. function popnvariance(const data : PExtended; Const N : Integer) : float;
  2460. begin
  2461. PopnVariance:=TotalVariance(Data,N)/N;
  2462. end;
  2463. procedure momentskewkurtosis(const data : array of Extended;
  2464. out m1,m2,m3,m4,skew,kurtosis : float); inline;
  2465. begin
  2466. momentskewkurtosis(PExtended(@Data[0]),High(Data)+1,m1,m2,m3,m4,skew,kurtosis);
  2467. end;
  2468. procedure momentskewkurtosis(
  2469. const data: pExtended;
  2470. Const N: Integer;
  2471. out m1: float;
  2472. out m2: float;
  2473. out m3: float;
  2474. out m4: float;
  2475. out skew: float;
  2476. out kurtosis: float
  2477. );
  2478. procedure CalcDevSums2to4(data: PExtended; N: SizeInt; m1: float; out m2to4: TMoments2to4);
  2479. var
  2480. tm2, tm3, tm4, dev, dev2: float;
  2481. i: SizeInt;
  2482. m2to4Part0, m2to4Part1: TMoments2to4;
  2483. begin
  2484. if N >= RecursiveSumThreshold then
  2485. begin
  2486. CalcDevSums2to4(data, SizeUint(N) div 2, m1, m2to4Part0);
  2487. CalcDevSums2to4(data + SizeUint(N) div 2, N - SizeUint(N) div 2, m1, m2to4Part1);
  2488. for i := Low(TMoments2to4) to High(TMoments2to4) do
  2489. m2to4[i] := m2to4Part0[i] + m2to4Part1[i];
  2490. end
  2491. else
  2492. begin
  2493. tm2 := 0;
  2494. tm3 := 0;
  2495. tm4 := 0;
  2496. for i := 0 to N - 1 do
  2497. begin
  2498. dev := data[i] - m1;
  2499. dev2 := sqr(dev);
  2500. tm2 := tm2 + dev2;
  2501. tm3 := tm3 + dev2 * dev;
  2502. tm4 := tm4 + sqr(dev2);
  2503. end;
  2504. m2to4[2] := tm2;
  2505. m2to4[3] := tm3;
  2506. m2to4[4] := tm4;
  2507. end;
  2508. end;
  2509. var
  2510. reciprocalN: float;
  2511. m2to4: TMoments2to4;
  2512. begin
  2513. m1 := 0;
  2514. reciprocalN := 1/N;
  2515. m1 := reciprocalN * sum(data, N);
  2516. CalcDevSums2to4(data, N, m1, m2to4);
  2517. m2 := reciprocalN * m2to4[2];
  2518. m3 := reciprocalN * m2to4[3];
  2519. m4 := reciprocalN * m2to4[4];
  2520. skew := m3 / (sqrt(m2)*m2);
  2521. kurtosis := m4 / (m2 * m2);
  2522. end;
  2523. function norm(const data : array of Extended) : float; inline;
  2524. begin
  2525. norm:=Norm(PExtended(@data[0]),High(Data)+1);
  2526. end;
  2527. function norm(const data : PExtended; Const N : Integer) : float;
  2528. begin
  2529. norm:=sqrt(sumofsquares(data,N));
  2530. end;
  2531. {$endif FPC_HAS_TYPE_EXTENDED}
  2532. function MinIntValue(const Data: array of Integer): Integer;
  2533. var
  2534. I: SizeInt;
  2535. begin
  2536. Result := Data[Low(Data)];
  2537. For I := Succ(Low(Data)) To High(Data) Do
  2538. If Data[I] < Result Then Result := Data[I];
  2539. end;
  2540. function MaxIntValue(const Data: array of Integer): Integer;
  2541. var
  2542. I: SizeInt;
  2543. begin
  2544. Result := Data[Low(Data)];
  2545. For I := Succ(Low(Data)) To High(Data) Do
  2546. If Data[I] > Result Then Result := Data[I];
  2547. end;
  2548. function MinValue(const Data: array of Integer): Integer; inline;
  2549. begin
  2550. Result:=MinValue(Pinteger(@Data[0]),High(Data)+1)
  2551. end;
  2552. function MinValue(const Data: PInteger; Const N : Integer): Integer;
  2553. var
  2554. I: SizeInt;
  2555. begin
  2556. Result := Data[0];
  2557. For I := 1 To N-1 do
  2558. If Data[I] < Result Then Result := Data[I];
  2559. end;
  2560. function MaxValue(const Data: array of Integer): Integer; inline;
  2561. begin
  2562. Result:=MaxValue(PInteger(@Data[0]),High(Data)+1)
  2563. end;
  2564. function maxvalue(const data : PInteger; Const N : Integer) : Integer;
  2565. var
  2566. i : SizeInt;
  2567. begin
  2568. { get an initial value }
  2569. maxvalue:=data[0];
  2570. for i:=1 to N-1 do
  2571. if data[i]>maxvalue then
  2572. maxvalue:=data[i];
  2573. end;
  2574. {$ifdef FPC_HAS_TYPE_SINGLE}
  2575. function minvalue(const data : array of Single) : Single; inline;
  2576. begin
  2577. Result:=minvalue(PSingle(@data[0]),High(Data)+1);
  2578. end;
  2579. function minvalue(const data : PSingle; Const N : Integer) : Single;
  2580. var
  2581. i : SizeInt;
  2582. begin
  2583. { get an initial value }
  2584. minvalue:=data[0];
  2585. for i:=1 to N-1 do
  2586. if data[i]<minvalue then
  2587. minvalue:=data[i];
  2588. end;
  2589. function maxvalue(const data : array of Single) : Single; inline;
  2590. begin
  2591. Result:=maxvalue(PSingle(@data[0]),High(Data)+1);
  2592. end;
  2593. function maxvalue(const data : PSingle; Const N : Integer) : Single;
  2594. var
  2595. i : SizeInt;
  2596. begin
  2597. { get an initial value }
  2598. maxvalue:=data[0];
  2599. for i:=1 to N-1 do
  2600. if data[i]>maxvalue then
  2601. maxvalue:=data[i];
  2602. end;
  2603. {$endif FPC_HAS_TYPE_SINGLE}
  2604. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2605. function minvalue(const data : array of Double) : Double; inline;
  2606. begin
  2607. Result:=minvalue(PDouble(@data[0]),High(Data)+1);
  2608. end;
  2609. function minvalue(const data : PDouble; Const N : Integer) : Double;
  2610. var
  2611. i : SizeInt;
  2612. begin
  2613. { get an initial value }
  2614. minvalue:=data[0];
  2615. for i:=1 to N-1 do
  2616. if data[i]<minvalue then
  2617. minvalue:=data[i];
  2618. end;
  2619. function maxvalue(const data : array of Double) : Double; inline;
  2620. begin
  2621. Result:=maxvalue(PDouble(@data[0]),High(Data)+1);
  2622. end;
  2623. function maxvalue(const data : PDouble; Const N : Integer) : Double;
  2624. var
  2625. i : SizeInt;
  2626. begin
  2627. { get an initial value }
  2628. maxvalue:=data[0];
  2629. for i:=1 to N-1 do
  2630. if data[i]>maxvalue then
  2631. maxvalue:=data[i];
  2632. end;
  2633. {$endif FPC_HAS_TYPE_DOUBLE}
  2634. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2635. function minvalue(const data : array of Extended) : Extended; inline;
  2636. begin
  2637. Result:=minvalue(PExtended(@data[0]),High(Data)+1);
  2638. end;
  2639. function minvalue(const data : PExtended; Const N : Integer) : Extended;
  2640. var
  2641. i : SizeInt;
  2642. begin
  2643. { get an initial value }
  2644. minvalue:=data[0];
  2645. for i:=1 to N-1 do
  2646. if data[i]<minvalue then
  2647. minvalue:=data[i];
  2648. end;
  2649. function maxvalue(const data : array of Extended) : Extended; inline;
  2650. begin
  2651. Result:=maxvalue(PExtended(@data[0]),High(Data)+1);
  2652. end;
  2653. function maxvalue(const data : PExtended; Const N : Integer) : Extended;
  2654. var
  2655. i : SizeInt;
  2656. begin
  2657. { get an initial value }
  2658. maxvalue:=data[0];
  2659. for i:=1 to N-1 do
  2660. if data[i]>maxvalue then
  2661. maxvalue:=data[i];
  2662. end;
  2663. {$endif FPC_HAS_TYPE_EXTENDED}
  2664. function Min(a, b: Integer): Integer;inline;
  2665. begin
  2666. if a < b then
  2667. Result := a
  2668. else
  2669. Result := b;
  2670. end;
  2671. function Max(a, b: Integer): Integer;inline;
  2672. begin
  2673. if a > b then
  2674. Result := a
  2675. else
  2676. Result := b;
  2677. end;
  2678. {
  2679. function Min(a, b: Cardinal): Cardinal;inline;
  2680. begin
  2681. if a < b then
  2682. Result := a
  2683. else
  2684. Result := b;
  2685. end;
  2686. function Max(a, b: Cardinal): Cardinal;inline;
  2687. begin
  2688. if a > b then
  2689. Result := a
  2690. else
  2691. Result := b;
  2692. end;
  2693. }
  2694. function Min(a, b: Int64): Int64;inline;
  2695. begin
  2696. if a < b then
  2697. Result := a
  2698. else
  2699. Result := b;
  2700. end;
  2701. function Max(a, b: Int64): Int64;inline;
  2702. begin
  2703. if a > b then
  2704. Result := a
  2705. else
  2706. Result := b;
  2707. end;
  2708. function Min(a, b: QWord): QWord; inline;
  2709. begin
  2710. if a < b then
  2711. Result := a
  2712. else
  2713. Result := b;
  2714. end;
  2715. function Max(a, b: QWord): Qword;inline;
  2716. begin
  2717. if a > b then
  2718. Result := a
  2719. else
  2720. Result := b;
  2721. end;
  2722. {$ifdef FPC_HAS_TYPE_SINGLE}
  2723. function Min(a, b: Single): Single;inline;
  2724. begin
  2725. if a < b then
  2726. Result := a
  2727. else
  2728. Result := b;
  2729. end;
  2730. function Max(a, b: Single): Single;inline;
  2731. begin
  2732. if a > b then
  2733. Result := a
  2734. else
  2735. Result := b;
  2736. end;
  2737. {$endif FPC_HAS_TYPE_SINGLE}
  2738. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2739. function Min(a, b: Double): Double;inline;
  2740. begin
  2741. if a < b then
  2742. Result := a
  2743. else
  2744. Result := b;
  2745. end;
  2746. function Max(a, b: Double): Double;inline;
  2747. begin
  2748. if a > b then
  2749. Result := a
  2750. else
  2751. Result := b;
  2752. end;
  2753. {$endif FPC_HAS_TYPE_DOUBLE}
  2754. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2755. function Min(a, b: Extended): Extended;inline;
  2756. begin
  2757. if a < b then
  2758. Result := a
  2759. else
  2760. Result := b;
  2761. end;
  2762. function Max(a, b: Extended): Extended;inline;
  2763. begin
  2764. if a > b then
  2765. Result := a
  2766. else
  2767. Result := b;
  2768. end;
  2769. {$endif FPC_HAS_TYPE_EXTENDED}
  2770. function InRange(const AValue, AMin, AMax: Integer): Boolean;inline;
  2771. begin
  2772. Result:=(AValue>=AMin) and (AValue<=AMax);
  2773. end;
  2774. function InRange(const AValue, AMin, AMax: Int64): Boolean;inline;
  2775. begin
  2776. Result:=(AValue>=AMin) and (AValue<=AMax);
  2777. end;
  2778. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2779. function InRange(const AValue, AMin, AMax: Double): Boolean;inline;
  2780. begin
  2781. Result:=(AValue>=AMin) and (AValue<=AMax);
  2782. end;
  2783. {$endif FPC_HAS_TYPE_DOUBLE}
  2784. function EnsureRange(const AValue, AMin, AMax: Integer): Integer;inline;
  2785. begin
  2786. Result:=AValue;
  2787. If Result<AMin then
  2788. Result:=AMin;
  2789. if Result>AMax then
  2790. Result:=AMax;
  2791. end;
  2792. function EnsureRange(const AValue, AMin, AMax: Int64): Int64;inline;
  2793. begin
  2794. Result:=AValue;
  2795. If Result<AMin then
  2796. Result:=AMin;
  2797. if Result>AMax then
  2798. Result:=AMax;
  2799. end;
  2800. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2801. function EnsureRange(const AValue, AMin, AMax: Double): Double;inline;
  2802. begin
  2803. Result:=AValue;
  2804. If Result<AMin then
  2805. Result:=AMin;
  2806. if Result>AMax then
  2807. Result:=AMax;
  2808. end;
  2809. {$endif FPC_HAS_TYPE_DOUBLE}
  2810. Const
  2811. EZeroResolution = Extended(1E-16);
  2812. DZeroResolution = Double(1E-12);
  2813. SZeroResolution = Single(1E-4);
  2814. function IsZero(const A: Single; Epsilon: Single): Boolean;
  2815. begin
  2816. if (Epsilon=0) then
  2817. Epsilon:=SZeroResolution;
  2818. Result:=Abs(A)<=Epsilon;
  2819. end;
  2820. function IsZero(const A: Single): Boolean;inline;
  2821. begin
  2822. Result:=IsZero(A,single(SZeroResolution));
  2823. end;
  2824. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2825. function IsZero(const A: Double; Epsilon: Double): Boolean;
  2826. begin
  2827. if (Epsilon=0) then
  2828. Epsilon:=DZeroResolution;
  2829. Result:=Abs(A)<=Epsilon;
  2830. end;
  2831. function IsZero(const A: Double): Boolean;inline;
  2832. begin
  2833. Result:=IsZero(A,DZeroResolution);
  2834. end;
  2835. {$endif FPC_HAS_TYPE_DOUBLE}
  2836. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2837. function IsZero(const A: Extended; Epsilon: Extended): Boolean;
  2838. begin
  2839. if (Epsilon=0) then
  2840. Epsilon:=EZeroResolution;
  2841. Result:=Abs(A)<=Epsilon;
  2842. end;
  2843. function IsZero(const A: Extended): Boolean;inline;
  2844. begin
  2845. Result:=IsZero(A,EZeroResolution);
  2846. end;
  2847. {$endif FPC_HAS_TYPE_EXTENDED}
  2848. type
  2849. TSplitDouble = packed record
  2850. cards: Array[0..1] of cardinal;
  2851. end;
  2852. TSplitExtended = packed record
  2853. cards: Array[0..1] of cardinal;
  2854. w: word;
  2855. end;
  2856. function IsNan(const d : Single): Boolean; overload;
  2857. begin
  2858. result:=(longword(d) and $7fffffff)>$7f800000;
  2859. end;
  2860. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2861. function IsNan(const d : Double): Boolean;
  2862. var
  2863. fraczero, expMaximal: boolean;
  2864. begin
  2865. {$if defined(FPC_BIG_ENDIAN) or defined(FPC_DOUBLE_HILO_SWAPPED)}
  2866. expMaximal := ((TSplitDouble(d).cards[0] shr 20) and $7ff) = 2047;
  2867. fraczero:= (TSplitDouble(d).cards[0] and $fffff = 0) and
  2868. (TSplitDouble(d).cards[1] = 0);
  2869. {$else FPC_BIG_ENDIAN}
  2870. expMaximal := ((TSplitDouble(d).cards[1] shr 20) and $7ff) = 2047;
  2871. fraczero := (TSplitDouble(d).cards[1] and $fffff = 0) and
  2872. (TSplitDouble(d).cards[0] = 0);
  2873. {$endif FPC_BIG_ENDIAN}
  2874. Result:=expMaximal and not(fraczero);
  2875. end;
  2876. {$endif FPC_HAS_TYPE_DOUBLE}
  2877. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2878. function IsNan(const d : Extended): Boolean; overload;
  2879. var
  2880. fraczero, expMaximal: boolean;
  2881. begin
  2882. {$ifdef FPC_BIG_ENDIAN}
  2883. {$error no support for big endian extended type yet}
  2884. {$else FPC_BIG_ENDIAN}
  2885. expMaximal := (TSplitExtended(d).w and $7fff) = 32767;
  2886. fraczero := (TSplitExtended(d).cards[0] = 0) and
  2887. ((TSplitExtended(d).cards[1] and $7fffffff) = 0);
  2888. {$endif FPC_BIG_ENDIAN}
  2889. Result:=expMaximal and not(fraczero);
  2890. end;
  2891. {$endif FPC_HAS_TYPE_EXTENDED}
  2892. function IsInfinite(const d : Single): Boolean; overload;
  2893. begin
  2894. result:=(longword(d) and $7fffffff)=$7f800000;
  2895. end;
  2896. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2897. function IsInfinite(const d : Double): Boolean; overload;
  2898. var
  2899. fraczero, expMaximal: boolean;
  2900. begin
  2901. {$if defined(FPC_BIG_ENDIAN) or defined(FPC_DOUBLE_HILO_SWAPPED)}
  2902. expMaximal := ((TSplitDouble(d).cards[0] shr 20) and $7ff) = 2047;
  2903. fraczero:= (TSplitDouble(d).cards[0] and $fffff = 0) and
  2904. (TSplitDouble(d).cards[1] = 0);
  2905. {$else FPC_BIG_ENDIAN}
  2906. expMaximal := ((TSplitDouble(d).cards[1] shr 20) and $7ff) = 2047;
  2907. fraczero := (TSplitDouble(d).cards[1] and $fffff = 0) and
  2908. (TSplitDouble(d).cards[0] = 0);
  2909. {$endif FPC_BIG_ENDIAN}
  2910. Result:=expMaximal and fraczero;
  2911. end;
  2912. {$endif FPC_HAS_TYPE_DOUBLE}
  2913. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2914. function IsInfinite(const d : Extended): Boolean; overload;
  2915. var
  2916. fraczero, expMaximal: boolean;
  2917. begin
  2918. {$ifdef FPC_BIG_ENDIAN}
  2919. {$error no support for big endian extended type yet}
  2920. {$else FPC_BIG_ENDIAN}
  2921. expMaximal := (TSplitExtended(d).w and $7fff) = 32767;
  2922. fraczero := (TSplitExtended(d).cards[0] = 0) and
  2923. ((TSplitExtended(d).cards[1] and $7fffffff) = 0);
  2924. {$endif FPC_BIG_ENDIAN}
  2925. Result:=expMaximal and fraczero;
  2926. end;
  2927. {$endif FPC_HAS_TYPE_EXTENDED}
  2928. function copysign(x,y: float): float;
  2929. begin
  2930. {$if defined(FPC_HAS_TYPE_FLOAT128)}
  2931. {$error copysign not yet implemented for float128}
  2932. {$elseif defined(FPC_HAS_TYPE_EXTENDED)}
  2933. TSplitExtended(x).w:=(TSplitExtended(x).w and $7fff) or (TSplitExtended(y).w and $8000);
  2934. {$elseif defined(FPC_HAS_TYPE_DOUBLE)}
  2935. {$if defined(FPC_BIG_ENDIAN) or defined(FPC_DOUBLE_HILO_SWAPPED)}
  2936. TSplitDouble(x).cards[0]:=(TSplitDouble(x).cards[0] and $7fffffff) or (TSplitDouble(y).cards[0] and longword($80000000));
  2937. {$else}
  2938. TSplitDouble(x).cards[1]:=(TSplitDouble(x).cards[1] and $7fffffff) or (TSplitDouble(y).cards[1] and longword($80000000));
  2939. {$endif}
  2940. {$else}
  2941. longword(x):=longword(x and $7fffffff) or (longword(y) and longword($80000000));
  2942. {$endif}
  2943. result:=x;
  2944. end;
  2945. {$ifdef FPC_HAS_TYPE_EXTENDED}
  2946. function SameValue(const A, B: Extended; Epsilon: Extended): Boolean;
  2947. begin
  2948. if (Epsilon=0) then
  2949. Epsilon:=Max(Min(Abs(A),Abs(B))*EZeroResolution,EZeroResolution);
  2950. if (A>B) then
  2951. Result:=((A-B)<=Epsilon)
  2952. else
  2953. Result:=((B-A)<=Epsilon);
  2954. end;
  2955. function SameValue(const A, B: Extended): Boolean;inline;
  2956. begin
  2957. Result:=SameValue(A,B,0.0);
  2958. end;
  2959. {$endif FPC_HAS_TYPE_EXTENDED}
  2960. {$ifdef FPC_HAS_TYPE_DOUBLE}
  2961. function SameValue(const A, B: Double): Boolean;inline;
  2962. begin
  2963. Result:=SameValue(A,B,0.0);
  2964. end;
  2965. function SameValue(const A, B: Double; Epsilon: Double): Boolean;
  2966. begin
  2967. if (Epsilon=0) then
  2968. Epsilon:=Max(Min(Abs(A),Abs(B))*DZeroResolution,DZeroResolution);
  2969. if (A>B) then
  2970. Result:=((A-B)<=Epsilon)
  2971. else
  2972. Result:=((B-A)<=Epsilon);
  2973. end;
  2974. {$endif FPC_HAS_TYPE_DOUBLE}
  2975. function SameValue(const A, B: Single): Boolean;inline;
  2976. begin
  2977. Result:=SameValue(A,B,0);
  2978. end;
  2979. function SameValue(const A, B: Single; Epsilon: Single): Boolean;
  2980. begin
  2981. if (Epsilon=0) then
  2982. Epsilon:=Max(Min(Abs(A),Abs(B))*SZeroResolution,SZeroResolution);
  2983. if (A>B) then
  2984. Result:=((A-B)<=Epsilon)
  2985. else
  2986. Result:=((B-A)<=Epsilon);
  2987. end;
  2988. // Some CPUs probably allow a faster way of doing this in a single operation...
  2989. // There weshould define FPC_MATH_HAS_CPUDIVMOD in the header mathuh.inc and implement it using asm.
  2990. {$ifndef FPC_MATH_HAS_DIVMOD}
  2991. procedure DivMod(Dividend: LongInt; Divisor: Word; var Result, Remainder: Word);
  2992. begin
  2993. if Dividend < 0 then
  2994. begin
  2995. { Use DivMod with >=0 dividend }
  2996. Dividend:=-Dividend;
  2997. { The documented behavior of Pascal's div/mod operators and DivMod
  2998. on negative dividends is to return Result closer to zero and
  2999. a negative Remainder. Which means that we can just negate both
  3000. Result and Remainder, and all it's Ok. }
  3001. Result:=-(Dividend Div Divisor);
  3002. Remainder:=-(Dividend+(Result*Divisor));
  3003. end
  3004. else
  3005. begin
  3006. Result:=Dividend Div Divisor;
  3007. Remainder:=Dividend-(Result*Divisor);
  3008. end;
  3009. end;
  3010. procedure DivMod(Dividend: LongInt; Divisor: Word; var Result, Remainder: SmallInt);
  3011. begin
  3012. if Dividend < 0 then
  3013. begin
  3014. { Use DivMod with >=0 dividend }
  3015. Dividend:=-Dividend;
  3016. { The documented behavior of Pascal's div/mod operators and DivMod
  3017. on negative dividends is to return Result closer to zero and
  3018. a negative Remainder. Which means that we can just negate both
  3019. Result and Remainder, and all it's Ok. }
  3020. Result:=-(Dividend Div Divisor);
  3021. Remainder:=-(Dividend+(Result*Divisor));
  3022. end
  3023. else
  3024. begin
  3025. Result:=Dividend Div Divisor;
  3026. Remainder:=Dividend-(Result*Divisor);
  3027. end;
  3028. end;
  3029. procedure DivMod(Dividend: DWord; Divisor: DWord; var Result, Remainder: DWord);
  3030. begin
  3031. Result:=Dividend Div Divisor;
  3032. Remainder:=Dividend-(Result*Divisor);
  3033. end;
  3034. procedure DivMod(Dividend: LongInt; Divisor: LongInt; var Result, Remainder: LongInt);
  3035. begin
  3036. if Dividend < 0 then
  3037. begin
  3038. { Use DivMod with >=0 dividend }
  3039. Dividend:=-Dividend;
  3040. { The documented behavior of Pascal's div/mod operators and DivMod
  3041. on negative dividends is to return Result closer to zero and
  3042. a negative Remainder. Which means that we can just negate both
  3043. Result and Remainder, and all it's Ok. }
  3044. Result:=-(Dividend Div Divisor);
  3045. Remainder:=-(Dividend+(Result*Divisor));
  3046. end
  3047. else
  3048. begin
  3049. Result:=Dividend Div Divisor;
  3050. Remainder:=Dividend-(Result*Divisor);
  3051. end;
  3052. end;
  3053. {$endif FPC_MATH_HAS_DIVMOD}
  3054. { Floating point modulo}
  3055. {$ifdef FPC_HAS_TYPE_SINGLE}
  3056. function FMod(const a, b: Single): Single;inline;overload;
  3057. begin
  3058. result:= a-b * Int(a/b);
  3059. end;
  3060. {$endif FPC_HAS_TYPE_SINGLE}
  3061. {$ifdef FPC_HAS_TYPE_DOUBLE}
  3062. function FMod(const a, b: Double): Double;inline;overload;
  3063. begin
  3064. result:= a-b * Int(a/b);
  3065. end;
  3066. {$endif FPC_HAS_TYPE_DOUBLE}
  3067. {$ifdef FPC_HAS_TYPE_EXTENDED}
  3068. function FMod(const a, b: Extended): Extended;inline;overload;
  3069. begin
  3070. result:= a-b * Int(a/b);
  3071. end;
  3072. {$endif FPC_HAS_TYPE_EXTENDED}
  3073. operator mod(const a,b:float) c:float;inline;
  3074. begin
  3075. c:= a-b * Int(a/b);
  3076. if SameValue(abs(c),abs(b)) then
  3077. c:=0.0;
  3078. end;
  3079. function ifthen(val:boolean;const iftrue:integer; const iffalse:integer= 0) :integer;
  3080. begin
  3081. if val then result:=iftrue else result:=iffalse;
  3082. end;
  3083. function ifthen(val:boolean;const iftrue:int64 ; const iffalse:int64 = 0) :int64;
  3084. begin
  3085. if val then result:=iftrue else result:=iffalse;
  3086. end;
  3087. function ifthen(val:boolean;const iftrue:double ; const iffalse:double =0.0):double;
  3088. begin
  3089. if val then result:=iftrue else result:=iffalse;
  3090. end;
  3091. // dilemma here. asm can do the two comparisons in one go?
  3092. // but pascal is portable and can be inlined. Ah well, we need purepascal's anyway:
  3093. function CompareValue(const A, B : Integer): TValueRelationship;
  3094. begin
  3095. result:=GreaterThanValue;
  3096. if a=b then
  3097. result:=EqualsValue
  3098. else
  3099. if a<b then
  3100. result:=LessThanValue;
  3101. end;
  3102. function CompareValue(const A, B: Int64): TValueRelationship;
  3103. begin
  3104. result:=GreaterThanValue;
  3105. if a=b then
  3106. result:=EqualsValue
  3107. else
  3108. if a<b then
  3109. result:=LessThanValue;
  3110. end;
  3111. function CompareValue(const A, B: QWord): TValueRelationship;
  3112. begin
  3113. result:=GreaterThanValue;
  3114. if a=b then
  3115. result:=EqualsValue
  3116. else
  3117. if a<b then
  3118. result:=LessThanValue;
  3119. end;
  3120. {$ifdef FPC_HAS_TYPE_SINGLE}
  3121. function CompareValue(const A, B: Single; delta: Single = 0.0): TValueRelationship;
  3122. begin
  3123. result:=GreaterThanValue;
  3124. if abs(a-b)<=delta then
  3125. result:=EqualsValue
  3126. else
  3127. if a<b then
  3128. result:=LessThanValue;
  3129. end;
  3130. {$endif}
  3131. {$ifdef FPC_HAS_TYPE_DOUBLE}
  3132. function CompareValue(const A, B: Double; delta: Double = 0.0): TValueRelationship;
  3133. begin
  3134. result:=GreaterThanValue;
  3135. if abs(a-b)<=delta then
  3136. result:=EqualsValue
  3137. else
  3138. if a<b then
  3139. result:=LessThanValue;
  3140. end;
  3141. {$endif}
  3142. {$ifdef FPC_HAS_TYPE_EXTENDED}
  3143. function CompareValue (const A, B: Extended; delta: Extended = 0.0): TValueRelationship;
  3144. begin
  3145. result:=GreaterThanValue;
  3146. if abs(a-b)<=delta then
  3147. result:=EqualsValue
  3148. else
  3149. if a<b then
  3150. result:=LessThanValue;
  3151. end;
  3152. {$endif}
  3153. {$ifdef FPC_HAS_TYPE_DOUBLE}
  3154. function RoundTo(const AValue: Double; const Digits: TRoundToRange): Double;
  3155. var
  3156. RV : Double;
  3157. begin
  3158. RV:=IntPower(10,Digits);
  3159. Result:=Round(AValue/RV)*RV;
  3160. end;
  3161. {$endif}
  3162. {$ifdef FPC_HAS_TYPE_EXTENDED}
  3163. function RoundTo(const AVAlue: Extended; const Digits: TRoundToRange): Extended;
  3164. var
  3165. RV : Extended;
  3166. begin
  3167. RV:=IntPower(10,Digits);
  3168. Result:=Round(AValue/RV)*RV;
  3169. end;
  3170. {$endif}
  3171. {$ifdef FPC_HAS_TYPE_SINGLE}
  3172. function RoundTo(const AValue: Single; const Digits: TRoundToRange): Single;
  3173. var
  3174. RV : Single;
  3175. begin
  3176. RV:=IntPower(10,Digits);
  3177. Result:=Round(AValue/RV)*RV;
  3178. end;
  3179. {$endif}
  3180. {$ifdef FPC_HAS_TYPE_SINGLE}
  3181. function SimpleRoundTo(const AValue: Single; const Digits: TRoundToRange = -2): Single;
  3182. var
  3183. RV : Single;
  3184. begin
  3185. RV := IntPower(10, -Digits);
  3186. if AValue < 0 then
  3187. Result := Int((AValue*RV) - 0.5)/RV
  3188. else
  3189. Result := Int((AValue*RV) + 0.5)/RV;
  3190. end;
  3191. {$endif}
  3192. {$ifdef FPC_HAS_TYPE_DOUBLE}
  3193. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange = -2): Double;
  3194. var
  3195. RV : Double;
  3196. begin
  3197. RV := IntPower(10, -Digits);
  3198. if AValue < 0 then
  3199. Result := Int((AValue*RV) - 0.5)/RV
  3200. else
  3201. Result := Int((AValue*RV) + 0.5)/RV;
  3202. end;
  3203. {$endif}
  3204. {$ifdef FPC_HAS_TYPE_EXTENDED}
  3205. function SimpleRoundTo(const AValue: Extended; const Digits: TRoundToRange = -2): Extended;
  3206. var
  3207. RV : Extended;
  3208. begin
  3209. RV := IntPower(10, -Digits);
  3210. if AValue < 0 then
  3211. Result := Int((AValue*RV) - 0.5)/RV
  3212. else
  3213. Result := Int((AValue*RV) + 0.5)/RV;
  3214. end;
  3215. {$endif}
  3216. function RandomFrom(const AValues: array of Double): Double; overload;
  3217. begin
  3218. result:=AValues[random(High(AValues)+1)];
  3219. end;
  3220. function RandomFrom(const AValues: array of Integer): Integer; overload;
  3221. begin
  3222. result:=AValues[random(High(AValues)+1)];
  3223. end;
  3224. function RandomFrom(const AValues: array of Int64): Int64; overload;
  3225. begin
  3226. result:=AValues[random(High(AValues)+1)];
  3227. end;
  3228. {$if FPC_FULLVERSION >=30101}
  3229. generic function RandomFrom<T>(const AValues:array of T):T;
  3230. begin
  3231. result:=AValues[random(High(AValues)+1)];
  3232. end;
  3233. {$endif}
  3234. function FutureValue(ARate: Float; NPeriods: Integer;
  3235. APayment, APresentValue: Float; APaymentTime: TPaymentTime): Float;
  3236. var
  3237. q, qn, factor: Float;
  3238. begin
  3239. if ARate = 0 then
  3240. Result := -APresentValue - APayment * NPeriods
  3241. else begin
  3242. q := 1.0 + ARate;
  3243. qn := power(q, NPeriods);
  3244. factor := (qn - 1) / (q - 1);
  3245. if APaymentTime = ptStartOfPeriod then
  3246. factor := factor * q;
  3247. Result := -(APresentValue * qn + APayment*factor);
  3248. end;
  3249. end;
  3250. function InterestRate(NPeriods: Integer; APayment, APresentValue, AFutureValue: Float;
  3251. APaymentTime: TPaymentTime): Float;
  3252. { The interest rate cannot be calculated analytically. We solve the equation
  3253. numerically by means of the Newton method:
  3254. - guess value for the interest reate
  3255. - calculate at which interest rate the tangent of the curve fv(rate)
  3256. (straight line!) has the requested future vale.
  3257. - use this rate for the next iteration. }
  3258. const
  3259. DELTA = 0.001;
  3260. EPS = 1E-9; // required precision of interest rate (after typ. 6 iterations)
  3261. MAXIT = 20; // max iteration count to protect agains non-convergence
  3262. var
  3263. r1, r2, dr: Float;
  3264. fv1, fv2: Float;
  3265. iteration: Integer;
  3266. begin
  3267. iteration := 0;
  3268. r1 := 0.05; // inital guess
  3269. repeat
  3270. r2 := r1 + DELTA;
  3271. fv1 := FutureValue(r1, NPeriods, APayment, APresentValue, APaymentTime);
  3272. fv2 := FutureValue(r2, NPeriods, APayment, APresentValue, APaymentTime);
  3273. dr := (AFutureValue - fv1) / (fv2 - fv1) * delta; // tangent at fv(r)
  3274. r1 := r1 + dr; // next guess
  3275. inc(iteration);
  3276. until (abs(dr) < EPS) or (iteration >= MAXIT);
  3277. Result := r1;
  3278. end;
  3279. function NumberOfPeriods(ARate, APayment, APresentValue, AFutureValue: Float;
  3280. APaymentTime: TPaymentTime): Float;
  3281. { Solve the cash flow equation (1) for q^n and take the logarithm }
  3282. var
  3283. q, x1, x2: Float;
  3284. begin
  3285. if ARate = 0 then
  3286. Result := -(APresentValue + AFutureValue) / APayment
  3287. else begin
  3288. q := 1.0 + ARate;
  3289. if APaymentTime = ptStartOfPeriod then
  3290. APayment := APayment * q;
  3291. x1 := APayment - AFutureValue * ARate;
  3292. x2 := APayment + APresentValue * ARate;
  3293. if (x2 = 0) // we have to divide by x2
  3294. or (sign(x1) * sign(x2) < 0) // the argument of the log is negative
  3295. then
  3296. Result := Infinity
  3297. else begin
  3298. Result := ln(x1/x2) / ln(q);
  3299. end;
  3300. end;
  3301. end;
  3302. function Payment(ARate: Float; NPeriods: Integer;
  3303. APresentValue, AFutureValue: Float; APaymentTime: TPaymentTime): Float;
  3304. var
  3305. q, qn, factor: Float;
  3306. begin
  3307. if ARate = 0 then
  3308. Result := -(AFutureValue + APresentValue) / NPeriods
  3309. else begin
  3310. q := 1.0 + ARate;
  3311. qn := power(q, NPeriods);
  3312. factor := (qn - 1) / (q - 1);
  3313. if APaymentTime = ptStartOfPeriod then
  3314. factor := factor * q;
  3315. Result := -(AFutureValue + APresentValue * qn) / factor;
  3316. end;
  3317. end;
  3318. function PresentValue(ARate: Float; NPeriods: Integer;
  3319. APayment, AFutureValue: Float; APaymentTime: TPaymentTime): Float;
  3320. var
  3321. q, qn, factor: Float;
  3322. begin
  3323. if ARate = 0.0 then
  3324. Result := -AFutureValue - APayment * NPeriods
  3325. else begin
  3326. q := 1.0 + ARate;
  3327. qn := power(q, NPeriods);
  3328. factor := (qn - 1) / (q - 1);
  3329. if APaymentTime = ptStartOfPeriod then
  3330. factor := factor * q;
  3331. Result := -(AFutureValue + APayment*factor) / qn;
  3332. end;
  3333. end;
  3334. {$else}
  3335. implementation
  3336. {$endif FPUNONE}
  3337. end.