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