fp_log.c 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219
  1. /*
  2. fp_trig.c: floating-point math routines for the Linux-m68k
  3. floating point emulator.
  4. Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.
  5. I hereby give permission, free of charge, to copy, modify, and
  6. redistribute this software, in source or binary form, provided that
  7. the above copyright notice and the following disclaimer are included
  8. in all such copies.
  9. THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
  10. OR IMPLIED.
  11. */
  12. #include "fp_emu.h"
  13. static const struct fp_ext fp_one =
  14. {
  15. .exp = 0x3fff,
  16. };
  17. extern struct fp_ext *fp_fadd(struct fp_ext *dest, const struct fp_ext *src);
  18. extern struct fp_ext *fp_fdiv(struct fp_ext *dest, const struct fp_ext *src);
  19. struct fp_ext *
  20. fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
  21. {
  22. struct fp_ext tmp, src2;
  23. int i, exp;
  24. dprint(PINSTR, "fsqrt\n");
  25. fp_monadic_check(dest, src);
  26. if (IS_ZERO(dest))
  27. return dest;
  28. if (dest->sign) {
  29. fp_set_nan(dest);
  30. return dest;
  31. }
  32. if (IS_INF(dest))
  33. return dest;
  34. /*
  35. * sqrt(m) * 2^(p) , if e = 2*p
  36. * sqrt(m*2^e) =
  37. * sqrt(2*m) * 2^(p) , if e = 2*p + 1
  38. *
  39. * So we use the last bit of the exponent to decide whether to
  40. * use the m or 2*m.
  41. *
  42. * Since only the fractional part of the mantissa is stored and
  43. * the integer part is assumed to be one, we place a 1 or 2 into
  44. * the fixed point representation.
  45. */
  46. exp = dest->exp;
  47. dest->exp = 0x3FFF;
  48. if (!(exp & 1)) /* lowest bit of exponent is set */
  49. dest->exp++;
  50. fp_copy_ext(&src2, dest);
  51. /*
  52. * The taylor row around a for sqrt(x) is:
  53. * sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
  54. * With a=1 this gives:
  55. * sqrt(x) = 1 + 1/2*(x-1)
  56. * = 1/2*(1+x)
  57. */
  58. fp_fadd(dest, &fp_one);
  59. dest->exp--; /* * 1/2 */
  60. /*
  61. * We now apply the newton rule to the function
  62. * f(x) := x^2 - r
  63. * which has a null point on x = sqrt(r).
  64. *
  65. * It gives:
  66. * x' := x - f(x)/f'(x)
  67. * = x - (x^2 -r)/(2*x)
  68. * = x - (x - r/x)/2
  69. * = (2*x - x + r/x)/2
  70. * = (x + r/x)/2
  71. */
  72. for (i = 0; i < 9; i++) {
  73. fp_copy_ext(&tmp, &src2);
  74. fp_fdiv(&tmp, dest);
  75. fp_fadd(dest, &tmp);
  76. dest->exp--;
  77. }
  78. dest->exp += (exp - 0x3FFF) / 2;
  79. return dest;
  80. }
  81. struct fp_ext *
  82. fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
  83. {
  84. uprint("fetoxm1\n");
  85. fp_monadic_check(dest, src);
  86. return dest;
  87. }
  88. struct fp_ext *
  89. fp_fetox(struct fp_ext *dest, struct fp_ext *src)
  90. {
  91. uprint("fetox\n");
  92. fp_monadic_check(dest, src);
  93. return dest;
  94. }
  95. struct fp_ext *
  96. fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
  97. {
  98. uprint("ftwotox\n");
  99. fp_monadic_check(dest, src);
  100. return dest;
  101. }
  102. struct fp_ext *
  103. fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
  104. {
  105. uprint("ftentox\n");
  106. fp_monadic_check(dest, src);
  107. return dest;
  108. }
  109. struct fp_ext *
  110. fp_flogn(struct fp_ext *dest, struct fp_ext *src)
  111. {
  112. uprint("flogn\n");
  113. fp_monadic_check(dest, src);
  114. return dest;
  115. }
  116. struct fp_ext *
  117. fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
  118. {
  119. uprint("flognp1\n");
  120. fp_monadic_check(dest, src);
  121. return dest;
  122. }
  123. struct fp_ext *
  124. fp_flog10(struct fp_ext *dest, struct fp_ext *src)
  125. {
  126. uprint("flog10\n");
  127. fp_monadic_check(dest, src);
  128. return dest;
  129. }
  130. struct fp_ext *
  131. fp_flog2(struct fp_ext *dest, struct fp_ext *src)
  132. {
  133. uprint("flog2\n");
  134. fp_monadic_check(dest, src);
  135. return dest;
  136. }
  137. struct fp_ext *
  138. fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
  139. {
  140. dprint(PINSTR, "fgetexp\n");
  141. fp_monadic_check(dest, src);
  142. if (IS_INF(dest)) {
  143. fp_set_nan(dest);
  144. return dest;
  145. }
  146. if (IS_ZERO(dest))
  147. return dest;
  148. fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
  149. fp_normalize_ext(dest);
  150. return dest;
  151. }
  152. struct fp_ext *
  153. fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
  154. {
  155. dprint(PINSTR, "fgetman\n");
  156. fp_monadic_check(dest, src);
  157. if (IS_ZERO(dest))
  158. return dest;
  159. if (IS_INF(dest))
  160. return dest;
  161. dest->exp = 0x3FFF;
  162. return dest;
  163. }