sp_sqrt.c 1.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103
  1. // SPDX-License-Identifier: GPL-2.0-only
  2. /* IEEE754 floating point arithmetic
  3. * single precision square root
  4. */
  5. /*
  6. * MIPS floating point support
  7. * Copyright (C) 1994-2000 Algorithmics Ltd.
  8. */
  9. #include "ieee754sp.h"
  10. union ieee754sp ieee754sp_sqrt(union ieee754sp x)
  11. {
  12. int ix, s, q, m, t, i;
  13. unsigned int r;
  14. COMPXSP;
  15. /* take care of Inf and NaN */
  16. EXPLODEXSP;
  17. ieee754_clearcx();
  18. FLUSHXSP;
  19. /* x == INF or NAN? */
  20. switch (xc) {
  21. case IEEE754_CLASS_SNAN:
  22. return ieee754sp_nanxcpt(x);
  23. case IEEE754_CLASS_QNAN:
  24. /* sqrt(Nan) = Nan */
  25. return x;
  26. case IEEE754_CLASS_ZERO:
  27. /* sqrt(0) = 0 */
  28. return x;
  29. case IEEE754_CLASS_INF:
  30. if (xs) {
  31. /* sqrt(-Inf) = Nan */
  32. ieee754_setcx(IEEE754_INVALID_OPERATION);
  33. return ieee754sp_indef();
  34. }
  35. /* sqrt(+Inf) = Inf */
  36. return x;
  37. case IEEE754_CLASS_DNORM:
  38. case IEEE754_CLASS_NORM:
  39. if (xs) {
  40. /* sqrt(-x) = Nan */
  41. ieee754_setcx(IEEE754_INVALID_OPERATION);
  42. return ieee754sp_indef();
  43. }
  44. break;
  45. }
  46. ix = x.bits;
  47. /* normalize x */
  48. m = (ix >> 23);
  49. if (m == 0) { /* subnormal x */
  50. for (i = 0; (ix & 0x00800000) == 0; i++)
  51. ix <<= 1;
  52. m -= i - 1;
  53. }
  54. m -= 127; /* unbias exponent */
  55. ix = (ix & 0x007fffff) | 0x00800000;
  56. if (m & 1) /* odd m, double x to make it even */
  57. ix += ix;
  58. m >>= 1; /* m = [m/2] */
  59. /* generate sqrt(x) bit by bit */
  60. ix += ix;
  61. s = 0;
  62. q = 0; /* q = sqrt(x) */
  63. r = 0x01000000; /* r = moving bit from right to left */
  64. while (r != 0) {
  65. t = s + r;
  66. if (t <= ix) {
  67. s = t + r;
  68. ix -= t;
  69. q += r;
  70. }
  71. ix += ix;
  72. r >>= 1;
  73. }
  74. if (ix != 0) {
  75. ieee754_setcx(IEEE754_INEXACT);
  76. switch (ieee754_csr.rm) {
  77. case FPU_CSR_RU:
  78. q += 2;
  79. break;
  80. case FPU_CSR_RN:
  81. q += (q & 1);
  82. break;
  83. }
  84. }
  85. ix = (q >> 1) + 0x3f000000;
  86. ix += (m << 23);
  87. x.bits = ix;
  88. return x;
  89. }