sfrem.c 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277
  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) <bame@debian.org>
  7. */
  8. /*
  9. * BEGIN_DESC
  10. *
  11. * File:
  12. * @(#) pa/spmath/sfrem.c $Revision: 1.1 $
  13. *
  14. * Purpose:
  15. * Single Precision Floating-point Remainder
  16. *
  17. * External Interfaces:
  18. * sgl_frem(srcptr1,srcptr2,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 Precision Floating-point Remainder
  31. */
  32. int
  33. sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
  34. sgl_floating_point * dstptr, unsigned int *status)
  35. {
  36. register unsigned int opnd1, opnd2, result;
  37. register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
  38. register boolean roundup = FALSE;
  39. opnd1 = *srcptr1;
  40. opnd2 = *srcptr2;
  41. /*
  42. * check first operand for NaN's or infinity
  43. */
  44. if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
  45. if (Sgl_iszero_mantissa(opnd1)) {
  46. if (Sgl_isnotnan(opnd2)) {
  47. /* invalid since first operand is infinity */
  48. if (Is_invalidtrap_enabled())
  49. return(INVALIDEXCEPTION);
  50. Set_invalidflag();
  51. Sgl_makequietnan(result);
  52. *dstptr = result;
  53. return(NOEXCEPTION);
  54. }
  55. }
  56. else {
  57. /*
  58. * is NaN; signaling or quiet?
  59. */
  60. if (Sgl_isone_signaling(opnd1)) {
  61. /* trap if INVALIDTRAP enabled */
  62. if (Is_invalidtrap_enabled())
  63. return(INVALIDEXCEPTION);
  64. /* make NaN quiet */
  65. Set_invalidflag();
  66. Sgl_set_quiet(opnd1);
  67. }
  68. /*
  69. * is second operand a signaling NaN?
  70. */
  71. else if (Sgl_is_signalingnan(opnd2)) {
  72. /* trap if INVALIDTRAP enabled */
  73. if (Is_invalidtrap_enabled())
  74. return(INVALIDEXCEPTION);
  75. /* make NaN quiet */
  76. Set_invalidflag();
  77. Sgl_set_quiet(opnd2);
  78. *dstptr = opnd2;
  79. return(NOEXCEPTION);
  80. }
  81. /*
  82. * return quiet NaN
  83. */
  84. *dstptr = opnd1;
  85. return(NOEXCEPTION);
  86. }
  87. }
  88. /*
  89. * check second operand for NaN's or infinity
  90. */
  91. if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
  92. if (Sgl_iszero_mantissa(opnd2)) {
  93. /*
  94. * return first operand
  95. */
  96. *dstptr = opnd1;
  97. return(NOEXCEPTION);
  98. }
  99. /*
  100. * is NaN; signaling or quiet?
  101. */
  102. if (Sgl_isone_signaling(opnd2)) {
  103. /* trap if INVALIDTRAP enabled */
  104. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  105. /* make NaN quiet */
  106. Set_invalidflag();
  107. Sgl_set_quiet(opnd2);
  108. }
  109. /*
  110. * return quiet NaN
  111. */
  112. *dstptr = opnd2;
  113. return(NOEXCEPTION);
  114. }
  115. /*
  116. * check second operand for zero
  117. */
  118. if (Sgl_iszero_exponentmantissa(opnd2)) {
  119. /* invalid since second operand is zero */
  120. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  121. Set_invalidflag();
  122. Sgl_makequietnan(result);
  123. *dstptr = result;
  124. return(NOEXCEPTION);
  125. }
  126. /*
  127. * get sign of result
  128. */
  129. result = opnd1;
  130. /*
  131. * check for denormalized operands
  132. */
  133. if (opnd1_exponent == 0) {
  134. /* check for zero */
  135. if (Sgl_iszero_mantissa(opnd1)) {
  136. *dstptr = opnd1;
  137. return(NOEXCEPTION);
  138. }
  139. /* normalize, then continue */
  140. opnd1_exponent = 1;
  141. Sgl_normalize(opnd1,opnd1_exponent);
  142. }
  143. else {
  144. Sgl_clear_signexponent_set_hidden(opnd1);
  145. }
  146. if (opnd2_exponent == 0) {
  147. /* normalize, then continue */
  148. opnd2_exponent = 1;
  149. Sgl_normalize(opnd2,opnd2_exponent);
  150. }
  151. else {
  152. Sgl_clear_signexponent_set_hidden(opnd2);
  153. }
  154. /* find result exponent and divide step loop count */
  155. dest_exponent = opnd2_exponent - 1;
  156. stepcount = opnd1_exponent - opnd2_exponent;
  157. /*
  158. * check for opnd1/opnd2 < 1
  159. */
  160. if (stepcount < 0) {
  161. /*
  162. * check for opnd1/opnd2 > 1/2
  163. *
  164. * In this case n will round to 1, so
  165. * r = opnd1 - opnd2
  166. */
  167. if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
  168. Sgl_all(result) = ~Sgl_all(result); /* set sign */
  169. /* align opnd2 with opnd1 */
  170. Sgl_leftshiftby1(opnd2);
  171. Sgl_subtract(opnd2,opnd1,opnd2);
  172. /* now normalize */
  173. while (Sgl_iszero_hidden(opnd2)) {
  174. Sgl_leftshiftby1(opnd2);
  175. dest_exponent--;
  176. }
  177. Sgl_set_exponentmantissa(result,opnd2);
  178. goto testforunderflow;
  179. }
  180. /*
  181. * opnd1/opnd2 <= 1/2
  182. *
  183. * In this case n will round to zero, so
  184. * r = opnd1
  185. */
  186. Sgl_set_exponentmantissa(result,opnd1);
  187. dest_exponent = opnd1_exponent;
  188. goto testforunderflow;
  189. }
  190. /*
  191. * Generate result
  192. *
  193. * Do iterative subtract until remainder is less than operand 2.
  194. */
  195. while (stepcount-- > 0 && Sgl_all(opnd1)) {
  196. if (Sgl_isnotlessthan(opnd1,opnd2))
  197. Sgl_subtract(opnd1,opnd2,opnd1);
  198. Sgl_leftshiftby1(opnd1);
  199. }
  200. /*
  201. * Do last subtract, then determine which way to round if remainder
  202. * is exactly 1/2 of opnd2
  203. */
  204. if (Sgl_isnotlessthan(opnd1,opnd2)) {
  205. Sgl_subtract(opnd1,opnd2,opnd1);
  206. roundup = TRUE;
  207. }
  208. if (stepcount > 0 || Sgl_iszero(opnd1)) {
  209. /* division is exact, remainder is zero */
  210. Sgl_setzero_exponentmantissa(result);
  211. *dstptr = result;
  212. return(NOEXCEPTION);
  213. }
  214. /*
  215. * Check for cases where opnd1/opnd2 < n
  216. *
  217. * In this case the result's sign will be opposite that of
  218. * opnd1. The mantissa also needs some correction.
  219. */
  220. Sgl_leftshiftby1(opnd1);
  221. if (Sgl_isgreaterthan(opnd1,opnd2)) {
  222. Sgl_invert_sign(result);
  223. Sgl_subtract((opnd2<<1),opnd1,opnd1);
  224. }
  225. /* check for remainder being exactly 1/2 of opnd2 */
  226. else if (Sgl_isequal(opnd1,opnd2) && roundup) {
  227. Sgl_invert_sign(result);
  228. }
  229. /* normalize result's mantissa */
  230. while (Sgl_iszero_hidden(opnd1)) {
  231. dest_exponent--;
  232. Sgl_leftshiftby1(opnd1);
  233. }
  234. Sgl_set_exponentmantissa(result,opnd1);
  235. /*
  236. * Test for underflow
  237. */
  238. testforunderflow:
  239. if (dest_exponent <= 0) {
  240. /* trap if UNDERFLOWTRAP enabled */
  241. if (Is_underflowtrap_enabled()) {
  242. /*
  243. * Adjust bias of result
  244. */
  245. Sgl_setwrapped_exponent(result,dest_exponent,unfl);
  246. *dstptr = result;
  247. /* frem is always exact */
  248. return(UNDERFLOWEXCEPTION);
  249. }
  250. /*
  251. * denormalize result or set to signed zero
  252. */
  253. if (dest_exponent >= (1 - SGL_P)) {
  254. Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
  255. }
  256. else {
  257. Sgl_setzero_exponentmantissa(result);
  258. }
  259. }
  260. else Sgl_set_exponent(result,dest_exponent);
  261. *dstptr = result;
  262. return(NOEXCEPTION);
  263. }