math.inc 8.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308
  1. {
  2. $Id$
  3. This file is part of the Free Pascal run time library.
  4. Copyright (c) 1999-2000 by the Free Pascal development team
  5. Implementation of mathamatical Routines (only for real)
  6. See the file COPYING.FPC, included in this distribution,
  7. for details about the copyright.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  11. **********************************************************************}
  12. {****************************************************************************
  13. EXTENDED data type routines
  14. ****************************************************************************}
  15. function pi : extended;[internproc:in_pi];
  16. function abs(d : extended) : extended;[internproc:in_abs_extended];
  17. function sqr(d : extended) : extended;[internproc:in_sqr_extended];
  18. function sqrt(d : extended) : extended;[internproc:in_sqrt_extended];
  19. function arctan(d : extended) : extended;[internproc:in_arctan_extended];
  20. function ln(d : extended) : extended;[internproc:in_ln_extended];
  21. function sin(d : extended) : extended;[internproc:in_sin_extended];
  22. function cos(d : extended) : extended;[internproc:in_cos_extended];
  23. function exp(d : extended) : extended;assembler;[internconst:in_const_exp];
  24. asm
  25. // comes from DJ GPP
  26. fldt d
  27. fldl2e
  28. fmulp %st(1)
  29. fstcw .LCW1
  30. fstcw .LCW2
  31. andw $0xf3ff,.LCW2
  32. orw $0x0400,.LCW2
  33. fldcw .LCW2
  34. fld %st(0)
  35. frndint
  36. fldcw .LCW1
  37. fxch %st(1)
  38. fsub %st(1),%st
  39. f2xm1
  40. fld1
  41. faddp %st(1)
  42. fscale
  43. fstp %st(1)
  44. jmp .LCW3
  45. // store some help data in the data segment
  46. .data
  47. .LCW1:
  48. .word 0
  49. .LCW2:
  50. .word 0
  51. .text
  52. .LCW3:
  53. end;
  54. function frac(d : extended) : extended;assembler;[internconst:in_const_frac];
  55. asm
  56. subl $16,%esp
  57. fnstcw -4(%ebp)
  58. fwait
  59. movw -4(%ebp),%cx
  60. orw $0x0c3f,%cx
  61. movw %cx,-8(%ebp)
  62. fldcw -8(%ebp)
  63. fwait
  64. fldt d
  65. frndint
  66. fldt d
  67. fsub %st(1)
  68. fstp %st(1)
  69. fclex
  70. fldcw -4(%ebp)
  71. end ['ECX'];
  72. function int(d : extended) : extended;assembler;[internconst:in_const_int];
  73. asm
  74. subl $16,%esp
  75. fnstcw -4(%ebp)
  76. fwait
  77. movw -4(%ebp),%cx
  78. orw $0x0c3f,%cx
  79. movw %cx,-8(%ebp)
  80. fldcw -8(%ebp)
  81. fwait
  82. fldt d
  83. frndint
  84. fclex
  85. fldcw -4(%ebp)
  86. end ['ECX'];
  87. function trunc(d : extended) : longint;assembler;[internconst:in_const_trunc];
  88. asm
  89. subl $16,%esp
  90. fnstcw -4(%ebp)
  91. fwait
  92. movw -4(%ebp),%cx
  93. orw $0x0c3f,%cx
  94. movw %cx,-8(%ebp)
  95. fldcw -8(%ebp)
  96. fwait
  97. fldt d
  98. fistpl -8(%ebp)
  99. movl -8(%ebp),%eax
  100. fldcw -4(%ebp)
  101. end ['EAX','ECX'];
  102. function round(d : extended) : longint;assembler;[internconst:in_const_round];
  103. asm
  104. subl $8,%esp
  105. fnstcw -4(%ebp)
  106. fwait
  107. movw $0x1372,-8(%ebp)
  108. fldcw -8(%ebp)
  109. fwait
  110. fldt d
  111. fistpl -8(%ebp)
  112. movl -8(%ebp),%eax
  113. fldcw -4(%ebp)
  114. end ['EAX','ECX'];
  115. function power(bas,expo : extended) : extended;
  116. begin
  117. if bas=0 then
  118. begin
  119. if expo<>0 then
  120. power:=0.0
  121. else
  122. HandleError(207);
  123. end
  124. else if expo=0 then
  125. power:=1
  126. else
  127. { bas < 0 is not allowed }
  128. if bas<0 then
  129. handleerror(207)
  130. else
  131. power:=exp(ln(bas)*expo);
  132. end;
  133. {****************************************************************************
  134. Longint data type routines
  135. ****************************************************************************}
  136. function power(bas,expo : longint) : longint;
  137. begin
  138. if bas=0 then
  139. begin
  140. if expo<>0 then
  141. power:=0
  142. else
  143. HandleError(207);
  144. end
  145. else if expo=0 then
  146. power:=1
  147. else
  148. begin
  149. if bas<0 then
  150. begin
  151. if odd(expo) then
  152. power:=-round(exp(ln(-bas)*expo))
  153. else
  154. power:=round(exp(ln(-bas)*expo));
  155. end
  156. else
  157. power:=round(exp(ln(bas)*expo));
  158. end;
  159. end;
  160. {****************************************************************************
  161. Fixed data type routines
  162. ****************************************************************************}
  163. {$ifdef HASFIXED} { Not yet allowed }
  164. function sqrt(d : fixed) : fixed;
  165. begin
  166. asm
  167. movl d,%eax
  168. movl %eax,%ebx
  169. movl %eax,%ecx
  170. jecxz .L_kl
  171. xorl %esi,%esi
  172. .L_it:
  173. xorl %edx,%edx
  174. idivl %ebx
  175. addl %ebx,%eax
  176. shrl $1,%eax
  177. subl %eax,%esi
  178. cmpl $1,%esi
  179. jbe .L_kl
  180. movl %eax,%esi
  181. movl %eax,%ebx
  182. movl %ecx,%eax
  183. jmp .L_it
  184. .L_kl:
  185. shl $8,%eax
  186. leave
  187. ret $4
  188. end;
  189. end;
  190. function int(d : fixed) : fixed;
  191. {*****************************************************************}
  192. { Returns the integral part of d }
  193. {*****************************************************************}
  194. begin
  195. int:=d and $ffff0000; { keep only upper bits }
  196. end;
  197. function trunc(d : fixed) : longint;
  198. {*****************************************************************}
  199. { Returns the Truncated integral part of d }
  200. {*****************************************************************}
  201. begin
  202. trunc:=longint(integer(d shr 16)); { keep only upper 16 bits }
  203. end;
  204. function frac(d : fixed) : fixed;
  205. {*****************************************************************}
  206. { Returns the Fractional part of d }
  207. {*****************************************************************}
  208. begin
  209. frac:=d AND $ffff; { keep only decimal parts - lower 16 bits }
  210. end;
  211. function abs(d : fixed) : fixed;
  212. {*****************************************************************}
  213. { Returns the Absolute value of d }
  214. {*****************************************************************}
  215. begin
  216. asm
  217. movl d,%eax
  218. rol $16,%eax { Swap high & low word.}
  219. {Absolute value: Invert all bits and increment when <0 .}
  220. cwd { When ax<0, dx contains $ffff}
  221. xorw %dx,%ax { Inverts all bits when dx=$ffff.}
  222. subw %dx,%ax { Increments when dx=$ffff.}
  223. rol $16,%eax { Swap high & low word.}
  224. leave
  225. ret $4
  226. end;
  227. end;
  228. function sqr(d : fixed) : fixed;
  229. {*****************************************************************}
  230. { Returns the Absolute squared value of d }
  231. {*****************************************************************}
  232. begin
  233. {16-bit precision needed, not 32 =)}
  234. sqr := d*d;
  235. { sqr := (d SHR 8 * d) SHR 8; }
  236. end;
  237. function Round(x: fixed): longint;
  238. {*****************************************************************}
  239. { Returns the Rounded value of d as a longint }
  240. {*****************************************************************}
  241. var
  242. lowf:integer;
  243. highf:integer;
  244. begin
  245. lowf:=x and $ffff; { keep decimal part ... }
  246. highf :=integer(x shr 16);
  247. if lowf > 5 then
  248. highf:=highf+1
  249. else
  250. if lowf = 5 then
  251. begin
  252. { here we must check the sign ... }
  253. { if greater or equal to zero, then }
  254. { greater value will be found by adding }
  255. { one... }
  256. if highf >= 0 then
  257. Highf:=Highf+1;
  258. end;
  259. Round:= longint(highf);
  260. end;
  261. {$endif HASFIXED}
  262. {
  263. $Log$
  264. Revision 1.3 2000-07-14 10:33:10 michael
  265. + Conditionals fixed
  266. Revision 1.2 2000/07/13 11:33:41 michael
  267. + removed logs
  268. }