math.inc 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355
  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. {$ifdef INTERNCONSTINTF}
  21. {$define FPC_SYSTEM_HAS_PI}
  22. function fpc_pi_real : valreal;compilerproc;
  23. begin
  24. { Function is handled internal in the compiler }
  25. runerror(207);
  26. result:=0;
  27. end;
  28. {$define FPC_SYSTEM_HAS_ABS}
  29. function fpc_abs_real(d : valreal) : valreal;compilerproc;
  30. begin
  31. { Function is handled internal in the compiler }
  32. runerror(207);
  33. result:=0;
  34. end;
  35. {$define FPC_SYSTEM_HAS_SQR}
  36. function fpc_sqr_real(d : valreal) : valreal;compilerproc;
  37. begin
  38. { Function is handled internal in the compiler }
  39. runerror(207);
  40. result:=0;
  41. end;
  42. {$else}
  43. {$define FPC_SYSTEM_HAS_PI}
  44. function pi : double;[internproc:fpc_in_pi];
  45. {$define FPC_SYSTEM_HAS_ABS}
  46. function abs(d : extended) : extended;[internproc:fpc_in_abs_real];
  47. {$define FPC_SYSTEM_HAS_SQR}
  48. function sqr(d : extended) : extended;[internproc:fpc_in_sqr_real];
  49. {$endif ndef INTERNCONSTINTF}
  50. const
  51. factor: double = double(int64(1) shl 32);
  52. factor2: double = double(int64(1) shl 31);
  53. {$ifndef FPC_SYSTEM_HAS_TRUNC}
  54. {$define FPC_SYSTEM_HAS_TRUNC}
  55. {$ifdef INTERNCONSTINTF}
  56. function fpc_trunc_real(d : valreal) : int64;assembler;compilerproc;
  57. {$else}
  58. function trunc(d : extended) : int64;assembler;[internconst:fpc_in_const_trunc];
  59. {$endif}
  60. { input: d in fr1 }
  61. { output: result in r3 }
  62. assembler;
  63. var
  64. temp: packed record
  65. case byte of
  66. 0: (l1,l2: longint);
  67. 1: (d: double);
  68. end;
  69. asm
  70. // store d in temp
  71. stfd f1,temp
  72. // extract sign bit (record in cr0)
  73. lwz r3,temp
  74. rlwinm. r3,r3,1,31,31
  75. // make d positive
  76. fabs f1,f1
  77. // load 2^32 in f2
  78. {$ifndef macos}
  79. lis r4,factor@ha
  80. lfd f2,factor@l(r4)
  81. {$else}
  82. lwz r4,factor(r2)
  83. lfd f2,0(r4)
  84. {$endif}
  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. lis r4,longint_to_real_helper@ha
  101. lfd f0,longint_to_real_helper@l(r4)
  102. {$else}
  103. lwz r4,longint_to_real_helper(r2)
  104. lfd f0,0(r4)
  105. {$endif}
  106. lfd f3,temp
  107. fsub f3,f3,f0
  108. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  109. fnmsub f4,f3,f2,f1
  110. // now, convert to unsigned 32 bit
  111. // load 2^31 in f2
  112. {$ifndef macos}
  113. lis r4,factor2@ha
  114. lfd f2,factor2@l(r4)
  115. {$else}
  116. lwz r4,factor2(r2)
  117. lfd f2,0(r4)
  118. {$endif}
  119. // subtract 2^31
  120. fsub f3,f4,f2
  121. // was the value > 2^31?
  122. fcmpu cr1,f4,f2
  123. // use diff if >= 2^31
  124. fsel f4,f3,f3,f4
  125. // next part same as conversion to signed integer word
  126. fctiwz f4,f4
  127. stfd f4,temp
  128. lwz r4,temp+4
  129. // add 2^31 if value was >=2^31
  130. blt cr1, .LTruncNoAdd
  131. xoris r4,r4,0x8000
  132. .LTruncNoAdd:
  133. // negate value if it was negative to start with
  134. beq cr0,.LTruncPositive
  135. subfic r4,r4,0
  136. subfze r3,r3
  137. .LTruncPositive:
  138. end;
  139. {$endif not FPC_SYSTEM_HAS_TRUNC}
  140. (*
  141. {$ifndef FPC_SYSTEM_HAS_ROUND}
  142. {$define FPC_SYSTEM_HAS_ROUND}
  143. {$ifdef hascompilerproc}
  144. function round(d : extended) : int64;{$ifndef INTERNCONSTINTF}[internconst:fpc_in_const_round, external name 'FPC_ROUND'];{$endif}
  145. function fpc_round(d : extended) : int64;assembler;[public, alias:'FPC_ROUND'];{$ifdef hascompilerproc}compilerproc;{$endif hascompilerproc}
  146. {$else}
  147. function round(d : extended) : int64;assembler;{$ifndef INTERNCONSTINTF}[internconst:fpc_in_const_round];{$endif}
  148. {$endif hascompilerproc}
  149. { exactly the same as trunc, except that one fctiwz has become fctiw }
  150. { input: d in fr1 }
  151. { output: result in r3 }
  152. assembler;
  153. var
  154. temp: packed record
  155. case byte of
  156. 0: (l1,l2: longint);
  157. 1: (d: double);
  158. end;
  159. asm
  160. // store d in temp
  161. stfd f1, temp
  162. // extract sign bit (record in cr0)
  163. lwz r4,temp
  164. rlwinm. r4,r4,1,31,31
  165. // make d positive
  166. fabs f1,f1
  167. // load 2^32 in f2
  168. {$ifndef macos}
  169. lis r4,factor@ha
  170. lfd f2,factor@l(r4)
  171. {$else}
  172. lwz r4,factor(r2)
  173. lfd f2,0(r4)
  174. {$endif}
  175. // check if value is < 0
  176. // f3 := d / 2^32;
  177. fdiv f3,f1,f2
  178. // round
  179. fctiwz f4,f3
  180. // store
  181. stfd f4,temp
  182. // and load into r4
  183. lwz r3,temp+4
  184. // convert back to float
  185. lis r0,0x4330
  186. stw r0,temp
  187. xoris r0,r3,0x8000
  188. stw r0,temp+4
  189. {$ifndef macos}
  190. lis r4,longint_to_real_helper@ha
  191. lfd f0,longint_to_real_helper@l(r4)
  192. {$else}
  193. lwz r4,longint_to_real_helper(r2)
  194. lfd f0,0(r4)
  195. {$endif}
  196. lfd f3,temp
  197. fsub f3,f3,f0
  198. // f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
  199. fnmsub f4,f3,f2,f1
  200. // now, convert to unsigned 32 bit
  201. // load 2^31 in f2
  202. {$ifndef macos}
  203. lis r4,factor2@ha
  204. lfd f2,factor2@l(r4)
  205. {$else}
  206. lwz r4,factor2(r2)
  207. lfd f2,0(r4)
  208. {$endif}
  209. // subtract 2^31
  210. fsub f3,f4,f2
  211. // was the value > 2^31?
  212. fcmpu cr1,f4,f2
  213. // use diff if >= 2^31
  214. fsel f4,f3,f3,f4
  215. // next part same as conversion to signed integer word
  216. fctiw f4,f4
  217. stfd f4,temp
  218. lwz r4,temp+4
  219. // add 2^31 if value was >=2^31
  220. blt cr1, .LRoundNoAdd
  221. xoris r4,r4,0x8000
  222. .LRoundNoAdd:
  223. // negate value if it was negative to start with
  224. beq cr0,.LRoundPositive
  225. subfic r4,r4,0
  226. subfze r3,r3
  227. .LRoundPositive:
  228. end;
  229. {$endif not FPC_SYSTEM_HAS_ROUND}
  230. *)
  231. {****************************************************************************
  232. Int to real helpers
  233. ****************************************************************************}
  234. {$define FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
  235. function fpc_int64_to_double(i: int64): double; compilerproc;
  236. assembler;
  237. { input: high(i) in r4, low(i) in r3 }
  238. { output: double(i) in f0 }
  239. var
  240. temp: packed record
  241. case byte of
  242. 0: (l1,l2: cardinal);
  243. 1: (d: double);
  244. end;
  245. asm
  246. lis r0,0x4330
  247. stw r0,temp
  248. xoris r3,r3,0x8000
  249. stw r3,temp+4
  250. {$ifndef macos}
  251. lis r3,longint_to_real_helper@ha
  252. lfd f1,longint_to_real_helper@l(r3)
  253. {$else}
  254. lwz r3,longint_to_real_helper(r2)
  255. lfd f1,0(r3)
  256. {$endif}
  257. lfd f0,temp
  258. stw r4,temp+4
  259. fsub f0,f0,f1
  260. {$ifndef macos}
  261. lis r4,cardinal_to_real_helper@ha
  262. lfd f1,cardinal_to_real_helper@l(r4)
  263. lis r4,int_to_real_factor@ha
  264. lfd f3,temp
  265. lfd f2,int_to_real_factor@l(r4)
  266. {$else}
  267. lwz r4,cardinal_to_real_helper(r2)
  268. lwz r3,int_to_real_factor(r2)
  269. lfd f3,temp
  270. lfd f1,0(r4)
  271. lfd f2,0(r3)
  272. {$endif}
  273. fsub f3,f3,f1
  274. fmadd f1,f0,f2,f3
  275. end;
  276. {$define FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
  277. function fpc_qword_to_double(q: qword): double; compilerproc;
  278. assembler;
  279. { input: high(q) in r4, low(q) in r3 }
  280. { output: double(q) in f0 }
  281. var
  282. temp: packed record
  283. case byte of
  284. 0: (l1,l2: cardinal);
  285. 1: (d: double);
  286. end;
  287. asm
  288. lis r0,0x4330
  289. stw r0,temp
  290. stw r3,temp+4
  291. lfd f0,temp
  292. {$ifndef macos}
  293. lis r3,cardinal_to_real_helper@ha
  294. lfd f1,cardinal_to_real_helper@l(r3)
  295. {$else}
  296. lwz r3,longint_to_real_helper(r2)
  297. lfd f1,0(r3)
  298. {$endif}
  299. stw r4,temp+4
  300. fsub f0,f0,f1
  301. lfd f3,temp
  302. {$ifndef macos}
  303. lis r4,int_to_real_factor@ha
  304. lfd f2,int_to_real_factor@l(r4)
  305. {$else}
  306. lwz r4,int_to_real_factor(r2)
  307. lfd f2,0(r4)
  308. {$endif}
  309. fsub f3,f3,f1
  310. fmadd f1,f0,f2,f3
  311. end;
  312. {
  313. $Log$
  314. Revision 1.39 2005-02-14 17:13:31 peter
  315. * truncate log
  316. }