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