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