1 /* $NetBSD: sfadd.c,v 1.4 2007/02/22 05:46:30 thorpej 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 <sys/cdefs.h> 45 __KERNEL_RCSID(0, "$NetBSD: sfadd.c,v 1.4 2007/02/22 05:46:30 thorpej Exp $"); 46 47 #include "../spmath/float.h" 48 #include "../spmath/sgl_float.h" 49 50 /* 51 * Single_add: add two single precision values. 52 */ 53 int 54 sgl_fadd(leftptr, rightptr, dstptr, status) 55 sgl_floating_point *leftptr, *rightptr, *dstptr; 56 unsigned int *status; 57 { 58 register unsigned int left, right, result, extent; 59 register unsigned int signless_upper_left, signless_upper_right, save; 60 61 62 register int result_exponent, right_exponent, diff_exponent; 63 register int sign_save, jumpsize; 64 register int inexact = false; 65 register int underflowtrap; 66 67 /* Create local copies of the numbers */ 68 left = *leftptr; 69 right = *rightptr; 70 71 /* A zero "save" helps discover equal operands (for later), * 72 * and is used in swapping operands (if needed). */ 73 Sgl_xortointp1(left,right,/*to*/save); 74 75 /* 76 * check first operand for NaN's or infinity 77 */ 78 if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT) 79 { 80 if (Sgl_iszero_mantissa(left)) 81 { 82 if (Sgl_isnotnan(right)) 83 { 84 if (Sgl_isinfinity(right) && save!=0) 85 { 86 /* 87 * invalid since operands are opposite signed infinity's 88 */ 89 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 90 Set_invalidflag(); 91 Sgl_makequietnan(result); 92 *dstptr = result; 93 return(NOEXCEPTION); 94 } 95 /* 96 * return infinity 97 */ 98 *dstptr = left; 99 return(NOEXCEPTION); 100 } 101 } 102 else 103 { 104 /* 105 * is NaN; signaling or quiet? 106 */ 107 if (Sgl_isone_signaling(left)) 108 { 109 /* trap if INVALIDTRAP enabled */ 110 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 111 /* make NaN quiet */ 112 Set_invalidflag(); 113 Sgl_set_quiet(left); 114 } 115 /* 116 * is second operand a signaling NaN? 117 */ 118 else if (Sgl_is_signalingnan(right)) 119 { 120 /* trap if INVALIDTRAP enabled */ 121 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 122 /* make NaN quiet */ 123 Set_invalidflag(); 124 Sgl_set_quiet(right); 125 *dstptr = right; 126 return(NOEXCEPTION); 127 } 128 /* 129 * return quiet NaN 130 */ 131 *dstptr = left; 132 return(NOEXCEPTION); 133 } 134 } /* End left NaN or Infinity processing */ 135 /* 136 * check second operand for NaN's or infinity 137 */ 138 if (Sgl_isinfinity_exponent(right)) 139 { 140 if (Sgl_iszero_mantissa(right)) 141 { 142 /* return infinity */ 143 *dstptr = right; 144 return(NOEXCEPTION); 145 } 146 /* 147 * is NaN; signaling or quiet? 148 */ 149 if (Sgl_isone_signaling(right)) 150 { 151 /* trap if INVALIDTRAP enabled */ 152 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 153 /* make NaN quiet */ 154 Set_invalidflag(); 155 Sgl_set_quiet(right); 156 } 157 /* 158 * return quiet NaN 159 */ 160 *dstptr = right; 161 return(NOEXCEPTION); 162 } /* End right NaN or Infinity processing */ 163 164 /* Invariant: Must be dealing with finite numbers */ 165 166 /* Compare operands by removing the sign */ 167 Sgl_copytoint_exponentmantissa(left,signless_upper_left); 168 Sgl_copytoint_exponentmantissa(right,signless_upper_right); 169 170 /* sign difference selects add or sub operation. */ 171 if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right)) 172 { 173 /* Set the left operand to the larger one by XOR swap * 174 * First finish the first word using "save" */ 175 Sgl_xorfromintp1(save,right,/*to*/right); 176 Sgl_xorfromintp1(save,left,/*to*/left); 177 result_exponent = Sgl_exponent(left); 178 } 179 /* Invariant: left is not smaller than right. */ 180 181 if((right_exponent = Sgl_exponent(right)) == 0) 182 { 183 /* Denormalized operands. First look for zeroes */ 184 if(Sgl_iszero_mantissa(right)) 185 { 186 /* right is zero */ 187 if(Sgl_iszero_exponentmantissa(left)) 188 { 189 /* Both operands are zeros */ 190 if(Is_rounding_mode(ROUNDMINUS)) 191 { 192 Sgl_or_signs(left,/*with*/right); 193 } 194 else 195 { 196 Sgl_and_signs(left,/*with*/right); 197 } 198 } 199 else 200 { 201 /* Left is not a zero and must be the result. Trapped 202 * underflows are signaled if left is denormalized. Result 203 * is always exact. */ 204 if( (result_exponent == 0) && Is_underflowtrap_enabled() ) 205 { 206 /* need to normalize results mantissa */ 207 sign_save = Sgl_signextendedsign(left); 208 Sgl_leftshiftby1(left); 209 Sgl_normalize(left,result_exponent); 210 Sgl_set_sign(left,/*using*/sign_save); 211 Sgl_setwrapped_exponent(left,result_exponent,unfl); 212 *dstptr = left; 213 return(UNDERFLOWEXCEPTION); 214 } 215 } 216 *dstptr = left; 217 return(NOEXCEPTION); 218 } 219 220 /* Neither are zeroes */ 221 Sgl_clear_sign(right); /* Exponent is already cleared */ 222 if(result_exponent == 0 ) 223 { 224 /* Both operands are denormalized. The result must be exact 225 * and is simply calculated. A sum could become normalized and a 226 * difference could cancel to a true zero. */ 227 if( (/*signed*/int) save < 0 ) 228 { 229 Sgl_subtract(left,/*minus*/right,/*into*/result); 230 if(Sgl_iszero_mantissa(result)) 231 { 232 if(Is_rounding_mode(ROUNDMINUS)) 233 { 234 Sgl_setone_sign(result); 235 } 236 else 237 { 238 Sgl_setzero_sign(result); 239 } 240 *dstptr = result; 241 return(NOEXCEPTION); 242 } 243 } 244 else 245 { 246 Sgl_addition(left,right,/*into*/result); 247 if(Sgl_isone_hidden(result)) 248 { 249 *dstptr = result; 250 return(NOEXCEPTION); 251 } 252 } 253 if(Is_underflowtrap_enabled()) 254 { 255 /* need to normalize result */ 256 sign_save = Sgl_signextendedsign(result); 257 Sgl_leftshiftby1(result); 258 Sgl_normalize(result,result_exponent); 259 Sgl_set_sign(result,/*using*/sign_save); 260 Sgl_setwrapped_exponent(result,result_exponent,unfl); 261 *dstptr = result; 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; 415 return(NOEXCEPTION); /* Sign bit is already set */ 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 ...add 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