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