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