math.inc 12 KB

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