math.inc 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468
  1. {
  2. $Id$
  3. This file is part of the Free Pascal run time library.
  4. Copyright (c) 2000 by Jonas Maebe and other members of the
  5. Free Pascal development team
  6. Implementation of mathematical Routines (only for real)
  7. See the file COPYING.FPC, included in this distribution,
  8. for details about the copyright.
  9. This program is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  12. **********************************************************************}
  13. const
  14. longint_to_real_helper: int64 = $4330000080000000;
  15. cardinal_to_real_helper: int64 = $4330000000000000;
  16. int_to_real_factor: double = double(high(cardinal))+1.0;
  17. {****************************************************************************
  18. EXTENDED data type routines
  19. ****************************************************************************}
  20. {$ifdef INTERNCONSTINTF}
  21. {$define FPC_SYSTEM_HAS_PI}
  22. function fpc_pi_real : valreal;compilerproc;
  23. begin
  24. { Function is handled internal in the compiler }
  25. runerror(207);
  26. result:=0;
  27. end;
  28. {$define FPC_SYSTEM_HAS_ABS}
  29. function fpc_abs_real(d : valreal) : valreal;compilerproc;
  30. begin
  31. { Function is handled internal in the compiler }
  32. runerror(207);
  33. result:=0;
  34. end;
  35. {$define FPC_SYSTEM_HAS_SQR}
  36. function fpc_sqr_real(d : valreal) : valreal;compilerproc;
  37. begin
  38. { Function is handled internal in the compiler }
  39. runerror(207);
  40. result:=0;
  41. end;
  42. {$else}
  43. {$define FPC_SYSTEM_HAS_PI}
  44. function pi : double;[internproc:fpc_in_pi];
  45. {$define FPC_SYSTEM_HAS_ABS}
  46. function abs(d : extended) : extended;[internproc:fpc_in_abs_extended];
  47. {$define FPC_SYSTEM_HAS_SQR}
  48. function sqr(d : extended) : extended;[internproc:fpc_in_sqr_extended];
  49. {$endif ndef INTERNCONSTINTF}
  50. const
  51. factor: double = double(int64(1) shl 32);
  52. factor2: double = double(int64(1) shl 31);
  53. {$ifndef FPC_SYSTEM_HAS_TRUNC}
  54. {$define FPC_SYSTEM_HAS_TRUNC}
  55. {$ifdef INTERNCONSTINTF}
  56. function fpc_trunc_real(d : valreal) : int64;assembler;compilerproc;
  57. {$else}
  58. function trunc(d : extended) : int64;assembler;[internconst:fpc_in_const_trunc];
  59. {$endif}
  60. { input: d in fr1 }
  61. { output: result in r3 }
  62. assembler;
  63. var
  64. temp: packed record
  65. case byte of
  66. 0: (l1,l2: longint);
  67. 1: (d: double);
  68. end;
  69. asm
  70. // store d in temp
  71. stfd f1,temp
  72. // extract sign bit (record in cr0)
  73. lwz r3,temp
  74. rlwinm. r3,r3,1,31,31
  75. // make d positive
  76. fabs f1,f1
  77. // load 2^32 in f2
  78. {$ifndef macos}
  79. lis r4,factor@ha
  80. lfd f2,factor@l(r4)
  81. {$else}
  82. lwz r4,factor(r2)
  83. lfd f2,0(r4)
  84. {$endif}
  85. // check if value is < 0
  86. // f3 := d / 2^32;
  87. fdiv f3,f1,f2
  88. // round
  89. fctiwz f4,f3
  90. // store
  91. stfd f4,temp
  92. // and load into r4
  93. lwz r3,temp+4
  94. // convert back to float
  95. lis r0,0x4330
  96. stw r0,temp
  97. xoris r0,r3,0x8000
  98. stw r0,temp+4
  99. {$ifndef macos}
  100. lis r4,longint_to_real_helper@ha
  101. lfd f0,longint_to_real_helper@l(r4)
  102. {$else}
  103. lwz r4,longint_to_real_helper(r2)
  104. lfd f0,0(r4)
  105. {$endif}
  106. lfd f3,temp
  107. fsub f3,f3,f0
  108. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  109. fnmsub f4,f3,f2,f1
  110. // now, convert to unsigned 32 bit
  111. // load 2^31 in f2
  112. {$ifndef macos}
  113. lis r4,factor2@ha
  114. lfd f2,factor2@l(r4)
  115. {$else}
  116. lwz r4,factor2(r2)
  117. lfd f2,0(r4)
  118. {$endif}
  119. // subtract 2^31
  120. fsub f3,f4,f2
  121. // was the value > 2^31?
  122. fcmpu cr1,f4,f2
  123. // use diff if >= 2^31
  124. fsel f4,f3,f3,f4
  125. // next part same as conversion to signed integer word
  126. fctiwz f4,f4
  127. stfd f4,temp
  128. lwz r4,temp+4
  129. // add 2^31 if value was >=2^31
  130. blt cr1, .LTruncNoAdd
  131. xoris r4,r4,0x8000
  132. .LTruncNoAdd:
  133. // negate value if it was negative to start with
  134. beq cr0,.LTruncPositive
  135. subfic r4,r4,0
  136. subfze r3,r3
  137. .LTruncPositive:
  138. end;
  139. {$endif not FPC_SYSTEM_HAS_TRUNC}
  140. (*
  141. {$ifndef FPC_SYSTEM_HAS_ROUND}
  142. {$define FPC_SYSTEM_HAS_ROUND}
  143. {$ifdef hascompilerproc}
  144. function round(d : extended) : int64;{$ifndef INTERNCONSTINTF}[internconst:fpc_in_const_round, external name 'FPC_ROUND'];{$endif}
  145. function fpc_round(d : extended) : int64;assembler;[public, alias:'FPC_ROUND'];{$ifdef hascompilerproc}compilerproc;{$endif hascompilerproc}
  146. {$else}
  147. function round(d : extended) : int64;assembler;{$ifndef INTERNCONSTINTF}[internconst:fpc_in_const_round];{$endif}
  148. {$endif hascompilerproc}
  149. { exactly the same as trunc, except that one fctiwz has become fctiw }
  150. { input: d in fr1 }
  151. { output: result in r3 }
  152. assembler;
  153. var
  154. temp: packed record
  155. case byte of
  156. 0: (l1,l2: longint);
  157. 1: (d: double);
  158. end;
  159. asm
  160. // store d in temp
  161. stfd f1, temp
  162. // extract sign bit (record in cr0)
  163. lwz r4,temp
  164. rlwinm. r4,r4,1,31,31
  165. // make d positive
  166. fabs f1,f1
  167. // load 2^32 in f2
  168. {$ifndef macos}
  169. lis r4,factor@ha
  170. lfd f2,factor@l(r4)
  171. {$else}
  172. lwz r4,factor(r2)
  173. lfd f2,0(r4)
  174. {$endif}
  175. // check if value is < 0
  176. // f3 := d / 2^32;
  177. fdiv f3,f1,f2
  178. // round
  179. fctiwz f4,f3
  180. // store
  181. stfd f4,temp
  182. // and load into r4
  183. lwz r3,temp+4
  184. // convert back to float
  185. lis r0,0x4330
  186. stw r0,temp
  187. xoris r0,r3,0x8000
  188. stw r0,temp+4
  189. {$ifndef macos}
  190. lis r4,longint_to_real_helper@ha
  191. lfd f0,longint_to_real_helper@l(r4)
  192. {$else}
  193. lwz r4,longint_to_real_helper(r2)
  194. lfd f0,0(r4)
  195. {$endif}
  196. lfd f3,temp
  197. fsub f3,f3,f0
  198. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  199. fnmsub f4,f3,f2,f1
  200. // now, convert to unsigned 32 bit
  201. // load 2^31 in f2
  202. {$ifndef macos}
  203. lis r4,factor2@ha
  204. lfd f2,factor2@l(r4)
  205. {$else}
  206. lwz r4,factor2(r2)
  207. lfd f2,0(r4)
  208. {$endif}
  209. // subtract 2^31
  210. fsub f3,f4,f2
  211. // was the value > 2^31?
  212. fcmpu cr1,f4,f2
  213. // use diff if >= 2^31
  214. fsel f4,f3,f3,f4
  215. // next part same as conversion to signed integer word
  216. fctiw f4,f4
  217. stfd f4,temp
  218. lwz r4,temp+4
  219. // add 2^31 if value was >=2^31
  220. blt cr1, .LRoundNoAdd
  221. xoris r4,r4,0x8000
  222. .LRoundNoAdd:
  223. // negate value if it was negative to start with
  224. beq cr0,.LRoundPositive
  225. subfic r4,r4,0
  226. subfze r3,r3
  227. .LRoundPositive:
  228. end;
  229. {$endif not FPC_SYSTEM_HAS_ROUND}
  230. *)
  231. {****************************************************************************
  232. Int to real helpers
  233. ****************************************************************************}
  234. {$define FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  235. function fpc_int64_to_double(i: int64): double; compilerproc;
  236. assembler;
  237. { input: high(i) in r4, low(i) in r3 }
  238. { output: double(i) in f0 }
  239. var
  240. temp: packed record
  241. case byte of
  242. 0: (l1,l2: cardinal);
  243. 1: (d: double);
  244. end;
  245. asm
  246. lis r0,0x4330
  247. stw r0,temp
  248. xoris r3,r3,0x8000
  249. stw r3,temp+4
  250. {$ifndef macos}
  251. lis r3,longint_to_real_helper@ha
  252. lfd f1,longint_to_real_helper@l(r3)
  253. {$else}
  254. lwz r3,longint_to_real_helper(r2)
  255. lfd f1,0(r3)
  256. {$endif}
  257. lfd f0,temp
  258. stw r4,temp+4
  259. fsub f0,f0,f1
  260. {$ifndef macos}
  261. lis r4,cardinal_to_real_helper@ha
  262. lfd f1,cardinal_to_real_helper@l(r4)
  263. lis r4,int_to_real_factor@ha
  264. lfd f3,temp
  265. lfd f2,int_to_real_factor@l(r4)
  266. {$else}
  267. lwz r4,cardinal_to_real_helper(r2)
  268. lwz r3,int_to_real_factor(r2)
  269. lfd f3,temp
  270. lfd f1,0(r4)
  271. lfd f2,0(r3)
  272. {$endif}
  273. fsub f3,f3,f1
  274. fmadd f1,f0,f2,f3
  275. end;
  276. {$define FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
  277. function fpc_qword_to_double(q: qword): double; compilerproc;
  278. assembler;
  279. { input: high(q) in r4, low(q) in r3 }
  280. { output: double(q) in f0 }
  281. var
  282. temp: packed record
  283. case byte of
  284. 0: (l1,l2: cardinal);
  285. 1: (d: double);
  286. end;
  287. asm
  288. lis r0,0x4330
  289. stw r0,temp
  290. stw r3,temp+4
  291. lfd f0,temp
  292. {$ifndef macos}
  293. lis r3,cardinal_to_real_helper@ha
  294. lfd f1,cardinal_to_real_helper@l(r3)
  295. {$else}
  296. lwz r3,longint_to_real_helper(r2)
  297. lfd f1,0(r3)
  298. {$endif}
  299. stw r4,temp+4
  300. fsub f0,f0,f1
  301. lfd f3,temp
  302. {$ifndef macos}
  303. lis r4,int_to_real_factor@ha
  304. lfd f2,int_to_real_factor@l(r4)
  305. {$else}
  306. lwz r4,int_to_real_factor(r2)
  307. lfd f2,0(r4)
  308. {$endif}
  309. fsub f3,f3,f1
  310. fmadd f1,f0,f2,f3
  311. end;
  312. {
  313. $Log$
  314. Revision 1.37 2004-11-21 15:35:23 peter
  315. * float routines all use internproc and compilerproc helpers
  316. Revision 1.36 2004/11/20 15:49:21 jonas
  317. * some compilation fixes for powerpc after all the internconst and
  318. internproc changes, still crashes with internalerror(88) for ppc1
  319. on real2str.inc(193,39)
  320. Revision 1.35 2004/11/07 19:09:02 jonas
  321. * use generic round, the ppc-specific one is buggy
  322. Revision 1.34 2004/10/09 21:00:46 jonas
  323. + cgenmath with libc math functions. Faster than the routines in genmath
  324. and also have full double support (exp() only has support for values in
  325. the single range in genmath, for example). Used in FPC_USE_LIBC is
  326. defined
  327. * several fixes to allow compilation with -dHASINLINE, but internalerrors
  328. because of missing support for inlining assembler code
  329. Revision 1.33 2004/02/09 20:21:06 olle
  330. * fixed global variable access in asm
  331. Revision 1.32 2003/12/07 19:55:37 jonas
  332. - reverted previous patch, solved with the new assembler reader
  333. (which didn't understand the new syntax)
  334. Revision 1.30 2003/11/15 19:01:27 florian
  335. * fixed rtl to work with the integrated fpc ppc assembler reader
  336. Revision 1.29 2003/09/04 16:07:31 florian
  337. * fixed qword_to_double conversion on powerpc
  338. Revision 1.28 2003/09/03 14:09:37 florian
  339. * arm fixes to the common rtl code
  340. * some generic math code fixed
  341. * ...
  342. Revision 1.27 2003/08/08 22:02:05 olle
  343. * small bugfix macos
  344. Revision 1.26 2003/06/14 12:41:08 jonas
  345. * fixed compilation problems (removed unnecessary modified registers
  346. lists from procedures)
  347. Revision 1.25 2003/05/31 20:22:06 jonas
  348. * fixed 64 bit results of trunc and round
  349. Revision 1.24 2003/05/30 23:56:41 florian
  350. * fixed parameter passing for int64
  351. Revision 1.23 2003/05/24 13:39:32 jonas
  352. * fsqrt is an optional instruction in the ppc architecture and isn't
  353. implemented by any current ppc afaik, so use the generic sqrt routine
  354. instead (adapted so it works with compilerproc)
  355. Revision 1.22 2003/05/16 16:04:33 jonas
  356. * fixed round() (almost the same as trunc)
  357. Revision 1.21 2003/05/11 18:09:45 jonas
  358. * fixed qword and int64 to double conversion
  359. Revision 1.20 2003/05/02 15:12:19 jonas
  360. - removed empty ppc-specific frac()
  361. + added correct generic frac() implementation for doubles (translated
  362. from glibc code)
  363. Revision 1.19 2003/04/26 20:36:24 jonas
  364. * trunc now also supports int64 (no NaN's etc though)
  365. Revision 1.18 2003/04/26 17:20:16 florian
  366. * fixed trunc, now it's working at least for longint range
  367. Revision 1.17 2003/04/23 21:28:21 peter
  368. * fpc_round added, needed for int64 currency
  369. Revision 1.16 2003/01/16 11:29:11 olle
  370. * changed access of globals to be indirect via TOC
  371. Revision 1.15 2003/01/15 01:09:04 florian
  372. * changed power(...) prototype to int64
  373. Revision 1.14 2002/11/28 11:04:16 olle
  374. * macos: refs to globals in asm adapted to macos
  375. Revision 1.13 2002/10/21 18:08:28 jonas
  376. * round has int64 instead of longint result
  377. Revision 1.12 2002/09/08 13:00:21 jonas
  378. * made pi an internproc instead of internconst
  379. Revision 1.11 2002/09/07 16:01:26 peter
  380. * old logs removed and tabs fixed
  381. Revision 1.10 2002/08/18 22:11:10 florian
  382. * fixed remaining assembler errors
  383. Revision 1.9 2002/08/18 21:37:48 florian
  384. * several errors in inline assembler fixed
  385. Revision 1.8 2002/08/10 17:14:36 jonas
  386. * various fixes, mostly changing the names of the modifies registers to
  387. upper case since that seems to be required by the compiler
  388. Revision 1.7 2002/07/31 16:58:12 jonas
  389. * fixed conversion from int64/qword to double errors
  390. Revision 1.6 2002/07/29 21:28:17 florian
  391. * several fixes to get further with linux/ppc system unit compilation
  392. Revision 1.5 2002/07/28 21:39:29 florian
  393. * made abs a compiler proc if it is generic
  394. Revision 1.4 2002/07/28 20:43:49 florian
  395. * several fixes for linux/powerpc
  396. * several fixes to MT
  397. }