1 /* $NetBSD: dfadd.c,v 1.1 2002/06/05 01:04:24 fredette Exp $ */ 2 3 /* $OpenBSD: dfadd.c,v 1.4 2001/03/29 03:58:17 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 * pmk1.1 27 */ 28 /* 29 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 30 * 31 * To anyone who acknowledges that this file is provided "AS IS" 32 * without any express or implied warranty: 33 * permission to use, copy, modify, and distribute this file 34 * for any purpose is hereby granted without fee, provided that 35 * the above copyright notice and this notice appears in all 36 * copies, and that the name of Hewlett-Packard Company not be 37 * used in advertising or publicity pertaining to distribution 38 * of the software without specific, written prior permission. 39 * Hewlett-Packard Company makes no representations about the 40 * suitability of this software for any purpose. 41 */ 42 43 #include "../spmath/float.h" 44 #include "../spmath/dbl_float.h" 45 46 /* 47 * Double_add: add two double precision values. 48 */ 49 int 50 dbl_fadd(leftptr, rightptr, dstptr, status) 51 dbl_floating_point *leftptr, *rightptr, *dstptr; 52 unsigned int *status; 53 { 54 register unsigned int signless_upper_left, signless_upper_right, save; 55 register unsigned int leftp1, leftp2, rightp1, rightp2, extent; 56 register unsigned int resultp1 = 0, resultp2 = 0; 57 58 register int result_exponent, right_exponent, diff_exponent; 59 register int sign_save, jumpsize; 60 register int inexact = FALSE; 61 register int underflowtrap; 62 63 /* Create local copies of the numbers */ 64 Dbl_copyfromptr(leftptr,leftp1,leftp2); 65 Dbl_copyfromptr(rightptr,rightp1,rightp2); 66 67 /* A zero "save" helps discover equal operands (for later), * 68 * and is used in swapping operands (if needed). */ 69 Dbl_xortointp1(leftp1,rightp1,/*to*/save); 70 71 /* 72 * check first operand for NaN's or infinity 73 */ 74 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT) 75 { 76 if (Dbl_iszero_mantissa(leftp1,leftp2)) 77 { 78 if (Dbl_isnotnan(rightp1,rightp2)) 79 { 80 if (Dbl_isinfinity(rightp1,rightp2) && save!=0) 81 { 82 /* 83 * invalid since operands are opposite signed infinity's 84 */ 85 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 86 Set_invalidflag(); 87 Dbl_makequietnan(resultp1,resultp2); 88 Dbl_copytoptr(resultp1,resultp2,dstptr); 89 return(NOEXCEPTION); 90 } 91 /* 92 * return infinity 93 */ 94 Dbl_copytoptr(leftp1,leftp2,dstptr); 95 return(NOEXCEPTION); 96 } 97 } 98 else 99 { 100 /* 101 * is NaN; signaling or quiet? 102 */ 103 if (Dbl_isone_signaling(leftp1)) 104 { 105 /* trap if INVALIDTRAP enabled */ 106 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 107 /* make NaN quiet */ 108 Set_invalidflag(); 109 Dbl_set_quiet(leftp1); 110 } 111 /* 112 * is second operand a signaling NaN? 113 */ 114 else if (Dbl_is_signalingnan(rightp1)) 115 { 116 /* trap if INVALIDTRAP enabled */ 117 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 118 /* make NaN quiet */ 119 Set_invalidflag(); 120 Dbl_set_quiet(rightp1); 121 Dbl_copytoptr(rightp1,rightp2,dstptr); 122 return(NOEXCEPTION); 123 } 124 /* 125 * return quiet NaN 126 */ 127 Dbl_copytoptr(leftp1,leftp2,dstptr); 128 return(NOEXCEPTION); 129 } 130 } /* End left NaN or Infinity processing */ 131 /* 132 * check second operand for NaN's or infinity 133 */ 134 if (Dbl_isinfinity_exponent(rightp1)) 135 { 136 if (Dbl_iszero_mantissa(rightp1,rightp2)) 137 { 138 /* return infinity */ 139 Dbl_copytoptr(rightp1,rightp2,dstptr); 140 return(NOEXCEPTION); 141 } 142 /* 143 * is NaN; signaling or quiet? 144 */ 145 if (Dbl_isone_signaling(rightp1)) 146 { 147 /* trap if INVALIDTRAP enabled */ 148 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 149 /* make NaN quiet */ 150 Set_invalidflag(); 151 Dbl_set_quiet(rightp1); 152 } 153 /* 154 * return quiet NaN 155 */ 156 Dbl_copytoptr(rightp1,rightp2,dstptr); 157 return(NOEXCEPTION); 158 } /* End right NaN or Infinity processing */ 159 160 /* Invariant: Must be dealing with finite numbers */ 161 162 /* Compare operands by removing the sign */ 163 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left); 164 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right); 165 166 /* sign difference selects add or sub operation. */ 167 if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right)) 168 { 169 /* Set the left operand to the larger one by XOR swap * 170 * First finish the first word using "save" */ 171 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1); 172 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1); 173 Dbl_swap_lower(leftp2,rightp2); 174 result_exponent = Dbl_exponent(leftp1); 175 } 176 /* Invariant: left is not smaller than right. */ 177 178 if((right_exponent = Dbl_exponent(rightp1)) == 0) 179 { 180 /* Denormalized operands. First look for zeroes */ 181 if(Dbl_iszero_mantissa(rightp1,rightp2)) 182 { 183 /* right is zero */ 184 if(Dbl_iszero_exponentmantissa(leftp1,leftp2)) 185 { 186 /* Both operands are zeros */ 187 if(Is_rounding_mode(ROUNDMINUS)) 188 { 189 Dbl_or_signs(leftp1,/*with*/rightp1); 190 } 191 else 192 { 193 Dbl_and_signs(leftp1,/*with*/rightp1); 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 = Dbl_signextendedsign(leftp1); 205 Dbl_leftshiftby1(leftp1,leftp2); 206 Dbl_normalize(leftp1,leftp2,result_exponent); 207 Dbl_set_sign(leftp1,/*using*/sign_save); 208 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl); 209 Dbl_copytoptr(leftp1,leftp2,dstptr); 210 /* inexact = FALSE */ 211 return(UNDERFLOWEXCEPTION); 212 } 213 } 214 Dbl_copytoptr(leftp1,leftp2,dstptr); 215 return(NOEXCEPTION); 216 } 217 218 /* Neither are zeroes */ 219 Dbl_clear_sign(rightp1); /* Exponent is already cleared */ 220 if(result_exponent == 0 ) 221 { 222 /* Both operands are denormalized. The result must be exact 223 * and is simply calculated. A sum could become normalized and a 224 * difference could cancel to a true zero. */ 225 if( (/*signed*/int) save < 0 ) 226 { 227 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2, 228 /*into*/resultp1,resultp2); 229 if(Dbl_iszero_mantissa(resultp1,resultp2)) 230 { 231 if(Is_rounding_mode(ROUNDMINUS)) 232 { 233 Dbl_setone_sign(resultp1); 234 } 235 else 236 { 237 Dbl_setzero_sign(resultp1); 238 } 239 Dbl_copytoptr(resultp1,resultp2,dstptr); 240 return(NOEXCEPTION); 241 } 242 } 243 else 244 { 245 Dbl_addition(leftp1,leftp2,rightp1,rightp2, 246 /*into*/resultp1,resultp2); 247 if(Dbl_isone_hidden(resultp1)) 248 { 249 Dbl_copytoptr(resultp1,resultp2,dstptr); 250 return(NOEXCEPTION); 251 } 252 } 253 if(Is_underflowtrap_enabled()) 254 { 255 /* need to normalize result */ 256 sign_save = Dbl_signextendedsign(resultp1); 257 Dbl_leftshiftby1(resultp1,resultp2); 258 Dbl_normalize(resultp1,resultp2,result_exponent); 259 Dbl_set_sign(resultp1,/*using*/sign_save); 260 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl); 261 Dbl_copytoptr(resultp1,resultp2,dstptr); 262 /* inexact = FALSE */ 263 return(UNDERFLOWEXCEPTION); 264 } 265 Dbl_copytoptr(resultp1,resultp2,dstptr); 266 return(NOEXCEPTION); 267 } 268 right_exponent = 1; /* Set exponent to reflect different bias 269 * with denomalized numbers. */ 270 } 271 else 272 { 273 Dbl_clear_signexponent_set_hidden(rightp1); 274 } 275 Dbl_clear_exponent_set_hidden(leftp1); 276 diff_exponent = result_exponent - right_exponent; 277 278 /* 279 * Special case alignment of operands that would force alignment 280 * beyond the extent of the extension. A further optimization 281 * could special case this but only reduces the path length for this 282 * infrequent case. 283 */ 284 if(diff_exponent > DBL_THRESHOLD) 285 { 286 diff_exponent = DBL_THRESHOLD; 287 } 288 289 /* Align right operand by shifting to right */ 290 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent, 291 /*and lower to*/extent); 292 293 /* Treat sum and difference of the operands separately. */ 294 if( (/*signed*/int) save < 0 ) 295 { 296 /* 297 * Difference of the two operands. Their can be no overflow. A 298 * borrow can occur out of the hidden bit and force a post 299 * normalization phase. 300 */ 301 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2, 302 /*with*/extent,/*into*/resultp1,resultp2); 303 if(Dbl_iszero_hidden(resultp1)) 304 { 305 /* Handle normalization */ 306 /* A straight foward algorithm would now shift the result 307 * and extension left until the hidden bit becomes one. Not 308 * all of the extension bits need participate in the shift. 309 * Only the two most significant bits (round and guard) are 310 * needed. If only a single shift is needed then the guard 311 * bit becomes a significant low order bit and the extension 312 * must participate in the rounding. If more than a single 313 * shift is needed, then all bits to the right of the guard 314 * bit are zeros, and the guard bit may or may not be zero. */ 315 sign_save = Dbl_signextendedsign(resultp1); 316 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2); 317 318 /* Need to check for a zero result. The sign and exponent 319 * fields have already been zeroed. The more efficient test 320 * of the full object can be used. 321 */ 322 if(Dbl_iszero(resultp1,resultp2)) 323 /* Must have been "x-x" or "x+(-x)". */ 324 { 325 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1); 326 Dbl_copytoptr(resultp1,resultp2,dstptr); 327 return(NOEXCEPTION); 328 } 329 result_exponent--; 330 /* Look to see if normalization is finished. */ 331 if(Dbl_isone_hidden(resultp1)) 332 { 333 if(result_exponent==0) 334 { 335 /* Denormalized, exponent should be zero. Left operand * 336 * was normalized, so extent (guard, round) was zero */ 337 goto underflow; 338 } 339 else 340 { 341 /* No further normalization is needed. */ 342 Dbl_set_sign(resultp1,/*using*/sign_save); 343 Ext_leftshiftby1(extent); 344 goto round; 345 } 346 } 347 348 /* Check for denormalized, exponent should be zero. Left * 349 * operand was normalized, so extent (guard, round) was zero */ 350 if(!(underflowtrap = Is_underflowtrap_enabled()) && 351 result_exponent==0) goto underflow; 352 353 /* Shift extension to complete one bit of normalization and 354 * update exponent. */ 355 Ext_leftshiftby1(extent); 356 357 /* Discover first one bit to determine shift amount. Use a 358 * modified binary search. We have already shifted the result 359 * one position right and still not found a one so the remainder 360 * of the extension must be zero and simplifies rounding. */ 361 /* Scan bytes */ 362 while(Dbl_iszero_hiddenhigh7mantissa(resultp1)) 363 { 364 Dbl_leftshiftby8(resultp1,resultp2); 365 if((result_exponent -= 8) <= 0 && !underflowtrap) 366 goto underflow; 367 } 368 /* Now narrow it down to the nibble */ 369 if(Dbl_iszero_hiddenhigh3mantissa(resultp1)) 370 { 371 /* The lower nibble contains the normalizing one */ 372 Dbl_leftshiftby4(resultp1,resultp2); 373 if((result_exponent -= 4) <= 0 && !underflowtrap) 374 goto underflow; 375 } 376 /* Select case were first bit is set (already normalized) 377 * otherwise select the proper shift. */ 378 if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7) 379 { 380 /* Already normalized */ 381 if(result_exponent <= 0) goto underflow; 382 Dbl_set_sign(resultp1,/*using*/sign_save); 383 Dbl_set_exponent(resultp1,/*using*/result_exponent); 384 Dbl_copytoptr(resultp1,resultp2,dstptr); 385 return(NOEXCEPTION); 386 } 387 Dbl_sethigh4bits(resultp1,/*using*/sign_save); 388 switch(jumpsize) 389 { 390 case 1: 391 { 392 Dbl_leftshiftby3(resultp1,resultp2); 393 result_exponent -= 3; 394 break; 395 } 396 case 2: 397 case 3: 398 { 399 Dbl_leftshiftby2(resultp1,resultp2); 400 result_exponent -= 2; 401 break; 402 } 403 case 4: 404 case 5: 405 case 6: 406 case 7: 407 { 408 Dbl_leftshiftby1(resultp1,resultp2); 409 result_exponent -= 1; 410 break; 411 } 412 } 413 if(result_exponent > 0) 414 { 415 Dbl_set_exponent(resultp1,/*using*/result_exponent); 416 Dbl_copytoptr(resultp1,resultp2,dstptr); 417 return(NOEXCEPTION); /* Sign bit is already set */ 418 } 419 /* Fixup potential underflows */ 420 underflow: 421 if(Is_underflowtrap_enabled()) 422 { 423 Dbl_set_sign(resultp1,sign_save); 424 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl); 425 Dbl_copytoptr(resultp1,resultp2,dstptr); 426 /* inexact = FALSE */ 427 return(UNDERFLOWEXCEPTION); 428 } 429 /* 430 * Since we cannot get an inexact denormalized result, 431 * we can now return. 432 */ 433 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent); 434 Dbl_clear_signexponent(resultp1); 435 Dbl_set_sign(resultp1,sign_save); 436 Dbl_copytoptr(resultp1,resultp2,dstptr); 437 return(NOEXCEPTION); 438 } /* end if(hidden...)... */ 439 /* Fall through and round */ 440 } /* end if(save < 0)... */ 441 else 442 { 443 /* Add magnitudes */ 444 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2); 445 if(Dbl_isone_hiddenoverflow(resultp1)) 446 { 447 /* Prenormalization required. */ 448 Dbl_rightshiftby1_withextent(resultp2,extent,extent); 449 Dbl_arithrightshiftby1(resultp1,resultp2); 450 result_exponent++; 451 } /* end if hiddenoverflow... */ 452 } /* end else ...add magnitudes... */ 453 454 /* Round the result. If the extension is all zeros,then the result is 455 * exact. Otherwise round in the correct direction. No underflow is 456 * possible. If a postnormalization is necessary, then the mantissa is 457 * all zeros so no shift is needed. */ 458 round: 459 if(Ext_isnotzero(extent)) 460 { 461 inexact = TRUE; 462 switch(Rounding_mode()) 463 { 464 case ROUNDNEAREST: /* The default. */ 465 if(Ext_isone_sign(extent)) 466 { 467 /* at least 1/2 ulp */ 468 if(Ext_isnotzero_lower(extent) || 469 Dbl_isone_lowmantissap2(resultp2)) 470 { 471 /* either exactly half way and odd or more than 1/2ulp */ 472 Dbl_increment(resultp1,resultp2); 473 } 474 } 475 break; 476 477 case ROUNDPLUS: 478 if(Dbl_iszero_sign(resultp1)) 479 { 480 /* Round up positive results */ 481 Dbl_increment(resultp1,resultp2); 482 } 483 break; 484 485 case ROUNDMINUS: 486 if(Dbl_isone_sign(resultp1)) 487 { 488 /* Round down negative results */ 489 Dbl_increment(resultp1,resultp2); 490 } 491 492 case ROUNDZERO:; 493 /* truncate is simple */ 494 } /* end switch... */ 495 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++; 496 } 497 if(result_exponent == DBL_INFINITY_EXPONENT) 498 { 499 /* Overflow */ 500 if(Is_overflowtrap_enabled()) 501 { 502 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl); 503 Dbl_copytoptr(resultp1,resultp2,dstptr); 504 if (inexact) { 505 if (Is_inexacttrap_enabled()) 506 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 507 else 508 Set_inexactflag(); 509 } 510 return(OVERFLOWEXCEPTION); 511 } 512 else 513 { 514 inexact = TRUE; 515 Set_overflowflag(); 516 Dbl_setoverflow(resultp1,resultp2); 517 } 518 } 519 else Dbl_set_exponent(resultp1,result_exponent); 520 Dbl_copytoptr(resultp1,resultp2,dstptr); 521 if(inexact) { 522 if(Is_inexacttrap_enabled()) 523 return(INEXACTEXCEPTION); 524 else 525 Set_inexactflag(); 526 } 527 return(NOEXCEPTION); 528 } 529