sfmpy.c 9.7 KB

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