math.inc 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411
  1. {
  2. This file is part of the Free Pascal run time library.
  3. Copyright (c) 2000 by Jonas Maebe and other members of the
  4. Free Pascal development team
  5. Implementation of mathematical Routines (only for 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. const
  13. longint_to_real_helper: int64 = $4330000080000000;
  14. cardinal_to_real_helper: int64 = $4330000000000000;
  15. int_to_real_factor: double = double(high(cardinal))+1.0;
  16. {****************************************************************************
  17. EXTENDED data type routines
  18. ****************************************************************************}
  19. {$define FPC_SYSTEM_HAS_PI}
  20. function fpc_pi_real : valreal;compilerproc;
  21. begin
  22. { Function is handled internal in the compiler }
  23. runerror(207);
  24. result:=0;
  25. end;
  26. {$define FPC_SYSTEM_HAS_ABS}
  27. function fpc_abs_real(d : valreal) : valreal;compilerproc;
  28. begin
  29. { Function is handled internal in the compiler }
  30. runerror(207);
  31. result:=0;
  32. end;
  33. {$define FPC_SYSTEM_HAS_SQR}
  34. function fpc_sqr_real(d : valreal) : valreal;compilerproc;
  35. begin
  36. { Function is handled internal in the compiler }
  37. runerror(207);
  38. result:=0;
  39. end;
  40. const
  41. factor: double = double(int64(1) shl 32);
  42. factor2: double = double(int64(1) shl 31);
  43. {$ifndef FPC_SYSTEM_HAS_TRUNC}
  44. {$define FPC_SYSTEM_HAS_TRUNC}
  45. function fpc_trunc_real(d : valreal) : int64;assembler;compilerproc;
  46. { input: d in fr1 }
  47. { output: result in r3 }
  48. assembler;
  49. var
  50. temp: packed record
  51. case byte of
  52. 0: (l1,l2: longint);
  53. 1: (d: double);
  54. end;
  55. asm
  56. // store d in temp
  57. stfd f1,temp
  58. // extract sign bit (record in cr0)
  59. lwz r3,temp
  60. rlwinm. r3,r3,1,31,31
  61. // make d positive
  62. fabs f1,f1
  63. // load 2^32 in f2
  64. {$if not defined(macos) and not defined(aix)}
  65. {$ifdef FPC_PIC}
  66. {$ifdef darwin}
  67. mflr r0
  68. bcl 20,31,.Lpiclab
  69. .Lpiclab:
  70. mflr r5
  71. mtlr r0
  72. addis r4,r5,(factor-.Lpiclab)@ha
  73. lfd f2,(factor-.Lpiclab)@l(r4)
  74. {$else darwin}
  75. {$error Add pic code for linux/ppc32}
  76. {$endif darwin}
  77. {$else FPC_PIC}
  78. lis r4,factor@ha
  79. lfd f2,factor@l(r4)
  80. {$endif FPC_PIC}
  81. {$else not macos/aix}
  82. lwz r4,factor(r2)
  83. lfd f2,0(r4)
  84. {$endif not macos/aix}
  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. {$if not defined(macos) and not defined(aix)}
  100. {$ifdef FPC_PIC}
  101. {$ifdef darwin}
  102. addis r4,r5,(longint_to_real_helper-.Lpiclab)@ha
  103. lfd f0,(longint_to_real_helper-.Lpiclab)@l(r4)
  104. {$else darwin}
  105. {$error Add pic code for linux/ppc32}
  106. {$endif darwin}
  107. {$else FPC_PIC}
  108. lis r4,longint_to_real_helper@ha
  109. lfd f0,longint_to_real_helper@l(r4)
  110. {$endif FPC_PIC}
  111. {$else not macos/aix}
  112. lwz r4,longint_to_real_helper(r2)
  113. lfd f0,0(r4)
  114. {$endif not macos/aix}
  115. lfd f3,temp
  116. fsub f3,f3,f0
  117. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  118. fnmsub f4,f3,f2,f1
  119. // now, convert to unsigned 32 bit
  120. // load 2^31 in f2
  121. {$if not defined(macos) and not defined(aix)}
  122. {$ifdef FPC_PIC}
  123. {$ifdef darwin}
  124. addis r4,r5,(factor2-.Lpiclab)@ha
  125. lfd f2,(factor2-.Lpiclab)@l(r4)
  126. {$else darwin}
  127. {$error Add pic code for linux/ppc32}
  128. {$endif darwin}
  129. {$else FPC_PIC}
  130. lis r4,factor2@ha
  131. lfd f2,factor2@l(r4)
  132. {$endif FPC_PIC}
  133. {$else not macos/aix}
  134. lwz r4,factor2(r2)
  135. lfd f2,0(r4)
  136. {$endif not macos/aix}
  137. // subtract 2^31
  138. fsub f3,f4,f2
  139. // was the value > 2^31?
  140. fcmpu cr1,f4,f2
  141. // use diff if >= 2^31
  142. fsel f4,f3,f3,f4
  143. // next part same as conversion to signed integer word
  144. fctiwz f4,f4
  145. stfd f4,temp
  146. lwz r4,temp+4
  147. // add 2^31 if value was >=2^31
  148. blt cr1, .LTruncNoAdd
  149. xoris r4,r4,0x8000
  150. .LTruncNoAdd:
  151. // negate value if it was negative to start with
  152. beq cr0,.LTruncPositive
  153. subfic r4,r4,0
  154. subfze r3,r3
  155. .LTruncPositive:
  156. end;
  157. {$endif not FPC_SYSTEM_HAS_TRUNC}
  158. (*
  159. {$ifndef FPC_SYSTEM_HAS_ROUND}
  160. {$define FPC_SYSTEM_HAS_ROUND}
  161. function round(d : extended) : int64;
  162. function fpc_round(d : extended) : int64;assembler;[public, alias:'FPC_ROUND'];compilerproc;
  163. { exactly the same as trunc, except that one fctiwz has become fctiw }
  164. { input: d in fr1 }
  165. { output: result in r3 }
  166. assembler;
  167. var
  168. temp: packed record
  169. case byte of
  170. 0: (l1,l2: longint);
  171. 1: (d: double);
  172. end;
  173. asm
  174. // store d in temp
  175. stfd f1, temp
  176. // extract sign bit (record in cr0)
  177. lwz r4,temp
  178. rlwinm. r4,r4,1,31,31
  179. // make d positive
  180. fabs f1,f1
  181. // load 2^32 in f2
  182. {$if not defined(macos) and not defined(aix)}
  183. lis r4,factor@ha
  184. lfd f2,factor@l(r4)
  185. {$else}
  186. lwz r4,factor(r2)
  187. lfd f2,0(r4)
  188. {$endif}
  189. // check if value is < 0
  190. // f3 := d / 2^32;
  191. fdiv f3,f1,f2
  192. // round
  193. fctiwz f4,f3
  194. // store
  195. stfd f4,temp
  196. // and load into r4
  197. lwz r3,temp+4
  198. // convert back to float
  199. lis r0,0x4330
  200. stw r0,temp
  201. xoris r0,r3,0x8000
  202. stw r0,temp+4
  203. {$if not defined(macos) and not defined(aix)}
  204. lis r4,longint_to_real_helper@ha
  205. lfd f0,longint_to_real_helper@l(r4)
  206. {$else}
  207. lwz r4,longint_to_real_helper(r2)
  208. lfd f0,0(r4)
  209. {$endif}
  210. lfd f3,temp
  211. fsub f3,f3,f0
  212. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  213. fnmsub f4,f3,f2,f1
  214. // now, convert to unsigned 32 bit
  215. // load 2^31 in f2
  216. {$if not defined(macos) and not defined(aix)}
  217. lis r4,factor2@ha
  218. lfd f2,factor2@l(r4)
  219. {$else}
  220. lwz r4,factor2(r2)
  221. lfd f2,0(r4)
  222. {$endif}
  223. // subtract 2^31
  224. fsub f3,f4,f2
  225. // was the value > 2^31?
  226. fcmpu cr1,f4,f2
  227. // use diff if >= 2^31
  228. fsel f4,f3,f3,f4
  229. // next part same as conversion to signed integer word
  230. fctiw f4,f4
  231. stfd f4,temp
  232. lwz r4,temp+4
  233. // add 2^31 if value was >=2^31
  234. blt cr1, .LRoundNoAdd
  235. xoris r4,r4,0x8000
  236. .LRoundNoAdd:
  237. // negate value if it was negative to start with
  238. beq cr0,.LRoundPositive
  239. subfic r4,r4,0
  240. subfze r3,r3
  241. .LRoundPositive:
  242. end;
  243. {$endif not FPC_SYSTEM_HAS_ROUND}
  244. *)
  245. {****************************************************************************
  246. Int to real helpers
  247. ****************************************************************************}
  248. {$ifndef aix}
  249. { these helpers somehow don't seem to work on AIX/Power 7 }
  250. {$define FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  251. function fpc_int64_to_double(i: int64): double; compilerproc;
  252. assembler;
  253. { input: high(i) in r4, low(i) in r3 }
  254. { output: double(i) 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. xoris r3,r3,0x8000
  265. stw r3,temp+4
  266. {$if not defined(macos) and not defined(aix)}
  267. {$ifdef FPC_PIC}
  268. {$ifdef darwin}
  269. mflr r0
  270. bcl 20,31,.Lpiclab
  271. .Lpiclab:
  272. mflr r5
  273. mtlr r0
  274. addis r3,r5,(longint_to_real_helper-.Lpiclab)@ha
  275. lfd f1,(longint_to_real_helper-.Lpiclab)@l(r3)
  276. {$else darwin}
  277. {$error Add pic code for linux/ppc32}
  278. {$endif darwin}
  279. {$else FPC_PIC}
  280. lis r3,longint_to_real_helper@ha
  281. lfd f1,longint_to_real_helper@l(r3)
  282. {$endif FPC_PIC}
  283. {$else not macos/aix}
  284. lwz r3,longint_to_real_helper(r2)
  285. lfd f1,0(r3)
  286. {$endif not mac os}
  287. lfd f0,temp
  288. stw r4,temp+4
  289. fsub f0,f0,f1
  290. {$if not defined(macos) and not defined(aix)}
  291. {$ifdef FPC_PIC}
  292. {$ifdef darwin}
  293. addis r4,r5,(cardinal_to_real_helper-.Lpiclab)@ha
  294. lfd f1,(cardinal_to_real_helper-.Lpiclab)@l(r4)
  295. addis r4,r5,(int_to_real_factor-.Lpiclab)@ha
  296. lfd f3,temp
  297. lfd f2,(int_to_real_factor-.Lpiclab)@l(r4)
  298. {$else darwin}
  299. {$error Add pic code for linux/ppc32}
  300. {$endif darwin}
  301. {$else FPC_PIC}
  302. lis r4,cardinal_to_real_helper@ha
  303. lfd f1,cardinal_to_real_helper@l(r4)
  304. lis r4,int_to_real_factor@ha
  305. lfd f3,temp
  306. lfd f2,int_to_real_factor@l(r4)
  307. {$endif FPC_PIC}
  308. {$else not macos/aix}
  309. lwz r4,cardinal_to_real_helper(r2)
  310. lwz r3,int_to_real_factor(r2)
  311. lfd f3,temp
  312. lfd f1,0(r4)
  313. lfd f2,0(r3)
  314. {$endif not macos/aix}
  315. fsub f3,f3,f1
  316. fmadd f1,f0,f2,f3
  317. end;
  318. {$define FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
  319. function fpc_qword_to_double(q: qword): double; compilerproc;
  320. assembler;
  321. { input: high(q) in r4, low(q) in r3 }
  322. { output: double(q) in f0 }
  323. var
  324. temp: packed record
  325. case byte of
  326. 0: (l1,l2: cardinal);
  327. 1: (d: double);
  328. end;
  329. asm
  330. lis r0,0x4330
  331. stw r0,temp
  332. stw r3,temp+4
  333. lfd f0,temp
  334. {$if not defined(macos) and not defined(aix)}
  335. {$ifdef FPC_PIC}
  336. {$ifdef darwin}
  337. mflr r0
  338. bcl 20,31,.Lpiclab
  339. .Lpiclab:
  340. mflr r5
  341. mtlr r0
  342. addis r3,r5,(cardinal_to_real_helper-.Lpiclab)@ha
  343. lfd f1,(cardinal_to_real_helper-.Lpiclab)@l(r3)
  344. {$else darwin}
  345. {$error Add pic code for linux/ppc32}
  346. {$endif darwin}
  347. {$else FPC_PIC}
  348. lis r3,cardinal_to_real_helper@ha
  349. lfd f1,cardinal_to_real_helper@l(r3)
  350. {$endif FPC_PIC}
  351. {$else not macos/aix}
  352. lwz r3,longint_to_real_helper(r2)
  353. lfd f1,0(r3)
  354. {$endif not macos/aix}
  355. stw r4,temp+4
  356. fsub f0,f0,f1
  357. lfd f3,temp
  358. {$if not defined(macos) and not defined(aix)}
  359. {$ifdef FPC_PIC}
  360. {$ifdef darwin}
  361. addis r4,r5,(int_to_real_factor-.Lpiclab)@ha
  362. lfd f2,(int_to_real_factor-.Lpiclab)@l(r4)
  363. {$else darwin}
  364. {$error Add pic code for linux/ppc32}
  365. {$endif darwin}
  366. {$else FPC_PIC}
  367. lis r4,int_to_real_factor@ha
  368. lfd f2,int_to_real_factor@l(r4)
  369. {$endif FPC_PIC}
  370. {$else not macos/aix}
  371. lwz r4,int_to_real_factor(r2)
  372. lfd f2,0(r4)
  373. {$endif not macos/aix}
  374. fsub f3,f3,f1
  375. fmadd f1,f0,f2,f3
  376. end;
  377. {$endif}