math.inc 9.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326
  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 mathamatical 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. {****************************************************************************
  14. EXTENDED data type routines
  15. ****************************************************************************}
  16. {$define FPC_SYSTEM_HAS_PI}
  17. function pi : double;[internproc:in_pi];
  18. {$define FPC_SYSTEM_HAS_ABS}
  19. function abs(d : extended) : extended;[internproc:in_abs_extended];
  20. {$define FPC_SYSTEM_HAS_SQR}
  21. function sqr(d : extended) : extended;[internproc:in_sqr_extended];
  22. {$define FPC_SYSTEM_HAS_SQRT}
  23. function sqrt(d : extended) : extended;[internproc:in_sqrt_extended];
  24. {
  25. function arctan(d : extended) : extended;[internconst:in_arctan_extended];
  26. begin
  27. runerror(207);
  28. end;
  29. function ln(d : extended) : extended;[internconst:in_ln_extended];
  30. begin
  31. runerror(207);
  32. end;
  33. function sin(d : extended) : extended;[internconst: in_sin_extended];
  34. begin
  35. runerror(207);
  36. end;
  37. function cos(d : extended) : extended;[internconst:in_cos_extended];
  38. begin
  39. runerror(207);
  40. end;
  41. function exp(d : extended) : extended;[internconst:in_const_exp];
  42. begin
  43. runerror(207);
  44. end;
  45. function frac(d : extended) : extended;[internconst:in_const_frac];
  46. begin
  47. runerror(207);
  48. end;
  49. }
  50. {$define FPC_SYSTEM_HAS_INT}
  51. {$warning FIX ME}
  52. function int(d : extended) : extended;[internconst:in_const_int];
  53. begin
  54. runerror(207);
  55. end;
  56. {$define FPC_SYSTEM_HAS_TRUNC}
  57. {$warning FIX ME}
  58. function trunc(d : extended) : int64;assembler;[internconst:in_const_trunc];
  59. { input: d in fr1 }
  60. { output: result in r3 }
  61. assembler;
  62. var
  63. temp: packed record
  64. case byte of
  65. 0: (l1,l2: longint);
  66. 1: (d: double);
  67. end;
  68. asm
  69. fctiwz f1,f1
  70. stfd f1,temp
  71. lwz r3,temp
  72. lwz r4,4+temp
  73. end ['R3','F1'];
  74. {$define FPC_SYSTEM_HAS_ROUND}
  75. function round(d : extended) : int64;assembler;[internconst:in_const_round];
  76. { input: d in fr1 }
  77. { output: result in r3 }
  78. assembler;
  79. var
  80. temp : packed record
  81. case byte of
  82. 0: (l1,l2: longint);
  83. 1: (d: double);
  84. end;
  85. asm
  86. fctiw f1,f1
  87. stfd f1,temp
  88. lwz r3,temp
  89. lwz r4,4+temp
  90. end ['R3','F1'];
  91. {$define FPC_SYSTEM_HAS_POWER}
  92. function power(bas,expo : extended) : extended;
  93. begin
  94. if bas=0 then
  95. begin
  96. if expo<>0 then
  97. power:=0.0
  98. else
  99. HandleError(207);
  100. end
  101. else if expo=0 then
  102. power:=1
  103. else
  104. { bas < 0 is not allowed }
  105. if bas<0 then
  106. handleerror(207)
  107. else
  108. power:=exp(ln(bas)*expo);
  109. end;
  110. {****************************************************************************
  111. Longint data type routines
  112. ****************************************************************************}
  113. {$define FPC_SYSTEM_HAS_POWER_INT64}
  114. function power(bas,expo : int64) : int64;
  115. begin
  116. if bas=0 then
  117. begin
  118. if expo<>0 then
  119. power:=0
  120. else
  121. HandleError(207);
  122. end
  123. else if expo=0 then
  124. power:=1
  125. else
  126. begin
  127. if bas<0 then
  128. begin
  129. if odd(expo) then
  130. power:=-round(exp(ln(-bas)*expo))
  131. else
  132. power:=round(exp(ln(-bas)*expo));
  133. end
  134. else
  135. power:=round(exp(ln(bas)*expo));
  136. end;
  137. end;
  138. {****************************************************************************
  139. Helper routines to support old TP styled reals
  140. ****************************************************************************}
  141. { warning: the following converts a little-endian TP-style real }
  142. { to a big-endian double. So don't byte-swap the TP real! }
  143. {$define FPC_SYSTEM_HAS_REAL2DOUBLE}
  144. function real2double(r : real48) : double;
  145. var
  146. res : array[0..7] of byte;
  147. exponent : word;
  148. begin
  149. { copy mantissa }
  150. res[6]:=0;
  151. res[5]:=r[1] shl 5;
  152. res[4]:=(r[1] shr 3) or (r[2] shl 5);
  153. res[3]:=(r[2] shr 3) or (r[3] shl 5);
  154. res[2]:=(r[3] shr 3) or (r[4] shl 5);
  155. res[1]:=(r[4] shr 3) or (r[5] and $7f) shl 5;
  156. res[0]:=(r[5] and $7f) shr 3;
  157. { copy exponent }
  158. { correct exponent: }
  159. exponent:=(word(r[0])+(1023-129));
  160. res[1]:=res[1] or ((exponent and $f) shl 4);
  161. res[0]:=exponent shr 4;
  162. { set sign }
  163. res[0]:=res[0] or (r[5] and $80);
  164. real2double:=double(res);
  165. end;
  166. {****************************************************************************
  167. Int to real helpers
  168. ****************************************************************************}
  169. const
  170. longint_to_real_helper: int64 = $4330000080000000;
  171. cardinal_to_real_helper: int64 = $430000000000000;
  172. int_to_real_factor: double = double(high(cardinal))+1.0;
  173. function fpc_int64_to_double(i: int64): double; compilerproc;
  174. assembler;
  175. { input: high(i) in r3, low(i) in r4 }
  176. { output: double(i) in f0 }
  177. var
  178. temp: packed record
  179. case byte of
  180. 0: (l1,l2: cardinal);
  181. 1: (d: double);
  182. end;
  183. asm
  184. lis r0,0x4330
  185. stw r0,temp
  186. xoris r3,r3,0x8000
  187. stw r3,4+temp
  188. {$ifndef macos}
  189. lis r3,longint_to_real_helper@ha
  190. lfd f1,longint_to_real_helper@l(r3)
  191. {$else}
  192. lfd f1,longint_to_real_helper(r2)
  193. {$endif}
  194. lfd f0,temp
  195. stw r4,4+temp
  196. fsub f0,f0,f1
  197. {$ifndef macos}
  198. lis r4,cardinal_to_real_helper@ha
  199. lfd f1,cardinal_to_real_helper@l(r4)
  200. lis r3,int_to_real_factor@ha
  201. lfd f3,temp
  202. lfd f2,int_to_real_factor@l(r3)
  203. {$else}
  204. lfd f1,cardinal_to_real_helper(r2)
  205. lfd f3,temp
  206. lfd f2,int_to_real_factor(r2)
  207. {$endif}
  208. fsub f3,f3,f1
  209. fmadd f1,f0,f2,f3
  210. end ['R0','R3','R4','F0','F1','F2','F3'];
  211. function fpc_qword_to_double(q: qword): double; compilerproc;
  212. assembler;
  213. { input: high(q) in r3, low(q) in r4 }
  214. { output: double(q) in f0 }
  215. var
  216. temp: packed record
  217. case byte of
  218. 0: (l1,l2: cardinal);
  219. 1: (d: double);
  220. end;
  221. asm
  222. lis r0,0x4330
  223. stw r0,temp
  224. stw r3,4+temp
  225. lfd f0,temp
  226. {$ifndef macos}
  227. lis r3,cardinal_to_real_helper@ha
  228. lfd f1,cardinal_to_real_helper@l(r3)
  229. {$else}
  230. lfd f1,cardinal_to_real_helper(r2)
  231. {$endif}
  232. stw r4,4+temp
  233. fsub f0,f0,f1
  234. lfd f3,temp
  235. {$ifndef macos}
  236. lis r3,int_to_real_factor@ha
  237. lfd f2,int_to_real_factor@l(r3)
  238. {$else}
  239. lfd f2,int_to_real_factor(r2)
  240. {$endif}
  241. fsub f3,f3,f1
  242. fmadd f1,f0,f2,f3
  243. end ['R0','R3','F0','F1','F2','F3'];
  244. {
  245. $Log$
  246. Revision 1.15 2003-01-15 01:09:04 florian
  247. * changed power(...) prototype to int64
  248. Revision 1.14 2002/11/28 11:04:16 olle
  249. * macos: refs to globals in asm adapted to macos
  250. Revision 1.13 2002/10/21 18:08:28 jonas
  251. * round has int64 instead of longint result
  252. Revision 1.12 2002/09/08 13:00:21 jonas
  253. * made pi an internproc instead of internconst
  254. Revision 1.11 2002/09/07 16:01:26 peter
  255. * old logs removed and tabs fixed
  256. Revision 1.10 2002/08/18 22:11:10 florian
  257. * fixed remaining assembler errors
  258. Revision 1.9 2002/08/18 21:37:48 florian
  259. * several errors in inline assembler fixed
  260. Revision 1.8 2002/08/10 17:14:36 jonas
  261. * various fixes, mostly changing the names of the modifies registers to
  262. upper case since that seems to be required by the compiler
  263. Revision 1.7 2002/07/31 16:58:12 jonas
  264. * fixed conversion from int64/qword to double errors
  265. Revision 1.6 2002/07/29 21:28:17 florian
  266. * several fixes to get further with linux/ppc system unit compilation
  267. Revision 1.5 2002/07/28 21:39:29 florian
  268. * made abs a compiler proc if it is generic
  269. Revision 1.4 2002/07/28 20:43:49 florian
  270. * several fixes for linux/powerpc
  271. * several fixes to MT
  272. }