sqrt.c 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185
  1. /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
  2. /*
  3. * ====================================================
  4. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  5. *
  6. * Developed at SunSoft, a Sun Microsystems, Inc. business.
  7. * Permission to use, copy, modify, and distribute this
  8. * software is freely granted, provided that this notice
  9. * is preserved.
  10. * ====================================================
  11. */
  12. /* sqrt(x)
  13. * Return correctly rounded sqrt.
  14. * ------------------------------------------
  15. * | Use the hardware sqrt if you have one |
  16. * ------------------------------------------
  17. * Method:
  18. * Bit by bit method using integer arithmetic. (Slow, but portable)
  19. * 1. Normalization
  20. * Scale x to y in [1,4) with even powers of 2:
  21. * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
  22. * sqrt(x) = 2^k * sqrt(y)
  23. * 2. Bit by bit computation
  24. * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
  25. * i 0
  26. * i+1 2
  27. * s = 2*q , and y = 2 * ( y - q ). (1)
  28. * i i i i
  29. *
  30. * To compute q from q , one checks whether
  31. * i+1 i
  32. *
  33. * -(i+1) 2
  34. * (q + 2 ) <= y. (2)
  35. * i
  36. * -(i+1)
  37. * If (2) is false, then q = q ; otherwise q = q + 2 .
  38. * i+1 i i+1 i
  39. *
  40. * With some algebric manipulation, it is not difficult to see
  41. * that (2) is equivalent to
  42. * -(i+1)
  43. * s + 2 <= y (3)
  44. * i i
  45. *
  46. * The advantage of (3) is that s and y can be computed by
  47. * i i
  48. * the following recurrence formula:
  49. * if (3) is false
  50. *
  51. * s = s , y = y ; (4)
  52. * i+1 i i+1 i
  53. *
  54. * otherwise,
  55. * -i -(i+1)
  56. * s = s + 2 , y = y - s - 2 (5)
  57. * i+1 i i+1 i i
  58. *
  59. * One may easily use induction to prove (4) and (5).
  60. * Note. Since the left hand side of (3) contain only i+2 bits,
  61. * it does not necessary to do a full (53-bit) comparison
  62. * in (3).
  63. * 3. Final rounding
  64. * After generating the 53 bits result, we compute one more bit.
  65. * Together with the remainder, we can decide whether the
  66. * result is exact, bigger than 1/2ulp, or less than 1/2ulp
  67. * (it will never equal to 1/2ulp).
  68. * The rounding mode can be detected by checking whether
  69. * huge + tiny is equal to huge, and whether huge - tiny is
  70. * equal to huge for some floating point number "huge" and "tiny".
  71. *
  72. * Special cases:
  73. * sqrt(+-0) = +-0 ... exact
  74. * sqrt(inf) = inf
  75. * sqrt(-ve) = NaN ... with invalid signal
  76. * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
  77. */
  78. #include "libm.h"
  79. static const double tiny = 1.0e-300;
  80. double sqrt(double x)
  81. {
  82. double z;
  83. int32_t sign = (int)0x80000000;
  84. int32_t ix0,s0,q,m,t,i;
  85. uint32_t r,t1,s1,ix1,q1;
  86. EXTRACT_WORDS(ix0, ix1, x);
  87. /* take care of Inf and NaN */
  88. if ((ix0&0x7ff00000) == 0x7ff00000) {
  89. return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
  90. }
  91. /* take care of zero */
  92. if (ix0 <= 0) {
  93. if (((ix0&~sign)|ix1) == 0)
  94. return x; /* sqrt(+-0) = +-0 */
  95. if (ix0 < 0)
  96. return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
  97. }
  98. /* normalize x */
  99. m = ix0>>20;
  100. if (m == 0) { /* subnormal x */
  101. while (ix0 == 0) {
  102. m -= 21;
  103. ix0 |= (ix1>>11);
  104. ix1 <<= 21;
  105. }
  106. for (i=0; (ix0&0x00100000) == 0; i++)
  107. ix0<<=1;
  108. m -= i - 1;
  109. ix0 |= ix1>>(32-i);
  110. ix1 <<= i;
  111. }
  112. m -= 1023; /* unbias exponent */
  113. ix0 = (ix0&0x000fffff)|0x00100000;
  114. if (m & 1) { /* odd m, double x to make it even */
  115. ix0 += ix0 + ((ix1&sign)>>31);
  116. ix1 += ix1;
  117. }
  118. m >>= 1; /* m = [m/2] */
  119. /* generate sqrt(x) bit by bit */
  120. ix0 += ix0 + ((ix1&sign)>>31);
  121. ix1 += ix1;
  122. q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
  123. r = 0x00200000; /* r = moving bit from right to left */
  124. while (r != 0) {
  125. t = s0 + r;
  126. if (t <= ix0) {
  127. s0 = t + r;
  128. ix0 -= t;
  129. q += r;
  130. }
  131. ix0 += ix0 + ((ix1&sign)>>31);
  132. ix1 += ix1;
  133. r >>= 1;
  134. }
  135. r = sign;
  136. while (r != 0) {
  137. t1 = s1 + r;
  138. t = s0;
  139. if (t < ix0 || (t == ix0 && t1 <= ix1)) {
  140. s1 = t1 + r;
  141. if ((t1&sign) == sign && (s1&sign) == 0)
  142. s0++;
  143. ix0 -= t;
  144. if (ix1 < t1)
  145. ix0--;
  146. ix1 -= t1;
  147. q1 += r;
  148. }
  149. ix0 += ix0 + ((ix1&sign)>>31);
  150. ix1 += ix1;
  151. r >>= 1;
  152. }
  153. /* use floating add to find out rounding direction */
  154. if ((ix0|ix1) != 0) {
  155. z = 1.0 - tiny; /* raise inexact flag */
  156. if (z >= 1.0) {
  157. z = 1.0 + tiny;
  158. if (q1 == (uint32_t)0xffffffff) {
  159. q1 = 0;
  160. q++;
  161. } else if (z > 1.0) {
  162. if (q1 == (uint32_t)0xfffffffe)
  163. q++;
  164. q1 += 2;
  165. } else
  166. q1 += q1 & 1;
  167. }
  168. }
  169. ix0 = (q>>1) + 0x3fe00000;
  170. ix1 = q1>>1;
  171. if (q&1)
  172. ix1 |= sign;
  173. ix0 += m << 20;
  174. INSERT_WORDS(z, ix0, ix1);
  175. return z;
  176. }