sfsub.c 14 KB

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