genmath.inc 43 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343
  1. {
  2. $Id$
  3. This file is part of the Free Pascal run time library.
  4. Copyright (c) 1999-2001 by Several contributors
  5. Generic mathemtical routines (on type real)
  6. See the file COPYING.FPC, included in this distribution,
  7. for details about the copyright.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  11. **********************************************************************}
  12. {*************************************************************************}
  13. { Credits }
  14. {*************************************************************************}
  15. { Copyright Abandoned, 1987, Fred Fish }
  16. { }
  17. { This previously copyrighted work has been placed into the }
  18. { public domain by the author (Fred Fish) and may be freely used }
  19. { for any purpose, private or commercial. I would appreciate }
  20. { it, as a courtesy, if this notice is left in all copies and }
  21. { derivative works. Thank you, and enjoy... }
  22. { }
  23. { The author makes no warranty of any kind with respect to this }
  24. { product and explicitly disclaims any implied warranties of }
  25. { merchantability or fitness for any particular purpose. }
  26. {-------------------------------------------------------------------------}
  27. { Copyright (c) 1992 Odent Jean Philippe }
  28. { }
  29. { The source can be modified as long as my name appears and some }
  30. { notes explaining the modifications done are included in the file. }
  31. {-------------------------------------------------------------------------}
  32. { Copyright (c) 1997 Carl Eric Codere }
  33. {-------------------------------------------------------------------------}
  34. {$goto on}
  35. type
  36. TabCoef = array[0..6] of Real;
  37. const
  38. PIO2 = 1.57079632679489661923; { pi/2 }
  39. PIO4 = 7.85398163397448309616E-1; { pi/4 }
  40. SQRT2 = 1.41421356237309504880; { sqrt(2) }
  41. SQRTH = 7.07106781186547524401E-1; { sqrt(2)/2 }
  42. LOG2E = 1.4426950408889634073599; { 1/log(2) }
  43. SQ2OPI = 7.9788456080286535587989E-1; { sqrt( 2/pi )}
  44. LOGE2 = 6.93147180559945309417E-1; { log(2) }
  45. LOGSQ2 = 3.46573590279972654709E-1; { log(2)/2 }
  46. THPIO4 = 2.35619449019234492885; { 3*pi/4 }
  47. TWOOPI = 6.36619772367581343075535E-1; { 2/pi }
  48. lossth = 1.073741824e9;
  49. MAXLOG = 8.8029691931113054295988E1; { log(2**127) }
  50. MINLOG = -8.872283911167299960540E1; { log(2**-128) }
  51. DP1 = 7.85398125648498535156E-1;
  52. DP2 = 3.77489470793079817668E-8;
  53. DP3 = 2.69515142907905952645E-15;
  54. const sincof : TabCoef = (
  55. 1.58962301576546568060E-10,
  56. -2.50507477628578072866E-8,
  57. 2.75573136213857245213E-6,
  58. -1.98412698295895385996E-4,
  59. 8.33333333332211858878E-3,
  60. -1.66666666666666307295E-1, 0);
  61. coscof : TabCoef = (
  62. -1.13585365213876817300E-11,
  63. 2.08757008419747316778E-9,
  64. -2.75573141792967388112E-7,
  65. 2.48015872888517045348E-5,
  66. -1.38888888888730564116E-3,
  67. 4.16666666666665929218E-2, 0);
  68. { also necessary for Int() on systems with 64bit floats (JM) }
  69. type
  70. {$ifdef ENDIAN_LITTLE}
  71. float64 = packed record
  72. low: longint;
  73. high: longint;
  74. end;
  75. {$else}
  76. float64 = packed record
  77. high: longint;
  78. low: longint;
  79. end;
  80. {$endif}
  81. {$ifndef FPC_SYSTEM_HAS_TRUNC}
  82. type
  83. float32 = longint;
  84. flag = byte;
  85. Function extractFloat64Frac0(const a: float64): longint;
  86. Begin
  87. extractFloat64Frac0 := a.high and $000FFFFF;
  88. End;
  89. Function extractFloat64Frac1(const a: float64): longint;
  90. Begin
  91. extractFloat64Frac1 := a.low;
  92. End;
  93. Function extractFloat64Exp(const a: float64): smallint;
  94. Begin
  95. extractFloat64Exp:= ( a.high shr 20 ) AND $7FF;
  96. End;
  97. Function extractFloat64Frac(const a: float64): int64;
  98. Begin
  99. extractFloat64Frac:=int64(a) and $000FFFFFFFFFFFFF;
  100. End;
  101. Function extractFloat64Sign(const a: float64) : flag;
  102. Begin
  103. extractFloat64Sign := a.high shr 31;
  104. End;
  105. Procedure shortShift64Left(a0:longint; a1:longint; count:smallint; VAR z0Ptr:longint; VAR z1Ptr:longint );
  106. Begin
  107. z1Ptr := a1 shl count;
  108. if count = 0 then
  109. z0Ptr := a0
  110. else
  111. z0Ptr := ( a0 shl count ) OR ( a1 shr ( ( - count ) AND 31 ) );
  112. End;
  113. function float64_to_int32_round_to_zero(a: float64 ): longint;
  114. Var
  115. aSign: flag;
  116. aExp, shiftCount: smallint;
  117. aSig0, aSig1, absZ, aSigExtra: longint;
  118. z: longint;
  119. Begin
  120. aSig1 := extractFloat64Frac1( a );
  121. aSig0 := extractFloat64Frac0( a );
  122. aExp := extractFloat64Exp( a );
  123. aSign := extractFloat64Sign( a );
  124. shiftCount := aExp - $413;
  125. if 0<=shiftCount then
  126. Begin
  127. if (aExp=$7FF) and ((aSig0 or aSig1)<>0) then
  128. HandleError(207);
  129. shortShift64Left(aSig0 OR $00100000, aSig1, shiftCount, absZ, aSigExtra );
  130. End
  131. else
  132. Begin
  133. if aExp<$3FF then
  134. begin
  135. float64_to_int32_round_to_zero := 0;
  136. exit;
  137. end;
  138. aSig0 := aSig0 or $00100000;
  139. aSigExtra := ( aSig0 shl ( shiftCount and 31 ) ) OR aSig1;
  140. absZ := aSig0 shr ( - shiftCount );
  141. End;
  142. if aSign<>0 then
  143. z:=-absZ
  144. else
  145. z:=absZ;
  146. if ((aSign<>0) xor (z<0)) AND (z<>0) then
  147. HandleError(207);
  148. float64_to_int32_round_to_zero := z;
  149. End;
  150. {$ifndef VER1_0}
  151. function float64_to_int64_round_to_zero(a : float64) : int64;
  152. var
  153. aSign : flag;
  154. aExp, shiftCount : smallint;
  155. aSig : int64;
  156. z : int64;
  157. begin
  158. aSig:=extractFloat64Frac(a);
  159. aExp:=extractFloat64Exp(a);
  160. aSign:=extractFloat64Sign(a);
  161. if aExp<>0 then
  162. aSig:=aSig or $0010000000000000;
  163. shiftCount:= aExp-$433;
  164. if 0<=shiftCount then
  165. begin
  166. if aExp>=$43e then
  167. begin
  168. if int64(a)<>$C3E0000000000000 then
  169. HandleError(207);
  170. { pascal doesn't know Inf for int64 }
  171. HandleError(207);
  172. end;
  173. z:=aSig shl shiftCount;
  174. end
  175. else
  176. begin
  177. if aExp<$3fe then
  178. begin
  179. result:=0;
  180. exit;
  181. end;
  182. z:=aSig shr -shiftCount;
  183. {
  184. if (aSig shl (shiftCount and 63))<>0 then
  185. float_exception_flags |= float_flag_inexact;
  186. }
  187. end;
  188. if aSign<>0 then
  189. z:=-z;
  190. result:=z;
  191. end;
  192. {$endif VER1_0}
  193. Function ExtractFloat32Frac(a : Float32) : longint;
  194. Begin
  195. ExtractFloat32Frac := A AND $007FFFFF;
  196. End;
  197. Function extractFloat32Exp( a: float32 ): smallint;
  198. Begin
  199. extractFloat32Exp := (a shr 23) AND $FF;
  200. End;
  201. Function extractFloat32Sign( a: float32 ): Flag;
  202. Begin
  203. extractFloat32Sign := a shr 31;
  204. End;
  205. Function float32_to_int32_round_to_zero( a: Float32 ): longint;
  206. Var
  207. aSign : flag;
  208. aExp, shiftCount : smallint;
  209. aSig : longint;
  210. z : longint;
  211. Begin
  212. aSig := extractFloat32Frac( a );
  213. aExp := extractFloat32Exp( a );
  214. aSign := extractFloat32Sign( a );
  215. shiftCount := aExp - $9E;
  216. if ( 0 <= shiftCount ) then
  217. Begin
  218. if ( a <> Float32($CF000000) ) then
  219. Begin
  220. if ( (aSign=0) or ( ( aExp = $FF ) and (aSig<>0) ) ) then
  221. Begin
  222. HandleError(207);
  223. exit;
  224. end;
  225. End;
  226. HandleError(207);
  227. exit;
  228. End
  229. else
  230. if ( aExp <= $7E ) then
  231. Begin
  232. float32_to_int32_round_to_zero := 0;
  233. exit;
  234. End;
  235. aSig := ( aSig or $00800000 ) shl 8;
  236. z := aSig shr ( - shiftCount );
  237. if ( aSign<>0 ) then z := - z;
  238. float32_to_int32_round_to_zero := z;
  239. End;
  240. function trunc(d : real) : int64;[internconst:in_const_trunc];
  241. var
  242. {$ifdef cpuarm}
  243. l: longint;
  244. {$endif cpuarm}
  245. f32 : float32;
  246. f64 : float64;
  247. Begin
  248. { in emulation mode the real is equal to a single }
  249. { otherwise in fpu mode, it is equal to a double }
  250. { extended is not supported yet. }
  251. if sizeof(D) > 8 then
  252. HandleError(255);
  253. if sizeof(D)=8 then
  254. begin
  255. move(d,f64,sizeof(f64));
  256. {$ifdef cpuarm}
  257. { the arm fpu has a strange opinion how a double has to be stored }
  258. l:=f64.low;
  259. f64.low:=f64.high;
  260. f64.high:=l;
  261. {$endif cpuarm}
  262. {$ifdef VER1_0}
  263. trunc:=float64_to_int32_round_to_zero(f64);
  264. {$else VER1_0}
  265. trunc:=float64_to_int64_round_to_zero(f64);
  266. {$endif VER1_0}
  267. end
  268. else
  269. begin
  270. move(d,f32,sizeof(f32));
  271. trunc:=float32_to_int32_round_to_zero(f32);
  272. end;
  273. end;
  274. {$endif}
  275. {$ifndef FPC_SYSTEM_HAS_INT}
  276. {$ifdef SUPPORT_DOUBLE}
  277. { straight Pascal translation of the code for __trunc() in }
  278. { the file sysdeps/libm-ieee754/s_trunc.c of glibc (JM) }
  279. function int(d: double): double;[internconst:in_const_int];
  280. var
  281. i0, j0: longint;
  282. i1: cardinal;
  283. sx: longint;
  284. f64 : float64;
  285. begin
  286. f64:=float64(d);
  287. {$ifdef cpuarm}
  288. { the arm fpu has a strange opinion how a double has to be stored }
  289. i0:=f64.low;
  290. f64.low:=f64.high;
  291. f64.high:=i0;
  292. {$endif cpuarm}
  293. i0 := f64.high;
  294. i1 := cardinal(f64.low);
  295. sx := i0 and $80000000;
  296. j0 := ((i0 shr 20) and $7ff) - $3ff;
  297. if (j0 < 20) then
  298. begin
  299. if (j0 < 0) then
  300. begin
  301. { the magnitude of the number is < 1 so the result is +-0. }
  302. f64.high := sx;
  303. f64.low := 0;
  304. end
  305. else
  306. begin
  307. f64.high := sx or (i0 and not($fffff shr j0));
  308. f64.low := 0;
  309. end
  310. end
  311. else if (j0 > 51) then
  312. begin
  313. if (j0 = $400) then
  314. { d is inf or NaN }
  315. exit(d + d); { don't know why they do this (JM) }
  316. end
  317. else
  318. begin
  319. f64.high := i0;
  320. f64.low := longint(i1 and not(cardinal($ffffffff) shr (j0 - 20)));
  321. end;
  322. {$ifdef cpuarm}
  323. { the arm fpu has a strange opinion how a double has to be stored }
  324. i0:=f64.low;
  325. f64.low:=f64.high;
  326. f64.high:=i0;
  327. {$endif cpuarm}
  328. result:=double(f64);
  329. end;
  330. {$else SUPPORT_DOUBLE}
  331. function int(d : real) : real;[internconst:in_const_int];
  332. begin
  333. { this will be correct since real = single in the case of }
  334. { the motorola version of the compiler... }
  335. int:=real(trunc(d));
  336. end;
  337. {$endif SUPPORT_DOUBLE}
  338. {$endif}
  339. {$ifndef FPC_SYSTEM_HAS_ABS}
  340. {$ifdef SUPPORT_DOUBLE}
  341. function abs(d : Double) : Double;[public,alias:'FPC_ABS_REAL'];
  342. begin
  343. if (d<0.0) then
  344. abs := -d
  345. else
  346. abs := d ;
  347. end;
  348. {$else}
  349. function abs(d : Real) : Real;[public,alias:'FPC_ABS_REAL'];
  350. begin
  351. if (d<0.0) then
  352. abs := -d
  353. else
  354. abs := d ;
  355. end;
  356. {$endif}
  357. {$ifdef hascompilerproc}
  358. function fpc_abs_real(d:Real):Real;compilerproc; external name 'FPC_ABS_REAL';
  359. {$endif hascompilerproc}
  360. {$endif not FPC_SYSTEM_HAS_ABS}
  361. {$ifndef SYSTEM_HAS_FREXP}
  362. function frexp(x:Real; var e:Integer ):Real;
  363. {* frexp() extracts the exponent from x. It returns an integer *}
  364. {* power of two to expnt and the significand between 0.5 and 1 *}
  365. {* to y. Thus x = y * 2**expn. *}
  366. begin
  367. e :=0;
  368. if (abs(x)<0.5) then
  369. While (abs(x)<0.5) do
  370. begin
  371. x := x*2;
  372. Dec(e);
  373. end
  374. else
  375. While (abs(x)>1) do
  376. begin
  377. x := x/2;
  378. Inc(e);
  379. end;
  380. frexp := x;
  381. end;
  382. {$endif not SYSTEM_HAS_FREXP}
  383. {$ifndef SYSTEM_HAS_LDEXP}
  384. function ldexp( x: Real; N: Integer):Real;
  385. {* ldexp() multiplies x by 2**n. *}
  386. var r : Real;
  387. begin
  388. R := 1;
  389. if N>0 then
  390. while N>0 do
  391. begin
  392. R:=R*2;
  393. Dec(N);
  394. end
  395. else
  396. while N<0 do
  397. begin
  398. R:=R/2;
  399. Inc(N);
  400. end;
  401. ldexp := x * R;
  402. end;
  403. {$endif not SYSTEM_HAS_LDEXP}
  404. function polevl(var x:Real; var Coef:TabCoef; N:Integer):Real;
  405. {*****************************************************************}
  406. { Evaluate polynomial }
  407. {*****************************************************************}
  408. { }
  409. { SYNOPSIS: }
  410. { }
  411. { int N; }
  412. { double x, y, coef[N+1], polevl[]; }
  413. { }
  414. { y = polevl( x, coef, N ); }
  415. { }
  416. { DESCRIPTION: }
  417. { }
  418. { Evaluates polynomial of degree N: }
  419. { }
  420. { 2 N }
  421. { y = C + C x + C x +...+ C x }
  422. { 0 1 2 N }
  423. { }
  424. { Coefficients are stored in reverse order: }
  425. { }
  426. { coef[0] = C , ..., coef[N] = C . }
  427. { N 0 }
  428. { }
  429. { The function p1evl() assumes that coef[N] = 1.0 and is }
  430. { omitted from the array. Its calling arguments are }
  431. { otherwise the same as polevl(). }
  432. { }
  433. { SPEED: }
  434. { }
  435. { In the interest of speed, there are no checks for out }
  436. { of bounds arithmetic. This routine is used by most of }
  437. { the functions in the library. Depending on available }
  438. { equipment features, the user may wish to rewrite the }
  439. { program in microcode or assembly language. }
  440. {*****************************************************************}
  441. var ans : Real;
  442. i : Integer;
  443. begin
  444. ans := Coef[0];
  445. for i:=1 to N do
  446. ans := ans * x + Coef[i];
  447. polevl:=ans;
  448. end;
  449. function p1evl(var x:Real; var Coef:TabCoef; N:Integer):Real;
  450. { }
  451. { Evaluate polynomial when coefficient of x is 1.0. }
  452. { Otherwise same as polevl. }
  453. { }
  454. var
  455. ans : Real;
  456. i : Integer;
  457. begin
  458. ans := x + Coef[0];
  459. for i:=1 to N-1 do
  460. ans := ans * x + Coef[i];
  461. p1evl := ans;
  462. end;
  463. {$ifndef FPC_SYSTEM_HAS_SQR}
  464. function sqr(d : Real) : Real;[internconst:in_const_sqr];
  465. begin
  466. sqr := d*d;
  467. end;
  468. {$endif}
  469. {$ifndef FPC_SYSTEM_HAS_PI}
  470. function pi : Real;[internconst:in_const_pi];
  471. begin
  472. pi := 3.1415926535897932385;
  473. end;
  474. {$endif}
  475. {$ifndef FPC_SYSTEM_HAS_SQRT}
  476. function sqrt(d:Real):Real;[internconst:in_const_sqrt]; [public, alias: 'FPC_SQRT_REAL'];
  477. {*****************************************************************}
  478. { Square root }
  479. {*****************************************************************}
  480. { }
  481. { SYNOPSIS: }
  482. { }
  483. { double x, y, sqrt(); }
  484. { }
  485. { y = sqrt( x ); }
  486. { }
  487. { DESCRIPTION: }
  488. { }
  489. { Returns the square root of x. }
  490. { }
  491. { Range reduction involves isolating the power of two of the }
  492. { argument and using a polynomial approximation to obtain }
  493. { a rough value for the square root. Then Heron's iteration }
  494. { is used three times to converge to an accurate value. }
  495. {*****************************************************************}
  496. var e : Integer;
  497. w,z : Real;
  498. begin
  499. if( d <= 0.0 ) then
  500. begin
  501. if( d < 0.0 ) then
  502. HandleError(207);
  503. sqrt := 0.0;
  504. end
  505. else
  506. begin
  507. w := d;
  508. { separate exponent and significand }
  509. z := frexp( d, e );
  510. { approximate square root of number between 0.5 and 1 }
  511. { relative error of approximation = 7.47e-3 }
  512. d := 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
  513. { adjust for odd powers of 2 }
  514. if odd(e) then
  515. d := d*SQRT2;
  516. { re-insert exponent }
  517. d := ldexp( d, (e div 2) );
  518. { Newton iterations: }
  519. d := 0.5*(d + w/d);
  520. d := 0.5*(d + w/d);
  521. d := 0.5*(d + w/d);
  522. d := 0.5*(d + w/d);
  523. d := 0.5*(d + w/d);
  524. d := 0.5*(d + w/d);
  525. sqrt := d;
  526. end;
  527. end;
  528. {$ifdef hascompilerproc}
  529. function fpc_sqrt_real(d:Real):Real;compilerproc; external name 'FPC_SQRT_REAL';
  530. {$endif hascompilerproc}
  531. {$endif}
  532. {$ifndef FPC_SYSTEM_HAS_EXP}
  533. function Exp(d:Real):Real;[internconst:in_const_exp];
  534. {*****************************************************************}
  535. { Exponential Function }
  536. {*****************************************************************}
  537. { }
  538. { SYNOPSIS: }
  539. { }
  540. { double x, y, exp(); }
  541. { }
  542. { y = exp( x ); }
  543. { }
  544. { DESCRIPTION: }
  545. { }
  546. { Returns e (2.71828...) raised to the x power. }
  547. { }
  548. { Range reduction is accomplished by separating the argument }
  549. { into an integer k and fraction f such that }
  550. { }
  551. { x k f }
  552. { e = 2 e. }
  553. { }
  554. { A Pade' form of degree 2/3 is used to approximate exp(f)- 1 }
  555. { in the basic range [-0.5 ln 2, 0.5 ln 2]. }
  556. {*****************************************************************}
  557. const P : TabCoef = (
  558. 1.26183092834458542160E-4,
  559. 3.02996887658430129200E-2,
  560. 1.00000000000000000000E0, 0, 0, 0, 0);
  561. Q : TabCoef = (
  562. 3.00227947279887615146E-6,
  563. 2.52453653553222894311E-3,
  564. 2.27266044198352679519E-1,
  565. 2.00000000000000000005E0, 0 ,0 ,0);
  566. C1 = 6.9335937500000000000E-1;
  567. C2 = 2.1219444005469058277E-4;
  568. var n : Integer;
  569. px, qx, xx : Real;
  570. begin
  571. if( d > MAXLOG) then
  572. HandleError(205)
  573. else
  574. if( d < MINLOG ) then
  575. begin
  576. HandleError(205);
  577. end
  578. else
  579. begin
  580. { Express e**x = e**g 2**n }
  581. { = e**g e**( n loge(2) ) }
  582. { = e**( g + n loge(2) ) }
  583. px := d * LOG2E;
  584. qx := Trunc( px + 0.5 ); { Trunc() truncates toward -infinity. }
  585. n := Trunc(qx);
  586. d := d - qx * C1;
  587. d := d + qx * C2;
  588. { rational approximation for exponential }
  589. { of the fractional part: }
  590. { e**x - 1 = 2x P(x**2)/( Q(x**2) - P(x**2) ) }
  591. xx := d * d;
  592. px := d * polevl( xx, P, 2 );
  593. d := px/( polevl( xx, Q, 3 ) - px );
  594. d := ldexp( d, 1 );
  595. d := d + 1.0;
  596. d := ldexp( d, n );
  597. Exp := d;
  598. end;
  599. end;
  600. {$endif}
  601. {$ifndef FPC_SYSTEM_HAS_ROUND}
  602. {$ifdef hascompilerproc}
  603. function round(d : Real) : int64;[internconst:in_const_round, external name 'FPC_ROUND'];
  604. function fpc_round(d : Real) : int64;[public, alias:'FPC_ROUND'];{$ifdef hascompilerproc}compilerproc;{$endif hascompilerproc}
  605. {$else}
  606. function round(d : Real) : int64;[internconst:in_const_round];
  607. {$endif hascompilerproc}
  608. var
  609. fr: Real;
  610. tr: Int64;
  611. Begin
  612. fr := abs(Frac(d));
  613. tr := Trunc(d);
  614. if fr > 0.5 then
  615. if d >= 0 then
  616. result:=tr+1
  617. else
  618. result:=tr-1
  619. else
  620. if fr < 0.5 then
  621. result:=tr
  622. else { fr = 0.5 }
  623. { check sign to decide ... }
  624. { as in Turbo Pascal... }
  625. if d >= 0.0 then
  626. result:=tr+1
  627. else
  628. result:=tr;
  629. end;
  630. {$endif}
  631. {$ifdef FPC_CURRENCY_IS_INT64}
  632. function trunc(c : currency) : int64;
  633. type
  634. tmyrec = record
  635. i: int64;
  636. end;
  637. begin
  638. result := int64(tmyrec(c)) div 10000
  639. end;
  640. function trunc(c : comp) : int64;
  641. begin
  642. result := c
  643. end;
  644. function round(c : currency) : int64;
  645. type
  646. tmyrec = record
  647. i: int64;
  648. end;
  649. var
  650. rem, absrem: longint;
  651. begin
  652. { (int64(tmyrec(c))(+/-)5000) div 10000 can overflow }
  653. result := int64(tmyrec(c)) div 10000;
  654. rem := int64(tmyrec(c)) - result * 10000;
  655. absrem := abs(rem);
  656. if (absrem > 5000) or
  657. ((absrem = 5000) and
  658. (rem > 0)) then
  659. if (rem > 0) then
  660. inc(result)
  661. else
  662. dec(result);
  663. end;
  664. function round(c : comp) : int64;
  665. begin
  666. result := c
  667. end;
  668. {$endif FPC_CURRENCY_IS_INT64}
  669. {$ifndef FPC_SYSTEM_HAS_LN}
  670. function Ln(d:Real):Real;[internconst:in_const_ln];
  671. {*****************************************************************}
  672. { Natural Logarithm }
  673. {*****************************************************************}
  674. { }
  675. { SYNOPSIS: }
  676. { }
  677. { double x, y, log(); }
  678. { }
  679. { y = ln( x ); }
  680. { }
  681. { DESCRIPTION: }
  682. { }
  683. { Returns the base e (2.718...) logarithm of x. }
  684. { }
  685. { The argument is separated into its exponent and fractional }
  686. { parts. If the exponent is between -1 and +1, the logarithm }
  687. { of the fraction is approximated by }
  688. { }
  689. { log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). }
  690. { }
  691. { Otherwise, setting z = 2(x-1)/x+1), }
  692. { }
  693. { log(x) = z + z**3 P(z)/Q(z). }
  694. { }
  695. {*****************************************************************}
  696. const P : TabCoef = (
  697. { Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  698. 1/sqrt(2) <= x < sqrt(2) }
  699. 4.58482948458143443514E-5,
  700. 4.98531067254050724270E-1,
  701. 6.56312093769992875930E0,
  702. 2.97877425097986925891E1,
  703. 6.06127134467767258030E1,
  704. 5.67349287391754285487E1,
  705. 1.98892446572874072159E1);
  706. Q : TabCoef = (
  707. 1.50314182634250003249E1,
  708. 8.27410449222435217021E1,
  709. 2.20664384982121929218E2,
  710. 3.07254189979530058263E2,
  711. 2.14955586696422947765E2,
  712. 5.96677339718622216300E1, 0);
  713. { Coefficients for log(x) = z + z**3 P(z)/Q(z),
  714. where z = 2(x-1)/(x+1)
  715. 1/sqrt(2) <= x < sqrt(2) }
  716. R : TabCoef = (
  717. -7.89580278884799154124E-1,
  718. 1.63866645699558079767E1,
  719. -6.41409952958715622951E1, 0, 0, 0, 0);
  720. S : TabCoef = (
  721. -3.56722798256324312549E1,
  722. 3.12093766372244180303E2,
  723. -7.69691943550460008604E2, 0, 0, 0, 0);
  724. var e : Integer;
  725. z, y : Real;
  726. Label Ldone;
  727. begin
  728. if( d <= 0.0 ) then
  729. HandleError(207);
  730. d := frexp( d, e );
  731. { logarithm using log(x) = z + z**3 P(z)/Q(z),
  732. where z = 2(x-1)/x+1) }
  733. if( (e > 2) or (e < -2) ) then
  734. begin
  735. if( d < SQRTH ) then
  736. begin
  737. { 2( 2x-1 )/( 2x+1 ) }
  738. Dec(e, 1);
  739. z := d - 0.5;
  740. y := 0.5 * z + 0.5;
  741. end
  742. else
  743. begin
  744. { 2 (x-1)/(x+1) }
  745. z := d - 0.5;
  746. z := z - 0.5;
  747. y := 0.5 * d + 0.5;
  748. end;
  749. d := z / y;
  750. { /* rational form */ }
  751. z := d*d;
  752. z := d + d * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
  753. goto ldone;
  754. end;
  755. { logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) }
  756. if( d < SQRTH ) then
  757. begin
  758. Dec(e, 1);
  759. d := ldexp( d, 1 ) - 1.0; { 2x - 1 }
  760. end
  761. else
  762. d := d - 1.0;
  763. { rational form }
  764. z := d*d;
  765. y := d * ( z * polevl( d, P, 6 ) / p1evl( d, Q, 6 ) );
  766. y := y - ldexp( z, -1 ); { y - 0.5 * z }
  767. z := d + y;
  768. ldone:
  769. { recombine with exponent term }
  770. if( e <> 0 ) then
  771. begin
  772. y := e;
  773. z := z - y * 2.121944400546905827679e-4;
  774. z := z + y * 0.693359375;
  775. end;
  776. Ln:= z;
  777. end;
  778. {$endif}
  779. {$ifndef FPC_SYSTEM_HAS_SIN}
  780. function Sin(d:Real):Real;[internconst:in_const_sin];
  781. {*****************************************************************}
  782. { Circular Sine }
  783. {*****************************************************************}
  784. { }
  785. { SYNOPSIS: }
  786. { }
  787. { double x, y, sin(); }
  788. { }
  789. { y = sin( x ); }
  790. { }
  791. { DESCRIPTION: }
  792. { }
  793. { Range reduction is into intervals of pi/4. The reduction }
  794. { error is nearly eliminated by contriving an extended }
  795. { precision modular arithmetic. }
  796. { }
  797. { Two polynomial approximating functions are employed. }
  798. { Between 0 and pi/4 the sine is approximated by }
  799. { x + x**3 P(x**2). }
  800. { Between pi/4 and pi/2 the cosine is represented as }
  801. { 1 - x**2 Q(x**2). }
  802. {*****************************************************************}
  803. var y, z, zz : Real;
  804. j, sign : Integer;
  805. begin
  806. { make argument positive but save the sign }
  807. sign := 1;
  808. if( d < 0 ) then
  809. begin
  810. d := -d;
  811. sign := -1;
  812. end;
  813. { above this value, approximate towards 0 }
  814. if( d > lossth ) then
  815. begin
  816. sin := 0.0;
  817. exit;
  818. end;
  819. y := Trunc( d/PIO4 ); { integer part of x/PIO4 }
  820. { strip high bits of integer part to prevent integer overflow }
  821. z := ldexp( y, -4 );
  822. z := Trunc(z); { integer part of y/8 }
  823. z := y - ldexp( z, 4 ); { y - 16 * (y/16) }
  824. j := Trunc(z); { convert to integer for tests on the phase angle }
  825. { map zeros to origin }
  826. { typecast is to avoid "can't determine which overloaded function }
  827. { to call" }
  828. if odd( longint(j) ) then
  829. begin
  830. inc(j);
  831. y := y + 1.0;
  832. end;
  833. j := j and 7; { octant modulo 360 degrees }
  834. { reflect in x axis }
  835. if( j > 3) then
  836. begin
  837. sign := -sign;
  838. dec(j, 4);
  839. end;
  840. { Extended precision modular arithmetic }
  841. z := ((d - y * DP1) - y * DP2) - y * DP3;
  842. zz := z * z;
  843. if( (j=1) or (j=2) ) then
  844. y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 )
  845. else
  846. { y = z + z * (zz * polevl( zz, sincof, 5 )); }
  847. y := z + z * z * z * polevl( zz, sincof, 5 );
  848. if(sign < 0) then
  849. y := -y;
  850. sin := y;
  851. end;
  852. {$endif}
  853. {$ifndef FPC_SYSTEM_HAS_COS}
  854. function Cos(d:Real):Real;[internconst:in_const_cos];
  855. {*****************************************************************}
  856. { Circular cosine }
  857. {*****************************************************************}
  858. { }
  859. { Circular cosine }
  860. { }
  861. { SYNOPSIS: }
  862. { }
  863. { double x, y, cos(); }
  864. { }
  865. { y = cos( x ); }
  866. { }
  867. { DESCRIPTION: }
  868. { }
  869. { Range reduction is into intervals of pi/4. The reduction }
  870. { error is nearly eliminated by contriving an extended }
  871. { precision modular arithmetic. }
  872. { }
  873. { Two polynomial approximating functions are employed. }
  874. { Between 0 and pi/4 the cosine is approximated by }
  875. { 1 - x**2 Q(x**2). }
  876. { Between pi/4 and pi/2 the sine is represented as }
  877. { x + x**3 P(x**2). }
  878. {*****************************************************************}
  879. var y, z, zz : Real;
  880. j, sign : Integer;
  881. i : LongInt;
  882. begin
  883. { make argument positive }
  884. sign := 1;
  885. if( d < 0 ) then
  886. d := -d;
  887. { above this value, round towards zero }
  888. if( d > lossth ) then
  889. begin
  890. cos := 0.0;
  891. exit;
  892. end;
  893. y := Trunc( d/PIO4 );
  894. z := ldexp( y, -4 );
  895. z := Trunc(z); { integer part of y/8 }
  896. z := y - ldexp( z, 4 ); { y - 16 * (y/16) }
  897. { integer and fractional part modulo one octant }
  898. i := Trunc(z);
  899. if odd( i ) then { map zeros to origin }
  900. begin
  901. inc(i);
  902. y := y + 1.0;
  903. end;
  904. j := i and 07;
  905. if( j > 3) then
  906. begin
  907. dec(j,4);
  908. sign := -sign;
  909. end;
  910. if( j > 1 ) then
  911. sign := -sign;
  912. { Extended precision modular arithmetic }
  913. z := ((d - y * DP1) - y * DP2) - y * DP3;
  914. zz := z * z;
  915. if( (j=1) or (j=2) ) then
  916. { y = z + z * (zz * polevl( zz, sincof, 5 )); }
  917. y := z + z * z * z * polevl( zz, sincof, 5 )
  918. else
  919. y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
  920. if(sign < 0) then
  921. y := -y;
  922. cos := y ;
  923. end;
  924. {$endif}
  925. {$ifndef FPC_SYSTEM_HAS_ARCTAN}
  926. function ArcTan(d:Real):Real;[internconst:in_const_arctan];
  927. {*****************************************************************}
  928. { Inverse circular tangent (arctangent) }
  929. {*****************************************************************}
  930. { }
  931. { SYNOPSIS: }
  932. { }
  933. { double x, y, atan(); }
  934. { }
  935. { y = atan( x ); }
  936. { }
  937. { DESCRIPTION: }
  938. { }
  939. { Returns radian angle between -pi/2 and +pi/2 whose tangent }
  940. { is x. }
  941. { }
  942. { Range reduction is from four intervals into the interval }
  943. { from zero to tan( pi/8 ). The approximant uses a rational }
  944. { function of degree 3/4 of the form x + x**3 P(x)/Q(x). }
  945. {*****************************************************************}
  946. const P : TabCoef = (
  947. -8.40980878064499716001E-1,
  948. -8.83860837023772394279E0,
  949. -2.18476213081316705724E1,
  950. -1.48307050340438946993E1, 0, 0, 0);
  951. Q : TabCoef = (
  952. 1.54974124675307267552E1,
  953. 6.27906555762653017263E1,
  954. 9.22381329856214406485E1,
  955. 4.44921151021319438465E1, 0, 0, 0);
  956. { tan( 3*pi/8 ) }
  957. T3P8 = 2.41421356237309504880;
  958. { tan( pi/8 ) }
  959. TP8 = 0.41421356237309504880;
  960. var y,z : Real;
  961. Sign : Integer;
  962. begin
  963. { make argument positive and save the sign }
  964. sign := 1;
  965. if( d < 0.0 ) then
  966. begin
  967. sign := -1;
  968. d := -d;
  969. end;
  970. { range reduction }
  971. if( d > T3P8 ) then
  972. begin
  973. y := PIO2;
  974. d := -( 1.0/d );
  975. end
  976. else if( d > TP8 ) then
  977. begin
  978. y := PIO4;
  979. d := (d-1.0)/(d+1.0);
  980. end
  981. else
  982. y := 0.0;
  983. { rational form in x**2 }
  984. z := d * d;
  985. y := y + ( polevl( z, P, 3 ) / p1evl( z, Q, 4 ) ) * z * d + d;
  986. if( sign < 0 ) then
  987. y := -y;
  988. Arctan := y;
  989. end;
  990. {$endif}
  991. {$ifndef FPC_SYSTEM_HAS_FRAC}
  992. function frac(d : Real) : Real;[internconst:in_const_frac];
  993. begin
  994. frac := d - Int(d);
  995. end;
  996. {$endif}
  997. {$ifndef FPC_SYSTEM_HAS_POWER}
  998. function power(bas,expo : real) : real;
  999. begin
  1000. if bas=0.0 then
  1001. begin
  1002. if expo<>0.0 then
  1003. power:=0.0
  1004. else
  1005. HandleError(207);
  1006. end
  1007. else if expo=0.0 then
  1008. power:=1
  1009. else
  1010. { bas < 0 is not allowed }
  1011. if bas<0.0 then
  1012. handleerror(207)
  1013. else
  1014. power:=exp(ln(bas)*expo);
  1015. end;
  1016. {$endif}
  1017. {$ifndef FPC_SYSTEM_HAS_POWER_INT64}
  1018. function power(bas,expo : int64) : int64;
  1019. begin
  1020. if bas=0 then
  1021. begin
  1022. if expo<>0 then
  1023. power:=0
  1024. else
  1025. HandleError(207);
  1026. end
  1027. else if expo=0 then
  1028. power:=1
  1029. else
  1030. begin
  1031. if bas<0 then
  1032. begin
  1033. if odd(expo) then
  1034. power:=-round(exp(ln(-bas)*expo))
  1035. else
  1036. power:=round(exp(ln(-bas)*expo));
  1037. end
  1038. else
  1039. power:=round(exp(ln(bas)*expo));
  1040. end;
  1041. end;
  1042. {$endif}
  1043. {$ifdef FPC_INCLUDE_SOFTWARE_INT64_TO_DOUBLE}
  1044. {$ifndef FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
  1045. function fpc_qword_to_double(q : qword): double; compilerproc;
  1046. begin
  1047. result:=dword(q and $ffffffff)+dword(q shr 32)*4294967296.0;
  1048. end;
  1049. {$endif FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  1050. {$ifndef FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  1051. function fpc_int64_to_double(i : int64): double; compilerproc;
  1052. begin
  1053. if i<0 then
  1054. result:=-double(qword(-i))
  1055. else
  1056. result:=qword(i);
  1057. end;
  1058. {$endif FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  1059. {$endif FPC_INCLUDE_SOFTWARE_INT64_TO_DOUBLE}
  1060. {$ifdef SUPPORT_DOUBLE}
  1061. {****************************************************************************
  1062. Helper routines to support old TP styled reals
  1063. ****************************************************************************}
  1064. {$ifndef FPC_SYSTEM_HAS_REAL2DOUBLE}
  1065. function real2double(r : real48) : double;
  1066. var
  1067. res : array[0..7] of byte;
  1068. exponent : word;
  1069. begin
  1070. { copy mantissa }
  1071. res[0]:=0;
  1072. res[1]:=r[1] shl 5;
  1073. res[2]:=(r[1] shr 3) or (r[2] shl 5);
  1074. res[3]:=(r[2] shr 3) or (r[3] shl 5);
  1075. res[4]:=(r[3] shr 3) or (r[4] shl 5);
  1076. res[5]:=(r[4] shr 3) or (r[5] and $7f) shl 5;
  1077. res[6]:=(r[5] and $7f) shr 3;
  1078. { copy exponent }
  1079. { correct exponent: }
  1080. exponent:=(word(r[0])+(1023-129));
  1081. res[6]:=res[6] or ((exponent and $f) shl 4);
  1082. res[7]:=exponent shr 4;
  1083. { set sign }
  1084. res[7]:=res[7] or (r[5] and $80);
  1085. real2double:=double(res);
  1086. end;
  1087. {$endif FPC_SYSTEM_HAS_REAL2DOUBLE}
  1088. {$endif SUPPORT_DOUBLE}
  1089. {
  1090. $Log$
  1091. Revision 1.27 2004-10-09 21:00:46 jonas
  1092. + cgenmath with libc math functions. Faster than the routines in genmath
  1093. and also have full double support (exp() only has support for values in
  1094. the single range in genmath, for example). Used in FPC_USE_LIBC is
  1095. defined
  1096. * several fixes to allow compilation with -dHASINLINE, but internalerrors
  1097. because of missing support for inlining assembler code
  1098. Revision 1.26 2004/10/03 14:09:39 florian
  1099. * fixed trunc for abs(value) < 1
  1100. Revision 1.25 2004/10/03 14:00:21 florian
  1101. + made generic trunc 64 bit aware
  1102. Revision 1.24 2004/05/31 20:25:04 peter
  1103. * removed warnings
  1104. Revision 1.23 2004/03/13 18:33:52 florian
  1105. * fixed some arm related real stuff
  1106. Revision 1.22 2004/03/11 22:39:53 florian
  1107. * arm startup code fixed
  1108. * made some generic math code more readable
  1109. Revision 1.21 2004/02/04 14:15:57 florian
  1110. * fixed generic system.int(...)
  1111. Revision 1.20 2004/01/24 18:15:58 florian
  1112. * fixed small bugs
  1113. * fixed some arm issues
  1114. Revision 1.18 2004/01/06 21:34:07 peter
  1115. * abs(double) added
  1116. * abs() alias
  1117. Revision 1.17 2004/01/02 17:19:04 jonas
  1118. * if currency = int64, FPC_CURRENCY_IS_INT64 is defined
  1119. + round and trunc for currency and comp if FPC_CURRENCY_IS_INT64 is
  1120. defined
  1121. * if currency = orddef, prefer currency -> int64/qword conversion over
  1122. currency -> float conversions
  1123. * optimized currency/currency if currency = orddef
  1124. * TODO: write FPC_DIV_CURRENCY and FPC_MUL_CURRENCY routines to prevent
  1125. precision loss if currency=int64 and bestreal = double
  1126. Revision 1.16 2003/12/08 19:44:11 jonas
  1127. * use HandleError instead of RunError so exception catching works
  1128. Revision 1.15 2003/09/03 14:09:37 florian
  1129. * arm fixes to the common rtl code
  1130. * some generic math code fixed
  1131. * ...
  1132. Revision 1.14 2003/05/24 13:39:32 jonas
  1133. * fsqrt is an optional instruction in the ppc architecture and isn't
  1134. implemented by any current ppc afaik, so use the generic sqrt routine
  1135. instead (adapted so it works with compilerproc)
  1136. Revision 1.13 2003/05/23 22:58:31 jonas
  1137. * added longint typecase to odd(smallint_var) call to avoid overload
  1138. problem
  1139. Revision 1.12 2003/05/02 15:12:19 jonas
  1140. - removed empty ppc-specific frac()
  1141. + added correct generic frac() implementation for doubles (translated
  1142. from glibc code)
  1143. Revision 1.11 2003/04/23 21:28:21 peter
  1144. * fpc_round added, needed for int64 currency
  1145. Revision 1.10 2003/01/15 00:45:17 peter
  1146. * use generic int64 power
  1147. Revision 1.9 2002/10/12 20:28:49 carl
  1148. * round returns int64
  1149. Revision 1.8 2002/10/07 15:15:02 florian
  1150. * fixed wrong commit
  1151. Revision 1.7 2002/10/07 15:10:45 florian
  1152. + variant wrappers for cmp operators added
  1153. Revision 1.6 2002/09/07 15:07:45 peter
  1154. * old logs removed and tabs fixed
  1155. Revision 1.5 2002/07/28 21:39:29 florian
  1156. * made abs a compiler proc if it is generic
  1157. Revision 1.4 2002/07/28 20:43:48 florian
  1158. * several fixes for linux/powerpc
  1159. * several fixes to MT
  1160. }