dfdiv.c 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387
  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/dfdiv.c $Revision: 1.1 $
  13. *
  14. * Purpose:
  15. * Double Precision Floating-point Divide
  16. *
  17. * External Interfaces:
  18. * dbl_fdiv(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 Divide
  31. */
  32. int
  33. dbl_fdiv (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 opnd3p1, opnd3p2, resultp1, resultp2;
  38. register int dest_exponent, count;
  39. register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  40. boolean is_tiny;
  41. Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
  42. Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
  43. /*
  44. * set sign bit of result
  45. */
  46. if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
  47. Dbl_setnegativezerop1(resultp1);
  48. else Dbl_setzerop1(resultp1);
  49. /*
  50. * check first operand for NaN's or infinity
  51. */
  52. if (Dbl_isinfinity_exponent(opnd1p1)) {
  53. if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  54. if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
  55. if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
  56. /*
  57. * invalid since both operands
  58. * are infinity
  59. */
  60. if (Is_invalidtrap_enabled())
  61. return(INVALIDEXCEPTION);
  62. Set_invalidflag();
  63. Dbl_makequietnan(resultp1,resultp2);
  64. Dbl_copytoptr(resultp1,resultp2,dstptr);
  65. return(NOEXCEPTION);
  66. }
  67. /*
  68. * return infinity
  69. */
  70. Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
  71. Dbl_copytoptr(resultp1,resultp2,dstptr);
  72. return(NOEXCEPTION);
  73. }
  74. }
  75. else {
  76. /*
  77. * is NaN; signaling or quiet?
  78. */
  79. if (Dbl_isone_signaling(opnd1p1)) {
  80. /* trap if INVALIDTRAP enabled */
  81. if (Is_invalidtrap_enabled())
  82. return(INVALIDEXCEPTION);
  83. /* make NaN quiet */
  84. Set_invalidflag();
  85. Dbl_set_quiet(opnd1p1);
  86. }
  87. /*
  88. * is second operand a signaling NaN?
  89. */
  90. else if (Dbl_is_signalingnan(opnd2p1)) {
  91. /* trap if INVALIDTRAP enabled */
  92. if (Is_invalidtrap_enabled())
  93. return(INVALIDEXCEPTION);
  94. /* make NaN quiet */
  95. Set_invalidflag();
  96. Dbl_set_quiet(opnd2p1);
  97. Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
  98. return(NOEXCEPTION);
  99. }
  100. /*
  101. * return quiet NaN
  102. */
  103. Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
  104. return(NOEXCEPTION);
  105. }
  106. }
  107. /*
  108. * check second operand for NaN's or infinity
  109. */
  110. if (Dbl_isinfinity_exponent(opnd2p1)) {
  111. if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
  112. /*
  113. * return zero
  114. */
  115. Dbl_setzero_exponentmantissa(resultp1,resultp2);
  116. Dbl_copytoptr(resultp1,resultp2,dstptr);
  117. return(NOEXCEPTION);
  118. }
  119. /*
  120. * is NaN; signaling or quiet?
  121. */
  122. if (Dbl_isone_signaling(opnd2p1)) {
  123. /* trap if INVALIDTRAP enabled */
  124. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  125. /* make NaN quiet */
  126. Set_invalidflag();
  127. Dbl_set_quiet(opnd2p1);
  128. }
  129. /*
  130. * return quiet NaN
  131. */
  132. Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
  133. return(NOEXCEPTION);
  134. }
  135. /*
  136. * check for division by zero
  137. */
  138. if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
  139. if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
  140. /* invalid since both operands are zero */
  141. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  142. Set_invalidflag();
  143. Dbl_makequietnan(resultp1,resultp2);
  144. Dbl_copytoptr(resultp1,resultp2,dstptr);
  145. return(NOEXCEPTION);
  146. }
  147. if (Is_divisionbyzerotrap_enabled())
  148. return(DIVISIONBYZEROEXCEPTION);
  149. Set_divisionbyzeroflag();
  150. Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
  151. Dbl_copytoptr(resultp1,resultp2,dstptr);
  152. return(NOEXCEPTION);
  153. }
  154. /*
  155. * Generate exponent
  156. */
  157. dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
  158. /*
  159. * Generate mantissa
  160. */
  161. if (Dbl_isnotzero_exponent(opnd1p1)) {
  162. /* set hidden bit */
  163. Dbl_clear_signexponent_set_hidden(opnd1p1);
  164. }
  165. else {
  166. /* check for zero */
  167. if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  168. Dbl_setzero_exponentmantissa(resultp1,resultp2);
  169. Dbl_copytoptr(resultp1,resultp2,dstptr);
  170. return(NOEXCEPTION);
  171. }
  172. /* is denormalized, want to normalize */
  173. Dbl_clear_signexponent(opnd1p1);
  174. Dbl_leftshiftby1(opnd1p1,opnd1p2);
  175. Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
  176. }
  177. /* opnd2 needs to have hidden bit set with msb in hidden bit */
  178. if (Dbl_isnotzero_exponent(opnd2p1)) {
  179. Dbl_clear_signexponent_set_hidden(opnd2p1);
  180. }
  181. else {
  182. /* is denormalized; want to normalize */
  183. Dbl_clear_signexponent(opnd2p1);
  184. Dbl_leftshiftby1(opnd2p1,opnd2p2);
  185. while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
  186. dest_exponent+=8;
  187. Dbl_leftshiftby8(opnd2p1,opnd2p2);
  188. }
  189. if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
  190. dest_exponent+=4;
  191. Dbl_leftshiftby4(opnd2p1,opnd2p2);
  192. }
  193. while (Dbl_iszero_hidden(opnd2p1)) {
  194. dest_exponent++;
  195. Dbl_leftshiftby1(opnd2p1,opnd2p2);
  196. }
  197. }
  198. /* Divide the source mantissas */
  199. /*
  200. * A non-restoring divide algorithm is used.
  201. */
  202. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  203. Dbl_setzero(opnd3p1,opnd3p2);
  204. for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
  205. Dbl_leftshiftby1(opnd1p1,opnd1p2);
  206. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  207. if (Dbl_iszero_sign(opnd1p1)) {
  208. Dbl_setone_lowmantissap2(opnd3p2);
  209. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  210. }
  211. else {
  212. Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
  213. }
  214. }
  215. if (count <= DBL_P) {
  216. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  217. Dbl_setone_lowmantissap2(opnd3p2);
  218. Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
  219. if (Dbl_iszero_hidden(opnd3p1)) {
  220. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  221. dest_exponent--;
  222. }
  223. }
  224. else {
  225. if (Dbl_iszero_hidden(opnd3p1)) {
  226. /* need to get one more bit of result */
  227. Dbl_leftshiftby1(opnd1p1,opnd1p2);
  228. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  229. if (Dbl_iszero_sign(opnd1p1)) {
  230. Dbl_setone_lowmantissap2(opnd3p2);
  231. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  232. }
  233. else {
  234. Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  235. }
  236. dest_exponent--;
  237. }
  238. if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
  239. stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
  240. }
  241. inexact = guardbit | stickybit;
  242. /*
  243. * round result
  244. */
  245. if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
  246. Dbl_clear_signexponent(opnd3p1);
  247. switch (Rounding_mode()) {
  248. case ROUNDPLUS:
  249. if (Dbl_iszero_sign(resultp1))
  250. Dbl_increment(opnd3p1,opnd3p2);
  251. break;
  252. case ROUNDMINUS:
  253. if (Dbl_isone_sign(resultp1))
  254. Dbl_increment(opnd3p1,opnd3p2);
  255. break;
  256. case ROUNDNEAREST:
  257. if (guardbit && (stickybit ||
  258. Dbl_isone_lowmantissap2(opnd3p2))) {
  259. Dbl_increment(opnd3p1,opnd3p2);
  260. }
  261. }
  262. if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
  263. }
  264. Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
  265. /*
  266. * Test for overflow
  267. */
  268. if (dest_exponent >= DBL_INFINITY_EXPONENT) {
  269. /* trap if OVERFLOWTRAP enabled */
  270. if (Is_overflowtrap_enabled()) {
  271. /*
  272. * Adjust bias of result
  273. */
  274. Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
  275. Dbl_copytoptr(resultp1,resultp2,dstptr);
  276. if (inexact)
  277. if (Is_inexacttrap_enabled())
  278. return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  279. else Set_inexactflag();
  280. return(OVERFLOWEXCEPTION);
  281. }
  282. Set_overflowflag();
  283. /* set result to infinity or largest number */
  284. Dbl_setoverflow(resultp1,resultp2);
  285. inexact = TRUE;
  286. }
  287. /*
  288. * Test for underflow
  289. */
  290. else if (dest_exponent <= 0) {
  291. /* trap if UNDERFLOWTRAP enabled */
  292. if (Is_underflowtrap_enabled()) {
  293. /*
  294. * Adjust bias of result
  295. */
  296. Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
  297. Dbl_copytoptr(resultp1,resultp2,dstptr);
  298. if (inexact)
  299. if (Is_inexacttrap_enabled())
  300. return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
  301. else Set_inexactflag();
  302. return(UNDERFLOWEXCEPTION);
  303. }
  304. /* Determine if should set underflow flag */
  305. is_tiny = TRUE;
  306. if (dest_exponent == 0 && inexact) {
  307. switch (Rounding_mode()) {
  308. case ROUNDPLUS:
  309. if (Dbl_iszero_sign(resultp1)) {
  310. Dbl_increment(opnd3p1,opnd3p2);
  311. if (Dbl_isone_hiddenoverflow(opnd3p1))
  312. is_tiny = FALSE;
  313. Dbl_decrement(opnd3p1,opnd3p2);
  314. }
  315. break;
  316. case ROUNDMINUS:
  317. if (Dbl_isone_sign(resultp1)) {
  318. Dbl_increment(opnd3p1,opnd3p2);
  319. if (Dbl_isone_hiddenoverflow(opnd3p1))
  320. is_tiny = FALSE;
  321. Dbl_decrement(opnd3p1,opnd3p2);
  322. }
  323. break;
  324. case ROUNDNEAREST:
  325. if (guardbit && (stickybit ||
  326. Dbl_isone_lowmantissap2(opnd3p2))) {
  327. Dbl_increment(opnd3p1,opnd3p2);
  328. if (Dbl_isone_hiddenoverflow(opnd3p1))
  329. is_tiny = FALSE;
  330. Dbl_decrement(opnd3p1,opnd3p2);
  331. }
  332. break;
  333. }
  334. }
  335. /*
  336. * denormalize result or set to signed zero
  337. */
  338. stickybit = inexact;
  339. Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
  340. stickybit,inexact);
  341. /* return rounded number */
  342. if (inexact) {
  343. switch (Rounding_mode()) {
  344. case ROUNDPLUS:
  345. if (Dbl_iszero_sign(resultp1)) {
  346. Dbl_increment(opnd3p1,opnd3p2);
  347. }
  348. break;
  349. case ROUNDMINUS:
  350. if (Dbl_isone_sign(resultp1)) {
  351. Dbl_increment(opnd3p1,opnd3p2);
  352. }
  353. break;
  354. case ROUNDNEAREST:
  355. if (guardbit && (stickybit ||
  356. Dbl_isone_lowmantissap2(opnd3p2))) {
  357. Dbl_increment(opnd3p1,opnd3p2);
  358. }
  359. break;
  360. }
  361. if (is_tiny) Set_underflowflag();
  362. }
  363. Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
  364. }
  365. else Dbl_set_exponent(resultp1,dest_exponent);
  366. Dbl_copytoptr(resultp1,resultp2,dstptr);
  367. /* check for inexact */
  368. if (inexact) {
  369. if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  370. else Set_inexactflag();
  371. }
  372. return(NOEXCEPTION);
  373. }