dfrem.c 8.1 KB

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