sp_sqrt.c 2.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114
  1. /* IEEE754 floating point arithmetic
  2. * single 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 "ieee754sp.h"
  22. union ieee754sp ieee754sp_sqrt(union ieee754sp x)
  23. {
  24. int ix, s, q, m, t, i;
  25. unsigned int r;
  26. COMPXSP;
  27. /* take care of Inf and NaN */
  28. EXPLODEXSP;
  29. ieee754_clearcx();
  30. FLUSHXSP;
  31. /* x == INF or NAN? */
  32. switch (xc) {
  33. case IEEE754_CLASS_SNAN:
  34. return ieee754sp_nanxcpt(x);
  35. case IEEE754_CLASS_QNAN:
  36. /* sqrt(Nan) = Nan */
  37. return x;
  38. case IEEE754_CLASS_ZERO:
  39. /* sqrt(0) = 0 */
  40. return x;
  41. case IEEE754_CLASS_INF:
  42. if (xs) {
  43. /* sqrt(-Inf) = Nan */
  44. ieee754_setcx(IEEE754_INVALID_OPERATION);
  45. return ieee754sp_indef();
  46. }
  47. /* sqrt(+Inf) = Inf */
  48. return x;
  49. case IEEE754_CLASS_DNORM:
  50. case IEEE754_CLASS_NORM:
  51. if (xs) {
  52. /* sqrt(-x) = Nan */
  53. ieee754_setcx(IEEE754_INVALID_OPERATION);
  54. return ieee754sp_indef();
  55. }
  56. break;
  57. }
  58. ix = x.bits;
  59. /* normalize x */
  60. m = (ix >> 23);
  61. if (m == 0) { /* subnormal x */
  62. for (i = 0; (ix & 0x00800000) == 0; i++)
  63. ix <<= 1;
  64. m -= i - 1;
  65. }
  66. m -= 127; /* unbias exponent */
  67. ix = (ix & 0x007fffff) | 0x00800000;
  68. if (m & 1) /* odd m, double x to make it even */
  69. ix += ix;
  70. m >>= 1; /* m = [m/2] */
  71. /* generate sqrt(x) bit by bit */
  72. ix += ix;
  73. q = s = 0; /* q = sqrt(x) */
  74. r = 0x01000000; /* r = moving bit from right to left */
  75. while (r != 0) {
  76. t = s + r;
  77. if (t <= ix) {
  78. s = t + r;
  79. ix -= t;
  80. q += r;
  81. }
  82. ix += ix;
  83. r >>= 1;
  84. }
  85. if (ix != 0) {
  86. ieee754_setcx(IEEE754_INEXACT);
  87. switch (ieee754_csr.rm) {
  88. case FPU_CSR_RU:
  89. q += 2;
  90. break;
  91. case FPU_CSR_RN:
  92. q += (q & 1);
  93. break;
  94. }
  95. }
  96. ix = (q >> 1) + 0x3f000000;
  97. ix += (m << 23);
  98. x.bits = ix;
  99. return x;
  100. }