2
0

math.inc 12 KB

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