math.inc 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409
  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. {$ifndef macos}
  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}
  82. lwz r4,factor(r2)
  83. lfd f2,0(r4)
  84. {$endif not macos}
  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. {$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}
  112. lwz r4,longint_to_real_helper(r2)
  113. lfd f0,0(r4)
  114. {$endif not macos}
  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. {$ifndef macos}
  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}
  134. lwz r4,factor2(r2)
  135. lfd f2,0(r4)
  136. {$endif not macos}
  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. {$ifndef macos}
  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. {$ifndef macos}
  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. {$ifndef macos}
  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. {$define FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  249. function fpc_int64_to_double(i: int64): double; compilerproc;
  250. assembler;
  251. { input: high(i) in r4, low(i) in r3 }
  252. { output: double(i) in f0 }
  253. var
  254. temp: packed record
  255. case byte of
  256. 0: (l1,l2: cardinal);
  257. 1: (d: double);
  258. end;
  259. asm
  260. lis r0,0x4330
  261. stw r0,temp
  262. xoris r3,r3,0x8000
  263. stw r3,temp+4
  264. {$ifndef macos}
  265. {$ifdef FPC_PIC}
  266. {$ifdef darwin}
  267. mflr r0
  268. bcl 20,31,.Lpiclab
  269. .Lpiclab:
  270. mflr r5
  271. mtlr r0
  272. addis r3,r5,(longint_to_real_helper-.Lpiclab)@ha
  273. lfd f1,(longint_to_real_helper-.Lpiclab)@l(r3)
  274. {$else darwin}
  275. {$error Add pic code for linux/ppc32}
  276. {$endif darwin}
  277. {$else FPC_PIC}
  278. lis r3,longint_to_real_helper@ha
  279. lfd f1,longint_to_real_helper@l(r3)
  280. {$endif FPC_PIC}
  281. {$else not macos}
  282. lwz r3,longint_to_real_helper(r2)
  283. lfd f1,0(r3)
  284. {$endif not mac os}
  285. lfd f0,temp
  286. stw r4,temp+4
  287. fsub f0,f0,f1
  288. {$ifndef macos}
  289. {$ifdef FPC_PIC}
  290. {$ifdef darwin}
  291. addis r4,r5,(cardinal_to_real_helper-.Lpiclab)@ha
  292. lfd f1,(cardinal_to_real_helper-.Lpiclab)@l(r4)
  293. addis r4,r5,(int_to_real_factor-.Lpiclab)@ha
  294. lfd f3,temp
  295. lfd f2,(int_to_real_factor-.Lpiclab)@l(r4)
  296. {$else darwin}
  297. {$error Add pic code for linux/ppc32}
  298. {$endif darwin}
  299. {$else FPC_PIC}
  300. lis r4,cardinal_to_real_helper@ha
  301. lfd f1,cardinal_to_real_helper@l(r4)
  302. lis r4,int_to_real_factor@ha
  303. lfd f3,temp
  304. lfd f2,int_to_real_factor@l(r4)
  305. {$endif FPC_PIC}
  306. {$else not macos}
  307. lwz r4,cardinal_to_real_helper(r2)
  308. lwz r3,int_to_real_factor(r2)
  309. lfd f3,temp
  310. lfd f1,0(r4)
  311. lfd f2,0(r3)
  312. {$endif not macos}
  313. fsub f3,f3,f1
  314. fmadd f1,f0,f2,f3
  315. end;
  316. {$define FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
  317. function fpc_qword_to_double(q: qword): double; compilerproc;
  318. assembler;
  319. { input: high(q) in r4, low(q) in r3 }
  320. { output: double(q) in f0 }
  321. var
  322. temp: packed record
  323. case byte of
  324. 0: (l1,l2: cardinal);
  325. 1: (d: double);
  326. end;
  327. asm
  328. lis r0,0x4330
  329. stw r0,temp
  330. stw r3,temp+4
  331. lfd f0,temp
  332. {$ifndef macos}
  333. {$ifdef FPC_PIC}
  334. {$ifdef darwin}
  335. mflr r0
  336. bcl 20,31,.Lpiclab
  337. .Lpiclab:
  338. mflr r5
  339. mtlr r0
  340. addis r3,r5,(cardinal_to_real_helper-.Lpiclab)@ha
  341. lfd f1,(cardinal_to_real_helper-.Lpiclab)@l(r3)
  342. {$else darwin}
  343. {$error Add pic code for linux/ppc32}
  344. {$endif darwin}
  345. {$else FPC_PIC}
  346. lis r3,cardinal_to_real_helper@ha
  347. lfd f1,cardinal_to_real_helper@l(r3)
  348. {$endif FPC_PIC}
  349. {$else not macos}
  350. lwz r3,longint_to_real_helper(r2)
  351. lfd f1,0(r3)
  352. {$endif not macos}
  353. stw r4,temp+4
  354. fsub f0,f0,f1
  355. lfd f3,temp
  356. {$ifndef macos}
  357. {$ifdef FPC_PIC}
  358. {$ifdef darwin}
  359. addis r4,r5,(int_to_real_factor-.Lpiclab)@ha
  360. lfd f2,(int_to_real_factor-.Lpiclab)@l(r4)
  361. {$else darwin}
  362. {$error Add pic code for linux/ppc32}
  363. {$endif darwin}
  364. {$else FPC_PIC}
  365. lis r4,int_to_real_factor@ha
  366. lfd f2,int_to_real_factor@l(r4)
  367. {$endif FPC_PIC}
  368. {$else not macos}
  369. lwz r4,int_to_real_factor(r2)
  370. lfd f2,0(r4)
  371. {$endif not macos}
  372. fsub f3,f3,f1
  373. fmadd f1,f0,f2,f3
  374. end;