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