math.inc 8.6 KB

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