math.inc 13 KB

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