sfsqrt.c 4.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174
  1. // SPDX-License-Identifier: GPL-2.0-or-later
  2. /*
  3. * Linux/PA-RISC Project (http://www.parisc-linux.org/)
  4. *
  5. * Floating-point emulation code
  6. * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <[email protected]>
  7. */
  8. /*
  9. * BEGIN_DESC
  10. *
  11. * File:
  12. * @(#) pa/spmath/sfsqrt.c $Revision: 1.1 $
  13. *
  14. * Purpose:
  15. * Single Floating-point Square Root
  16. *
  17. * External Interfaces:
  18. * sgl_fsqrt(srcptr,nullptr,dstptr,status)
  19. *
  20. * Internal Interfaces:
  21. *
  22. * Theory:
  23. * <<please update with a overview of the operation of this file>>
  24. *
  25. * END_DESC
  26. */
  27. #include "float.h"
  28. #include "sgl_float.h"
  29. /*
  30. * Single Floating-point Square Root
  31. */
  32. /*ARGSUSED*/
  33. unsigned int
  34. sgl_fsqrt(
  35. sgl_floating_point *srcptr,
  36. unsigned int *nullptr,
  37. sgl_floating_point *dstptr,
  38. unsigned int *status)
  39. {
  40. register unsigned int src, result;
  41. register int src_exponent;
  42. register unsigned int newbit, sum;
  43. register boolean guardbit = FALSE, even_exponent;
  44. src = *srcptr;
  45. /*
  46. * check source operand for NaN or infinity
  47. */
  48. if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
  49. /*
  50. * is signaling NaN?
  51. */
  52. if (Sgl_isone_signaling(src)) {
  53. /* trap if INVALIDTRAP enabled */
  54. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  55. /* make NaN quiet */
  56. Set_invalidflag();
  57. Sgl_set_quiet(src);
  58. }
  59. /*
  60. * Return quiet NaN or positive infinity.
  61. * Fall through to negative test if negative infinity.
  62. */
  63. if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
  64. *dstptr = src;
  65. return(NOEXCEPTION);
  66. }
  67. }
  68. /*
  69. * check for zero source operand
  70. */
  71. if (Sgl_iszero_exponentmantissa(src)) {
  72. *dstptr = src;
  73. return(NOEXCEPTION);
  74. }
  75. /*
  76. * check for negative source operand
  77. */
  78. if (Sgl_isone_sign(src)) {
  79. /* trap if INVALIDTRAP enabled */
  80. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  81. /* make NaN quiet */
  82. Set_invalidflag();
  83. Sgl_makequietnan(src);
  84. *dstptr = src;
  85. return(NOEXCEPTION);
  86. }
  87. /*
  88. * Generate result
  89. */
  90. if (src_exponent > 0) {
  91. even_exponent = Sgl_hidden(src);
  92. Sgl_clear_signexponent_set_hidden(src);
  93. }
  94. else {
  95. /* normalize operand */
  96. Sgl_clear_signexponent(src);
  97. src_exponent++;
  98. Sgl_normalize(src,src_exponent);
  99. even_exponent = src_exponent & 1;
  100. }
  101. if (even_exponent) {
  102. /* exponent is even */
  103. /* Add comment here. Explain why odd exponent needs correction */
  104. Sgl_leftshiftby1(src);
  105. }
  106. /*
  107. * Add comment here. Explain following algorithm.
  108. *
  109. * Trust me, it works.
  110. *
  111. */
  112. Sgl_setzero(result);
  113. newbit = 1 << SGL_P;
  114. while (newbit && Sgl_isnotzero(src)) {
  115. Sgl_addition(result,newbit,sum);
  116. if(sum <= Sgl_all(src)) {
  117. /* update result */
  118. Sgl_addition(result,(newbit<<1),result);
  119. Sgl_subtract(src,sum,src);
  120. }
  121. Sgl_rightshiftby1(newbit);
  122. Sgl_leftshiftby1(src);
  123. }
  124. /* correct exponent for pre-shift */
  125. if (even_exponent) {
  126. Sgl_rightshiftby1(result);
  127. }
  128. /* check for inexact */
  129. if (Sgl_isnotzero(src)) {
  130. if (!even_exponent && Sgl_islessthan(result,src))
  131. Sgl_increment(result);
  132. guardbit = Sgl_lowmantissa(result);
  133. Sgl_rightshiftby1(result);
  134. /* now round result */
  135. switch (Rounding_mode()) {
  136. case ROUNDPLUS:
  137. Sgl_increment(result);
  138. break;
  139. case ROUNDNEAREST:
  140. /* stickybit is always true, so guardbit
  141. * is enough to determine rounding */
  142. if (guardbit) {
  143. Sgl_increment(result);
  144. }
  145. break;
  146. }
  147. /* increment result exponent by 1 if mantissa overflowed */
  148. if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
  149. if (Is_inexacttrap_enabled()) {
  150. Sgl_set_exponent(result,
  151. ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
  152. *dstptr = result;
  153. return(INEXACTEXCEPTION);
  154. }
  155. else Set_inexactflag();
  156. }
  157. else {
  158. Sgl_rightshiftby1(result);
  159. }
  160. Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
  161. *dstptr = result;
  162. return(NOEXCEPTION);
  163. }