math.inc 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299
  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) : longint;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. function power(bas,expo : longint) : longint;
  114. begin
  115. if bas=0 then
  116. begin
  117. if expo<>0 then
  118. power:=0
  119. else
  120. HandleError(207);
  121. end
  122. else if expo=0 then
  123. power:=1
  124. else
  125. begin
  126. if bas<0 then
  127. begin
  128. if odd(expo) then
  129. power:=-round(exp(ln(-bas)*expo))
  130. else
  131. power:=round(exp(ln(-bas)*expo));
  132. end
  133. else
  134. power:=round(exp(ln(bas)*expo));
  135. end;
  136. end;
  137. {****************************************************************************
  138. Helper routines to support old TP styled reals
  139. ****************************************************************************}
  140. { warning: the following converts a little-endian TP-style real }
  141. { to a big-endian double. So don't byte-swap the TP real! }
  142. {$define FPC_SYSTEM_HAS_REAL2DOUBLE}
  143. function real2double(r : real48) : double;
  144. var
  145. res : array[0..7] of byte;
  146. exponent : word;
  147. begin
  148. { copy mantissa }
  149. res[6]:=0;
  150. res[5]:=r[1] shl 5;
  151. res[4]:=(r[1] shr 3) or (r[2] shl 5);
  152. res[3]:=(r[2] shr 3) or (r[3] shl 5);
  153. res[2]:=(r[3] shr 3) or (r[4] shl 5);
  154. res[1]:=(r[4] shr 3) or (r[5] and $7f) shl 5;
  155. res[0]:=(r[5] and $7f) shr 3;
  156. { copy exponent }
  157. { correct exponent: }
  158. exponent:=(word(r[0])+(1023-129));
  159. res[1]:=res[1] or ((exponent and $f) shl 4);
  160. res[0]:=exponent shr 4;
  161. { set sign }
  162. res[0]:=res[0] or (r[5] and $80);
  163. real2double:=double(res);
  164. end;
  165. {****************************************************************************
  166. Int to real helpers
  167. ****************************************************************************}
  168. const
  169. longint_to_real_helper: int64 = $4330000080000000;
  170. cardinal_to_real_helper: int64 = $430000000000000;
  171. int_to_real_factor: double = double(high(cardinal))+1.0;
  172. function fpc_int64_to_double(i: int64): double; compilerproc;
  173. assembler;
  174. { input: high(i) in r3, low(i) in r4 }
  175. { output: double(i) in f0 }
  176. var
  177. temp: packed record
  178. case byte of
  179. 0: (l1,l2: cardinal);
  180. 1: (d: double);
  181. end;
  182. asm
  183. lis r0,0x4330
  184. stw r0,temp
  185. xoris r3,r3,0x8000
  186. stw r3,4+temp
  187. lis r3,longint_to_real_helper@ha
  188. lfd f1,longint_to_real_helper@l(r3)
  189. lfd f0,temp
  190. stw r4,4+temp
  191. fsub f0,f0,f1
  192. lis r4,cardinal_to_real_helper@ha
  193. lfd f1,cardinal_to_real_helper@l(r4)
  194. lis r3,int_to_real_factor@ha
  195. lfd f3,temp
  196. lfd f2,int_to_real_factor@l(r3)
  197. fsub f3,f3,f1
  198. fmadd f1,f0,f2,f3
  199. end ['R0','R3','R4','F0','F1','F2','F3'];
  200. function fpc_qword_to_double(q: qword): double; compilerproc;
  201. assembler;
  202. { input: high(q) in r3, low(q) in r4 }
  203. { output: double(q) in f0 }
  204. var
  205. temp: packed record
  206. case byte of
  207. 0: (l1,l2: cardinal);
  208. 1: (d: double);
  209. end;
  210. asm
  211. lis r0,0x4330
  212. stw r0,temp
  213. stw r3,4+temp
  214. lfd f0,temp
  215. lis r3,cardinal_to_real_helper@ha
  216. lfd f1,cardinal_to_real_helper@l(r3)
  217. stw r4,4+temp
  218. fsub f0,f0,f1
  219. lfd f3,temp
  220. lis r3,int_to_real_factor@ha
  221. lfd f2,int_to_real_factor@l(r3)
  222. fsub f3,f3,f1
  223. fmadd f1,f0,f2,f3
  224. end ['R0','R3','F0','F1','F2','F3'];
  225. {
  226. $Log$
  227. Revision 1.12 2002-09-08 13:00:21 jonas
  228. * made pi an internproc instead of internconst
  229. Revision 1.11 2002/09/07 16:01:26 peter
  230. * old logs removed and tabs fixed
  231. Revision 1.10 2002/08/18 22:11:10 florian
  232. * fixed remaining assembler errors
  233. Revision 1.9 2002/08/18 21:37:48 florian
  234. * several errors in inline assembler fixed
  235. Revision 1.8 2002/08/10 17:14:36 jonas
  236. * various fixes, mostly changing the names of the modifies registers to
  237. upper case since that seems to be required by the compiler
  238. Revision 1.7 2002/07/31 16:58:12 jonas
  239. * fixed conversion from int64/qword to double errors
  240. Revision 1.6 2002/07/29 21:28:17 florian
  241. * several fixes to get further with linux/ppc system unit compilation
  242. Revision 1.5 2002/07/28 21:39:29 florian
  243. * made abs a compiler proc if it is generic
  244. Revision 1.4 2002/07/28 20:43:49 florian
  245. * several fixes for linux/powerpc
  246. * several fixes to MT
  247. }