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