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