math.inc 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322
  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. 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. {$ifndef macos}
  188. lis r3,longint_to_real_helper@ha
  189. lfd f1,longint_to_real_helper@l(r3)
  190. {$else}
  191. lfd f1,longint_to_real_helper(r2)
  192. {$endif}
  193. lfd f0,temp
  194. stw r4,4+temp
  195. fsub f0,f0,f1
  196. {$ifndef macos}
  197. lis r4,cardinal_to_real_helper@ha
  198. lfd f1,cardinal_to_real_helper@l(r4)
  199. lis r3,int_to_real_factor@ha
  200. lfd f3,temp
  201. lfd f2,int_to_real_factor@l(r3)
  202. {$else}
  203. lfd f1,cardinal_to_real_helper(r2)
  204. lfd f3,temp
  205. lfd f2,int_to_real_factor(r2)
  206. {$endif}
  207. fsub f3,f3,f1
  208. fmadd f1,f0,f2,f3
  209. end ['R0','R3','R4','F0','F1','F2','F3'];
  210. function fpc_qword_to_double(q: qword): double; compilerproc;
  211. assembler;
  212. { input: high(q) in r3, low(q) in r4 }
  213. { output: double(q) in f0 }
  214. var
  215. temp: packed record
  216. case byte of
  217. 0: (l1,l2: cardinal);
  218. 1: (d: double);
  219. end;
  220. asm
  221. lis r0,0x4330
  222. stw r0,temp
  223. stw r3,4+temp
  224. lfd f0,temp
  225. {$ifndef macos}
  226. lis r3,cardinal_to_real_helper@ha
  227. lfd f1,cardinal_to_real_helper@l(r3)
  228. {$else}
  229. lfd f1,cardinal_to_real_helper(r2)
  230. {$endif}
  231. stw r4,4+temp
  232. fsub f0,f0,f1
  233. lfd f3,temp
  234. {$ifndef macos}
  235. lis r3,int_to_real_factor@ha
  236. lfd f2,int_to_real_factor@l(r3)
  237. {$else}
  238. lfd f2,int_to_real_factor(r2)
  239. {$endif}
  240. fsub f3,f3,f1
  241. fmadd f1,f0,f2,f3
  242. end ['R0','R3','F0','F1','F2','F3'];
  243. {
  244. $Log$
  245. Revision 1.14 2002-11-28 11:04:16 olle
  246. * macos: refs to globals in asm adapted to macos
  247. Revision 1.13 2002/10/21 18:08:28 jonas
  248. * round has int64 instead of longint result
  249. Revision 1.12 2002/09/08 13:00:21 jonas
  250. * made pi an internproc instead of internconst
  251. Revision 1.11 2002/09/07 16:01:26 peter
  252. * old logs removed and tabs fixed
  253. Revision 1.10 2002/08/18 22:11:10 florian
  254. * fixed remaining assembler errors
  255. Revision 1.9 2002/08/18 21:37:48 florian
  256. * several errors in inline assembler fixed
  257. Revision 1.8 2002/08/10 17:14:36 jonas
  258. * various fixes, mostly changing the names of the modifies registers to
  259. upper case since that seems to be required by the compiler
  260. Revision 1.7 2002/07/31 16:58:12 jonas
  261. * fixed conversion from int64/qword to double errors
  262. Revision 1.6 2002/07/29 21:28:17 florian
  263. * several fixes to get further with linux/ppc system unit compilation
  264. Revision 1.5 2002/07/28 21:39:29 florian
  265. * made abs a compiler proc if it is generic
  266. Revision 1.4 2002/07/28 20:43:49 florian
  267. * several fixes for linux/powerpc
  268. * several fixes to MT
  269. }