sfadd.c 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505
  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/sfadd.c $Revision: 1.1 $
  13. *
  14. * Purpose:
  15. * Single_add: add two single precision values.
  16. *
  17. * External Interfaces:
  18. * sgl_fadd(leftptr, rightptr, 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_add: add two single precision values.
  31. */
  32. int
  33. sgl_fadd(
  34. sgl_floating_point *leftptr,
  35. sgl_floating_point *rightptr,
  36. sgl_floating_point *dstptr,
  37. unsigned int *status)
  38. {
  39. register unsigned int left, right, result, extent;
  40. register unsigned int signless_upper_left, signless_upper_right, save;
  41. register int result_exponent, right_exponent, diff_exponent;
  42. register int sign_save, jumpsize;
  43. register boolean inexact = FALSE;
  44. register boolean underflowtrap;
  45. /* Create local copies of the numbers */
  46. left = *leftptr;
  47. right = *rightptr;
  48. /* A zero "save" helps discover equal operands (for later), *
  49. * and is used in swapping operands (if needed). */
  50. Sgl_xortointp1(left,right,/*to*/save);
  51. /*
  52. * check first operand for NaN's or infinity
  53. */
  54. if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT)
  55. {
  56. if (Sgl_iszero_mantissa(left))
  57. {
  58. if (Sgl_isnotnan(right))
  59. {
  60. if (Sgl_isinfinity(right) && save!=0)
  61. {
  62. /*
  63. * invalid since operands are opposite signed infinity's
  64. */
  65. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  66. Set_invalidflag();
  67. Sgl_makequietnan(result);
  68. *dstptr = result;
  69. return(NOEXCEPTION);
  70. }
  71. /*
  72. * return infinity
  73. */
  74. *dstptr = left;
  75. return(NOEXCEPTION);
  76. }
  77. }
  78. else
  79. {
  80. /*
  81. * is NaN; signaling or quiet?
  82. */
  83. if (Sgl_isone_signaling(left))
  84. {
  85. /* trap if INVALIDTRAP enabled */
  86. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  87. /* make NaN quiet */
  88. Set_invalidflag();
  89. Sgl_set_quiet(left);
  90. }
  91. /*
  92. * is second operand a signaling NaN?
  93. */
  94. else if (Sgl_is_signalingnan(right))
  95. {
  96. /* trap if INVALIDTRAP enabled */
  97. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  98. /* make NaN quiet */
  99. Set_invalidflag();
  100. Sgl_set_quiet(right);
  101. *dstptr = right;
  102. return(NOEXCEPTION);
  103. }
  104. /*
  105. * return quiet NaN
  106. */
  107. *dstptr = left;
  108. return(NOEXCEPTION);
  109. }
  110. } /* End left NaN or Infinity processing */
  111. /*
  112. * check second operand for NaN's or infinity
  113. */
  114. if (Sgl_isinfinity_exponent(right))
  115. {
  116. if (Sgl_iszero_mantissa(right))
  117. {
  118. /* return infinity */
  119. *dstptr = right;
  120. return(NOEXCEPTION);
  121. }
  122. /*
  123. * is NaN; signaling or quiet?
  124. */
  125. if (Sgl_isone_signaling(right))
  126. {
  127. /* trap if INVALIDTRAP enabled */
  128. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  129. /* make NaN quiet */
  130. Set_invalidflag();
  131. Sgl_set_quiet(right);
  132. }
  133. /*
  134. * return quiet NaN
  135. */
  136. *dstptr = right;
  137. return(NOEXCEPTION);
  138. } /* End right NaN or Infinity processing */
  139. /* Invariant: Must be dealing with finite numbers */
  140. /* Compare operands by removing the sign */
  141. Sgl_copytoint_exponentmantissa(left,signless_upper_left);
  142. Sgl_copytoint_exponentmantissa(right,signless_upper_right);
  143. /* sign difference selects add or sub operation. */
  144. if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right))
  145. {
  146. /* Set the left operand to the larger one by XOR swap *
  147. * First finish the first word using "save" */
  148. Sgl_xorfromintp1(save,right,/*to*/right);
  149. Sgl_xorfromintp1(save,left,/*to*/left);
  150. result_exponent = Sgl_exponent(left);
  151. }
  152. /* Invariant: left is not smaller than right. */
  153. if((right_exponent = Sgl_exponent(right)) == 0)
  154. {
  155. /* Denormalized operands. First look for zeroes */
  156. if(Sgl_iszero_mantissa(right))
  157. {
  158. /* right is zero */
  159. if(Sgl_iszero_exponentmantissa(left))
  160. {
  161. /* Both operands are zeros */
  162. if(Is_rounding_mode(ROUNDMINUS))
  163. {
  164. Sgl_or_signs(left,/*with*/right);
  165. }
  166. else
  167. {
  168. Sgl_and_signs(left,/*with*/right);
  169. }
  170. }
  171. else
  172. {
  173. /* Left is not a zero and must be the result. Trapped
  174. * underflows are signaled if left is denormalized. Result
  175. * is always exact. */
  176. if( (result_exponent == 0) && Is_underflowtrap_enabled() )
  177. {
  178. /* need to normalize results mantissa */
  179. sign_save = Sgl_signextendedsign(left);
  180. Sgl_leftshiftby1(left);
  181. Sgl_normalize(left,result_exponent);
  182. Sgl_set_sign(left,/*using*/sign_save);
  183. Sgl_setwrapped_exponent(left,result_exponent,unfl);
  184. *dstptr = left;
  185. return(UNDERFLOWEXCEPTION);
  186. }
  187. }
  188. *dstptr = left;
  189. return(NOEXCEPTION);
  190. }
  191. /* Neither are zeroes */
  192. Sgl_clear_sign(right); /* Exponent is already cleared */
  193. if(result_exponent == 0 )
  194. {
  195. /* Both operands are denormalized. The result must be exact
  196. * and is simply calculated. A sum could become normalized and a
  197. * difference could cancel to a true zero. */
  198. if( (/*signed*/int) save < 0 )
  199. {
  200. Sgl_subtract(left,/*minus*/right,/*into*/result);
  201. if(Sgl_iszero_mantissa(result))
  202. {
  203. if(Is_rounding_mode(ROUNDMINUS))
  204. {
  205. Sgl_setone_sign(result);
  206. }
  207. else
  208. {
  209. Sgl_setzero_sign(result);
  210. }
  211. *dstptr = result;
  212. return(NOEXCEPTION);
  213. }
  214. }
  215. else
  216. {
  217. Sgl_addition(left,right,/*into*/result);
  218. if(Sgl_isone_hidden(result))
  219. {
  220. *dstptr = result;
  221. return(NOEXCEPTION);
  222. }
  223. }
  224. if(Is_underflowtrap_enabled())
  225. {
  226. /* need to normalize result */
  227. sign_save = Sgl_signextendedsign(result);
  228. Sgl_leftshiftby1(result);
  229. Sgl_normalize(result,result_exponent);
  230. Sgl_set_sign(result,/*using*/sign_save);
  231. Sgl_setwrapped_exponent(result,result_exponent,unfl);
  232. *dstptr = result;
  233. return(UNDERFLOWEXCEPTION);
  234. }
  235. *dstptr = result;
  236. return(NOEXCEPTION);
  237. }
  238. right_exponent = 1; /* Set exponent to reflect different bias
  239. * with denormalized numbers. */
  240. }
  241. else
  242. {
  243. Sgl_clear_signexponent_set_hidden(right);
  244. }
  245. Sgl_clear_exponent_set_hidden(left);
  246. diff_exponent = result_exponent - right_exponent;
  247. /*
  248. * Special case alignment of operands that would force alignment
  249. * beyond the extent of the extension. A further optimization
  250. * could special case this but only reduces the path length for this
  251. * infrequent case.
  252. */
  253. if(diff_exponent > SGL_THRESHOLD)
  254. {
  255. diff_exponent = SGL_THRESHOLD;
  256. }
  257. /* Align right operand by shifting to right */
  258. Sgl_right_align(/*operand*/right,/*shifted by*/diff_exponent,
  259. /*and lower to*/extent);
  260. /* Treat sum and difference of the operands separately. */
  261. if( (/*signed*/int) save < 0 )
  262. {
  263. /*
  264. * Difference of the two operands. Their can be no overflow. A
  265. * borrow can occur out of the hidden bit and force a post
  266. * normalization phase.
  267. */
  268. Sgl_subtract_withextension(left,/*minus*/right,/*with*/extent,/*into*/result);
  269. if(Sgl_iszero_hidden(result))
  270. {
  271. /* Handle normalization */
  272. /* A straightforward algorithm would now shift the result
  273. * and extension left until the hidden bit becomes one. Not
  274. * all of the extension bits need participate in the shift.
  275. * Only the two most significant bits (round and guard) are
  276. * needed. If only a single shift is needed then the guard
  277. * bit becomes a significant low order bit and the extension
  278. * must participate in the rounding. If more than a single
  279. * shift is needed, then all bits to the right of the guard
  280. * bit are zeros, and the guard bit may or may not be zero. */
  281. sign_save = Sgl_signextendedsign(result);
  282. Sgl_leftshiftby1_withextent(result,extent,result);
  283. /* Need to check for a zero result. The sign and exponent
  284. * fields have already been zeroed. The more efficient test
  285. * of the full object can be used.
  286. */
  287. if(Sgl_iszero(result))
  288. /* Must have been "x-x" or "x+(-x)". */
  289. {
  290. if(Is_rounding_mode(ROUNDMINUS)) Sgl_setone_sign(result);
  291. *dstptr = result;
  292. return(NOEXCEPTION);
  293. }
  294. result_exponent--;
  295. /* Look to see if normalization is finished. */
  296. if(Sgl_isone_hidden(result))
  297. {
  298. if(result_exponent==0)
  299. {
  300. /* Denormalized, exponent should be zero. Left operand *
  301. * was normalized, so extent (guard, round) was zero */
  302. goto underflow;
  303. }
  304. else
  305. {
  306. /* No further normalization is needed. */
  307. Sgl_set_sign(result,/*using*/sign_save);
  308. Ext_leftshiftby1(extent);
  309. goto round;
  310. }
  311. }
  312. /* Check for denormalized, exponent should be zero. Left *
  313. * operand was normalized, so extent (guard, round) was zero */
  314. if(!(underflowtrap = Is_underflowtrap_enabled()) &&
  315. result_exponent==0) goto underflow;
  316. /* Shift extension to complete one bit of normalization and
  317. * update exponent. */
  318. Ext_leftshiftby1(extent);
  319. /* Discover first one bit to determine shift amount. Use a
  320. * modified binary search. We have already shifted the result
  321. * one position right and still not found a one so the remainder
  322. * of the extension must be zero and simplifies rounding. */
  323. /* Scan bytes */
  324. while(Sgl_iszero_hiddenhigh7mantissa(result))
  325. {
  326. Sgl_leftshiftby8(result);
  327. if((result_exponent -= 8) <= 0 && !underflowtrap)
  328. goto underflow;
  329. }
  330. /* Now narrow it down to the nibble */
  331. if(Sgl_iszero_hiddenhigh3mantissa(result))
  332. {
  333. /* The lower nibble contains the normalizing one */
  334. Sgl_leftshiftby4(result);
  335. if((result_exponent -= 4) <= 0 && !underflowtrap)
  336. goto underflow;
  337. }
  338. /* Select case were first bit is set (already normalized)
  339. * otherwise select the proper shift. */
  340. if((jumpsize = Sgl_hiddenhigh3mantissa(result)) > 7)
  341. {
  342. /* Already normalized */
  343. if(result_exponent <= 0) goto underflow;
  344. Sgl_set_sign(result,/*using*/sign_save);
  345. Sgl_set_exponent(result,/*using*/result_exponent);
  346. *dstptr = result;
  347. return(NOEXCEPTION);
  348. }
  349. Sgl_sethigh4bits(result,/*using*/sign_save);
  350. switch(jumpsize)
  351. {
  352. case 1:
  353. {
  354. Sgl_leftshiftby3(result);
  355. result_exponent -= 3;
  356. break;
  357. }
  358. case 2:
  359. case 3:
  360. {
  361. Sgl_leftshiftby2(result);
  362. result_exponent -= 2;
  363. break;
  364. }
  365. case 4:
  366. case 5:
  367. case 6:
  368. case 7:
  369. {
  370. Sgl_leftshiftby1(result);
  371. result_exponent -= 1;
  372. break;
  373. }
  374. }
  375. if(result_exponent > 0)
  376. {
  377. Sgl_set_exponent(result,/*using*/result_exponent);
  378. *dstptr = result;
  379. return(NOEXCEPTION); /* Sign bit is already set */
  380. }
  381. /* Fixup potential underflows */
  382. underflow:
  383. if(Is_underflowtrap_enabled())
  384. {
  385. Sgl_set_sign(result,sign_save);
  386. Sgl_setwrapped_exponent(result,result_exponent,unfl);
  387. *dstptr = result;
  388. /* inexact = FALSE; */
  389. return(UNDERFLOWEXCEPTION);
  390. }
  391. /*
  392. * Since we cannot get an inexact denormalized result,
  393. * we can now return.
  394. */
  395. Sgl_right_align(result,/*by*/(1-result_exponent),extent);
  396. Sgl_clear_signexponent(result);
  397. Sgl_set_sign(result,sign_save);
  398. *dstptr = result;
  399. return(NOEXCEPTION);
  400. } /* end if(hidden...)... */
  401. /* Fall through and round */
  402. } /* end if(save < 0)... */
  403. else
  404. {
  405. /* Add magnitudes */
  406. Sgl_addition(left,right,/*to*/result);
  407. if(Sgl_isone_hiddenoverflow(result))
  408. {
  409. /* Prenormalization required. */
  410. Sgl_rightshiftby1_withextent(result,extent,extent);
  411. Sgl_arithrightshiftby1(result);
  412. result_exponent++;
  413. } /* end if hiddenoverflow... */
  414. } /* end else ...add magnitudes... */
  415. /* Round the result. If the extension is all zeros,then the result is
  416. * exact. Otherwise round in the correct direction. No underflow is
  417. * possible. If a postnormalization is necessary, then the mantissa is
  418. * all zeros so no shift is needed. */
  419. round:
  420. if(Ext_isnotzero(extent))
  421. {
  422. inexact = TRUE;
  423. switch(Rounding_mode())
  424. {
  425. case ROUNDNEAREST: /* The default. */
  426. if(Ext_isone_sign(extent))
  427. {
  428. /* at least 1/2 ulp */
  429. if(Ext_isnotzero_lower(extent) ||
  430. Sgl_isone_lowmantissa(result))
  431. {
  432. /* either exactly half way and odd or more than 1/2ulp */
  433. Sgl_increment(result);
  434. }
  435. }
  436. break;
  437. case ROUNDPLUS:
  438. if(Sgl_iszero_sign(result))
  439. {
  440. /* Round up positive results */
  441. Sgl_increment(result);
  442. }
  443. break;
  444. case ROUNDMINUS:
  445. if(Sgl_isone_sign(result))
  446. {
  447. /* Round down negative results */
  448. Sgl_increment(result);
  449. }
  450. case ROUNDZERO:;
  451. /* truncate is simple */
  452. } /* end switch... */
  453. if(Sgl_isone_hiddenoverflow(result)) result_exponent++;
  454. }
  455. if(result_exponent == SGL_INFINITY_EXPONENT)
  456. {
  457. /* Overflow */
  458. if(Is_overflowtrap_enabled())
  459. {
  460. Sgl_setwrapped_exponent(result,result_exponent,ovfl);
  461. *dstptr = result;
  462. if (inexact)
  463. if (Is_inexacttrap_enabled())
  464. return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  465. else Set_inexactflag();
  466. return(OVERFLOWEXCEPTION);
  467. }
  468. else
  469. {
  470. Set_overflowflag();
  471. inexact = TRUE;
  472. Sgl_setoverflow(result);
  473. }
  474. }
  475. else Sgl_set_exponent(result,result_exponent);
  476. *dstptr = result;
  477. if(inexact)
  478. if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  479. else Set_inexactflag();
  480. return(NOEXCEPTION);
  481. }