dp_sqrt.c 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  1. /* IEEE754 floating point arithmetic
  2. * double precision square root
  3. */
  4. /*
  5. * MIPS floating point support
  6. * Copyright (C) 1994-2000 Algorithmics Ltd.
  7. *
  8. * This program is free software; you can distribute it and/or modify it
  9. * under the terms of the GNU General Public License (Version 2) as
  10. * published by the Free Software Foundation.
  11. *
  12. * This program is distributed in the hope it will be useful, but WITHOUT
  13. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  15. * for more details.
  16. *
  17. * You should have received a copy of the GNU General Public License along
  18. * with this program; if not, write to the Free Software Foundation, Inc.,
  19. * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  20. */
  21. #include "ieee754dp.h"
  22. static const unsigned table[] = {
  23. 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592,
  24. 29598, 36145, 43202, 50740, 58733, 67158, 75992,
  25. 85215, 83599, 71378, 60428, 50647, 41945, 34246,
  26. 27478, 21581, 16499, 12183, 8588, 5674, 3403,
  27. 1742, 661, 130
  28. };
  29. union ieee754dp ieee754dp_sqrt(union ieee754dp x)
  30. {
  31. struct _ieee754_csr oldcsr;
  32. union ieee754dp y, z, t;
  33. unsigned scalx, yh;
  34. COMPXDP;
  35. EXPLODEXDP;
  36. ieee754_clearcx();
  37. FLUSHXDP;
  38. /* x == INF or NAN? */
  39. switch (xc) {
  40. case IEEE754_CLASS_SNAN:
  41. return ieee754dp_nanxcpt(x);
  42. case IEEE754_CLASS_QNAN:
  43. /* sqrt(Nan) = Nan */
  44. return x;
  45. case IEEE754_CLASS_ZERO:
  46. /* sqrt(0) = 0 */
  47. return x;
  48. case IEEE754_CLASS_INF:
  49. if (xs) {
  50. /* sqrt(-Inf) = Nan */
  51. ieee754_setcx(IEEE754_INVALID_OPERATION);
  52. return ieee754dp_indef();
  53. }
  54. /* sqrt(+Inf) = Inf */
  55. return x;
  56. case IEEE754_CLASS_DNORM:
  57. DPDNORMX;
  58. /* fall through */
  59. case IEEE754_CLASS_NORM:
  60. if (xs) {
  61. /* sqrt(-x) = Nan */
  62. ieee754_setcx(IEEE754_INVALID_OPERATION);
  63. return ieee754dp_indef();
  64. }
  65. break;
  66. }
  67. /* save old csr; switch off INX enable & flag; set RN rounding */
  68. oldcsr = ieee754_csr;
  69. ieee754_csr.mx &= ~IEEE754_INEXACT;
  70. ieee754_csr.sx &= ~IEEE754_INEXACT;
  71. ieee754_csr.rm = FPU_CSR_RN;
  72. /* adjust exponent to prevent overflow */
  73. scalx = 0;
  74. if (xe > 512) { /* x > 2**-512? */
  75. xe -= 512; /* x = x / 2**512 */
  76. scalx += 256;
  77. } else if (xe < -512) { /* x < 2**-512? */
  78. xe += 512; /* x = x * 2**512 */
  79. scalx -= 256;
  80. }
  81. y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
  82. /* magic initial approximation to almost 8 sig. bits */
  83. yh = y.bits >> 32;
  84. yh = (yh >> 1) + 0x1ff80000;
  85. yh = yh - table[(yh >> 15) & 31];
  86. y.bits = ((u64) yh << 32) | (y.bits & 0xffffffff);
  87. /* Heron's rule once with correction to improve to ~18 sig. bits */
  88. /* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */
  89. t = ieee754dp_div(x, y);
  90. y = ieee754dp_add(y, t);
  91. y.bits -= 0x0010000600000000LL;
  92. y.bits &= 0xffffffff00000000LL;
  93. /* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */
  94. /* t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */
  95. z = t = ieee754dp_mul(y, y);
  96. t.bexp += 0x001;
  97. t = ieee754dp_add(t, z);
  98. z = ieee754dp_mul(ieee754dp_sub(x, z), y);
  99. /* t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; */
  100. t = ieee754dp_div(z, ieee754dp_add(t, x));
  101. t.bexp += 0x001;
  102. y = ieee754dp_add(y, t);
  103. /* twiddle last bit to force y correctly rounded */
  104. /* set RZ, clear INEX flag */
  105. ieee754_csr.rm = FPU_CSR_RZ;
  106. ieee754_csr.sx &= ~IEEE754_INEXACT;
  107. /* t=x/y; ...chopped quotient, possibly inexact */
  108. t = ieee754dp_div(x, y);
  109. if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) {
  110. if (!(ieee754_csr.sx & IEEE754_INEXACT))
  111. /* t = t-ulp */
  112. t.bits -= 1;
  113. /* add inexact to result status */
  114. oldcsr.cx |= IEEE754_INEXACT;
  115. oldcsr.sx |= IEEE754_INEXACT;
  116. switch (oldcsr.rm) {
  117. case FPU_CSR_RU:
  118. y.bits += 1;
  119. /* drop through */
  120. case FPU_CSR_RN:
  121. t.bits += 1;
  122. break;
  123. }
  124. /* y=y+t; ...chopped sum */
  125. y = ieee754dp_add(y, t);
  126. /* adjust scalx for correctly rounded sqrt(x) */
  127. scalx -= 1;
  128. }
  129. /* py[n0]=py[n0]+scalx; ...scale back y */
  130. y.bexp += scalx;
  131. /* restore rounding mode, possibly set inexact */
  132. ieee754_csr = oldcsr;
  133. return y;
  134. }