sfdiv.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379
  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/sfdiv.c $Revision: 1.1 $
  13. *
  14. * Purpose:
  15. * Single Precision Floating-point Divide
  16. *
  17. * External Interfaces:
  18. * sgl_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 "sgl_float.h"
  29. /*
  30. * Single Precision Floating-point Divide
  31. */
  32. int
  33. sgl_fdiv (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
  34. sgl_floating_point * dstptr, unsigned int *status)
  35. {
  36. register unsigned int opnd1, opnd2, opnd3, result;
  37. register int dest_exponent, count;
  38. register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  39. boolean is_tiny;
  40. opnd1 = *srcptr1;
  41. opnd2 = *srcptr2;
  42. /*
  43. * set sign bit of result
  44. */
  45. if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
  46. else Sgl_setzero(result);
  47. /*
  48. * check first operand for NaN's or infinity
  49. */
  50. if (Sgl_isinfinity_exponent(opnd1)) {
  51. if (Sgl_iszero_mantissa(opnd1)) {
  52. if (Sgl_isnotnan(opnd2)) {
  53. if (Sgl_isinfinity(opnd2)) {
  54. /*
  55. * invalid since both operands
  56. * are infinity
  57. */
  58. if (Is_invalidtrap_enabled())
  59. return(INVALIDEXCEPTION);
  60. Set_invalidflag();
  61. Sgl_makequietnan(result);
  62. *dstptr = result;
  63. return(NOEXCEPTION);
  64. }
  65. /*
  66. * return infinity
  67. */
  68. Sgl_setinfinity_exponentmantissa(result);
  69. *dstptr = result;
  70. return(NOEXCEPTION);
  71. }
  72. }
  73. else {
  74. /*
  75. * is NaN; signaling or quiet?
  76. */
  77. if (Sgl_isone_signaling(opnd1)) {
  78. /* trap if INVALIDTRAP enabled */
  79. if (Is_invalidtrap_enabled())
  80. return(INVALIDEXCEPTION);
  81. /* make NaN quiet */
  82. Set_invalidflag();
  83. Sgl_set_quiet(opnd1);
  84. }
  85. /*
  86. * is second operand a signaling NaN?
  87. */
  88. else if (Sgl_is_signalingnan(opnd2)) {
  89. /* trap if INVALIDTRAP enabled */
  90. if (Is_invalidtrap_enabled())
  91. return(INVALIDEXCEPTION);
  92. /* make NaN quiet */
  93. Set_invalidflag();
  94. Sgl_set_quiet(opnd2);
  95. *dstptr = opnd2;
  96. return(NOEXCEPTION);
  97. }
  98. /*
  99. * return quiet NaN
  100. */
  101. *dstptr = opnd1;
  102. return(NOEXCEPTION);
  103. }
  104. }
  105. /*
  106. * check second operand for NaN's or infinity
  107. */
  108. if (Sgl_isinfinity_exponent(opnd2)) {
  109. if (Sgl_iszero_mantissa(opnd2)) {
  110. /*
  111. * return zero
  112. */
  113. Sgl_setzero_exponentmantissa(result);
  114. *dstptr = result;
  115. return(NOEXCEPTION);
  116. }
  117. /*
  118. * is NaN; signaling or quiet?
  119. */
  120. if (Sgl_isone_signaling(opnd2)) {
  121. /* trap if INVALIDTRAP enabled */
  122. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  123. /* make NaN quiet */
  124. Set_invalidflag();
  125. Sgl_set_quiet(opnd2);
  126. }
  127. /*
  128. * return quiet NaN
  129. */
  130. *dstptr = opnd2;
  131. return(NOEXCEPTION);
  132. }
  133. /*
  134. * check for division by zero
  135. */
  136. if (Sgl_iszero_exponentmantissa(opnd2)) {
  137. if (Sgl_iszero_exponentmantissa(opnd1)) {
  138. /* invalid since both operands are zero */
  139. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  140. Set_invalidflag();
  141. Sgl_makequietnan(result);
  142. *dstptr = result;
  143. return(NOEXCEPTION);
  144. }
  145. if (Is_divisionbyzerotrap_enabled())
  146. return(DIVISIONBYZEROEXCEPTION);
  147. Set_divisionbyzeroflag();
  148. Sgl_setinfinity_exponentmantissa(result);
  149. *dstptr = result;
  150. return(NOEXCEPTION);
  151. }
  152. /*
  153. * Generate exponent
  154. */
  155. dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
  156. /*
  157. * Generate mantissa
  158. */
  159. if (Sgl_isnotzero_exponent(opnd1)) {
  160. /* set hidden bit */
  161. Sgl_clear_signexponent_set_hidden(opnd1);
  162. }
  163. else {
  164. /* check for zero */
  165. if (Sgl_iszero_mantissa(opnd1)) {
  166. Sgl_setzero_exponentmantissa(result);
  167. *dstptr = result;
  168. return(NOEXCEPTION);
  169. }
  170. /* is denormalized; want to normalize */
  171. Sgl_clear_signexponent(opnd1);
  172. Sgl_leftshiftby1(opnd1);
  173. Sgl_normalize(opnd1,dest_exponent);
  174. }
  175. /* opnd2 needs to have hidden bit set with msb in hidden bit */
  176. if (Sgl_isnotzero_exponent(opnd2)) {
  177. Sgl_clear_signexponent_set_hidden(opnd2);
  178. }
  179. else {
  180. /* is denormalized; want to normalize */
  181. Sgl_clear_signexponent(opnd2);
  182. Sgl_leftshiftby1(opnd2);
  183. while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
  184. Sgl_leftshiftby8(opnd2);
  185. dest_exponent += 8;
  186. }
  187. if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
  188. Sgl_leftshiftby4(opnd2);
  189. dest_exponent += 4;
  190. }
  191. while(Sgl_iszero_hidden(opnd2)) {
  192. Sgl_leftshiftby1(opnd2);
  193. dest_exponent += 1;
  194. }
  195. }
  196. /* Divide the source mantissas */
  197. /*
  198. * A non_restoring divide algorithm is used.
  199. */
  200. Sgl_subtract(opnd1,opnd2,opnd1);
  201. Sgl_setzero(opnd3);
  202. for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
  203. Sgl_leftshiftby1(opnd1);
  204. Sgl_leftshiftby1(opnd3);
  205. if (Sgl_iszero_sign(opnd1)) {
  206. Sgl_setone_lowmantissa(opnd3);
  207. Sgl_subtract(opnd1,opnd2,opnd1);
  208. }
  209. else Sgl_addition(opnd1,opnd2,opnd1);
  210. }
  211. if (count <= SGL_P) {
  212. Sgl_leftshiftby1(opnd3);
  213. Sgl_setone_lowmantissa(opnd3);
  214. Sgl_leftshift(opnd3,SGL_P-count);
  215. if (Sgl_iszero_hidden(opnd3)) {
  216. Sgl_leftshiftby1(opnd3);
  217. dest_exponent--;
  218. }
  219. }
  220. else {
  221. if (Sgl_iszero_hidden(opnd3)) {
  222. /* need to get one more bit of result */
  223. Sgl_leftshiftby1(opnd1);
  224. Sgl_leftshiftby1(opnd3);
  225. if (Sgl_iszero_sign(opnd1)) {
  226. Sgl_setone_lowmantissa(opnd3);
  227. Sgl_subtract(opnd1,opnd2,opnd1);
  228. }
  229. else Sgl_addition(opnd1,opnd2,opnd1);
  230. dest_exponent--;
  231. }
  232. if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
  233. stickybit = Sgl_all(opnd1);
  234. }
  235. inexact = guardbit | stickybit;
  236. /*
  237. * round result
  238. */
  239. if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
  240. Sgl_clear_signexponent(opnd3);
  241. switch (Rounding_mode()) {
  242. case ROUNDPLUS:
  243. if (Sgl_iszero_sign(result))
  244. Sgl_increment_mantissa(opnd3);
  245. break;
  246. case ROUNDMINUS:
  247. if (Sgl_isone_sign(result))
  248. Sgl_increment_mantissa(opnd3);
  249. break;
  250. case ROUNDNEAREST:
  251. if (guardbit) {
  252. if (stickybit || Sgl_isone_lowmantissa(opnd3))
  253. Sgl_increment_mantissa(opnd3);
  254. }
  255. }
  256. if (Sgl_isone_hidden(opnd3)) dest_exponent++;
  257. }
  258. Sgl_set_mantissa(result,opnd3);
  259. /*
  260. * Test for overflow
  261. */
  262. if (dest_exponent >= SGL_INFINITY_EXPONENT) {
  263. /* trap if OVERFLOWTRAP enabled */
  264. if (Is_overflowtrap_enabled()) {
  265. /*
  266. * Adjust bias of result
  267. */
  268. Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
  269. *dstptr = result;
  270. if (inexact)
  271. if (Is_inexacttrap_enabled())
  272. return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  273. else Set_inexactflag();
  274. return(OVERFLOWEXCEPTION);
  275. }
  276. Set_overflowflag();
  277. /* set result to infinity or largest number */
  278. Sgl_setoverflow(result);
  279. inexact = TRUE;
  280. }
  281. /*
  282. * Test for underflow
  283. */
  284. else if (dest_exponent <= 0) {
  285. /* trap if UNDERFLOWTRAP enabled */
  286. if (Is_underflowtrap_enabled()) {
  287. /*
  288. * Adjust bias of result
  289. */
  290. Sgl_setwrapped_exponent(result,dest_exponent,unfl);
  291. *dstptr = result;
  292. if (inexact)
  293. if (Is_inexacttrap_enabled())
  294. return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
  295. else Set_inexactflag();
  296. return(UNDERFLOWEXCEPTION);
  297. }
  298. /* Determine if should set underflow flag */
  299. is_tiny = TRUE;
  300. if (dest_exponent == 0 && inexact) {
  301. switch (Rounding_mode()) {
  302. case ROUNDPLUS:
  303. if (Sgl_iszero_sign(result)) {
  304. Sgl_increment(opnd3);
  305. if (Sgl_isone_hiddenoverflow(opnd3))
  306. is_tiny = FALSE;
  307. Sgl_decrement(opnd3);
  308. }
  309. break;
  310. case ROUNDMINUS:
  311. if (Sgl_isone_sign(result)) {
  312. Sgl_increment(opnd3);
  313. if (Sgl_isone_hiddenoverflow(opnd3))
  314. is_tiny = FALSE;
  315. Sgl_decrement(opnd3);
  316. }
  317. break;
  318. case ROUNDNEAREST:
  319. if (guardbit && (stickybit ||
  320. Sgl_isone_lowmantissa(opnd3))) {
  321. Sgl_increment(opnd3);
  322. if (Sgl_isone_hiddenoverflow(opnd3))
  323. is_tiny = FALSE;
  324. Sgl_decrement(opnd3);
  325. }
  326. break;
  327. }
  328. }
  329. /*
  330. * denormalize result or set to signed zero
  331. */
  332. stickybit = inexact;
  333. Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
  334. /* return rounded number */
  335. if (inexact) {
  336. switch (Rounding_mode()) {
  337. case ROUNDPLUS:
  338. if (Sgl_iszero_sign(result)) {
  339. Sgl_increment(opnd3);
  340. }
  341. break;
  342. case ROUNDMINUS:
  343. if (Sgl_isone_sign(result)) {
  344. Sgl_increment(opnd3);
  345. }
  346. break;
  347. case ROUNDNEAREST:
  348. if (guardbit && (stickybit ||
  349. Sgl_isone_lowmantissa(opnd3))) {
  350. Sgl_increment(opnd3);
  351. }
  352. break;
  353. }
  354. if (is_tiny) Set_underflowflag();
  355. }
  356. Sgl_set_exponentmantissa(result,opnd3);
  357. }
  358. else Sgl_set_exponent(result,dest_exponent);
  359. *dstptr = result;
  360. /* check for inexact */
  361. if (inexact) {
  362. if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  363. else Set_inexactflag();
  364. }
  365. return(NOEXCEPTION);
  366. }