1 /* $OpenBSD: dfadd.c,v 1.4 2001/03/29 03:58:17 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 * pmk1.1 25 */ 26 /* 27 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 28 * 29 * To anyone who acknowledges that this file is provided "AS IS" 30 * without any express or implied warranty: 31 * permission to use, copy, modify, and distribute this file 32 * for any purpose is hereby granted without fee, provided that 33 * the above copyright notice and this notice appears in all 34 * copies, and that the name of Hewlett-Packard Company not be 35 * used in advertising or publicity pertaining to distribution 36 * of the software without specific, written prior permission. 37 * Hewlett-Packard Company makes no representations about the 38 * suitability of this software for any purpose. 39 */ 40 41 #include "../spmath/float.h" 42 #include "../spmath/dbl_float.h" 43 44 /* 45 * Double_add: add two double precision values. 46 */ 47 int 48 dbl_fadd(leftptr, rightptr, dstptr, status) 49 dbl_floating_point *leftptr, *rightptr, *dstptr; 50 unsigned int *status; 51 { 52 register unsigned int signless_upper_left, signless_upper_right, save; 53 register unsigned int leftp1, leftp2, rightp1, rightp2, extent; 54 register unsigned int resultp1 = 0, resultp2 = 0; 55 56 register int result_exponent, right_exponent, diff_exponent; 57 register int sign_save, jumpsize; 58 register int inexact = FALSE; 59 register int underflowtrap; 60 61 /* Create local copies of the numbers */ 62 Dbl_copyfromptr(leftptr,leftp1,leftp2); 63 Dbl_copyfromptr(rightptr,rightp1,rightp2); 64 65 /* A zero "save" helps discover equal operands (for later), * 66 * and is used in swapping operands (if needed). */ 67 Dbl_xortointp1(leftp1,rightp1,/*to*/save); 68 69 /* 70 * check first operand for NaN's or infinity 71 */ 72 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT) 73 { 74 if (Dbl_iszero_mantissa(leftp1,leftp2)) 75 { 76 if (Dbl_isnotnan(rightp1,rightp2)) 77 { 78 if (Dbl_isinfinity(rightp1,rightp2) && save!=0) 79 { 80 /* 81 * invalid since operands are opposite signed infinity's 82 */ 83 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 84 Set_invalidflag(); 85 Dbl_makequietnan(resultp1,resultp2); 86 Dbl_copytoptr(resultp1,resultp2,dstptr); 87 return(NOEXCEPTION); 88 } 89 /* 90 * return infinity 91 */ 92 Dbl_copytoptr(leftp1,leftp2,dstptr); 93 return(NOEXCEPTION); 94 } 95 } 96 else 97 { 98 /* 99 * is NaN; signaling or quiet? 100 */ 101 if (Dbl_isone_signaling(leftp1)) 102 { 103 /* trap if INVALIDTRAP enabled */ 104 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 105 /* make NaN quiet */ 106 Set_invalidflag(); 107 Dbl_set_quiet(leftp1); 108 } 109 /* 110 * is second operand a signaling NaN? 111 */ 112 else if (Dbl_is_signalingnan(rightp1)) 113 { 114 /* trap if INVALIDTRAP enabled */ 115 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 116 /* make NaN quiet */ 117 Set_invalidflag(); 118 Dbl_set_quiet(rightp1); 119 Dbl_copytoptr(rightp1,rightp2,dstptr); 120 return(NOEXCEPTION); 121 } 122 /* 123 * return quiet NaN 124 */ 125 Dbl_copytoptr(leftp1,leftp2,dstptr); 126 return(NOEXCEPTION); 127 } 128 } /* End left NaN or Infinity processing */ 129 /* 130 * check second operand for NaN's or infinity 131 */ 132 if (Dbl_isinfinity_exponent(rightp1)) 133 { 134 if (Dbl_iszero_mantissa(rightp1,rightp2)) 135 { 136 /* return infinity */ 137 Dbl_copytoptr(rightp1,rightp2,dstptr); 138 return(NOEXCEPTION); 139 } 140 /* 141 * is NaN; signaling or quiet? 142 */ 143 if (Dbl_isone_signaling(rightp1)) 144 { 145 /* trap if INVALIDTRAP enabled */ 146 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 147 /* make NaN quiet */ 148 Set_invalidflag(); 149 Dbl_set_quiet(rightp1); 150 } 151 /* 152 * return quiet NaN 153 */ 154 Dbl_copytoptr(rightp1,rightp2,dstptr); 155 return(NOEXCEPTION); 156 } /* End right NaN or Infinity processing */ 157 158 /* Invariant: Must be dealing with finite numbers */ 159 160 /* Compare operands by removing the sign */ 161 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left); 162 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right); 163 164 /* sign difference selects add or sub operation. */ 165 if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right)) 166 { 167 /* Set the left operand to the larger one by XOR swap * 168 * First finish the first word using "save" */ 169 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1); 170 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1); 171 Dbl_swap_lower(leftp2,rightp2); 172 result_exponent = Dbl_exponent(leftp1); 173 } 174 /* Invariant: left is not smaller than right. */ 175 176 if((right_exponent = Dbl_exponent(rightp1)) == 0) 177 { 178 /* Denormalized operands. First look for zeroes */ 179 if(Dbl_iszero_mantissa(rightp1,rightp2)) 180 { 181 /* right is zero */ 182 if(Dbl_iszero_exponentmantissa(leftp1,leftp2)) 183 { 184 /* Both operands are zeros */ 185 if(Is_rounding_mode(ROUNDMINUS)) 186 { 187 Dbl_or_signs(leftp1,/*with*/rightp1); 188 } 189 else 190 { 191 Dbl_and_signs(leftp1,/*with*/rightp1); 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 = Dbl_signextendedsign(leftp1); 203 Dbl_leftshiftby1(leftp1,leftp2); 204 Dbl_normalize(leftp1,leftp2,result_exponent); 205 Dbl_set_sign(leftp1,/*using*/sign_save); 206 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl); 207 Dbl_copytoptr(leftp1,leftp2,dstptr); 208 /* inexact = FALSE */ 209 return(UNDERFLOWEXCEPTION); 210 } 211 } 212 Dbl_copytoptr(leftp1,leftp2,dstptr); 213 return(NOEXCEPTION); 214 } 215 216 /* Neither are zeroes */ 217 Dbl_clear_sign(rightp1); /* Exponent is already cleared */ 218 if(result_exponent == 0 ) 219 { 220 /* Both operands are denormalized. The result must be exact 221 * and is simply calculated. A sum could become normalized and a 222 * difference could cancel to a true zero. */ 223 if( (/*signed*/int) save < 0 ) 224 { 225 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2, 226 /*into*/resultp1,resultp2); 227 if(Dbl_iszero_mantissa(resultp1,resultp2)) 228 { 229 if(Is_rounding_mode(ROUNDMINUS)) 230 { 231 Dbl_setone_sign(resultp1); 232 } 233 else 234 { 235 Dbl_setzero_sign(resultp1); 236 } 237 Dbl_copytoptr(resultp1,resultp2,dstptr); 238 return(NOEXCEPTION); 239 } 240 } 241 else 242 { 243 Dbl_addition(leftp1,leftp2,rightp1,rightp2, 244 /*into*/resultp1,resultp2); 245 if(Dbl_isone_hidden(resultp1)) 246 { 247 Dbl_copytoptr(resultp1,resultp2,dstptr); 248 return(NOEXCEPTION); 249 } 250 } 251 if(Is_underflowtrap_enabled()) 252 { 253 /* need to normalize result */ 254 sign_save = Dbl_signextendedsign(resultp1); 255 Dbl_leftshiftby1(resultp1,resultp2); 256 Dbl_normalize(resultp1,resultp2,result_exponent); 257 Dbl_set_sign(resultp1,/*using*/sign_save); 258 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl); 259 Dbl_copytoptr(resultp1,resultp2,dstptr); 260 /* inexact = FALSE */ 261 return(UNDERFLOWEXCEPTION); 262 } 263 Dbl_copytoptr(resultp1,resultp2,dstptr); 264 return(NOEXCEPTION); 265 } 266 right_exponent = 1; /* Set exponent to reflect different bias 267 * with denomalized numbers. */ 268 } 269 else 270 { 271 Dbl_clear_signexponent_set_hidden(rightp1); 272 } 273 Dbl_clear_exponent_set_hidden(leftp1); 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 > DBL_THRESHOLD) 283 { 284 diff_exponent = DBL_THRESHOLD; 285 } 286 287 /* Align right operand by shifting to right */ 288 Dbl_right_align(/*operand*/rightp1,rightp2,/*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 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2, 300 /*with*/extent,/*into*/resultp1,resultp2); 301 if(Dbl_iszero_hidden(resultp1)) 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 = Dbl_signextendedsign(resultp1); 314 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2); 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(Dbl_iszero(resultp1,resultp2)) 321 /* Must have been "x-x" or "x+(-x)". */ 322 { 323 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1); 324 Dbl_copytoptr(resultp1,resultp2,dstptr); 325 return(NOEXCEPTION); 326 } 327 result_exponent--; 328 /* Look to see if normalization is finished. */ 329 if(Dbl_isone_hidden(resultp1)) 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 Dbl_set_sign(resultp1,/*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(Dbl_iszero_hiddenhigh7mantissa(resultp1)) 361 { 362 Dbl_leftshiftby8(resultp1,resultp2); 363 if((result_exponent -= 8) <= 0 && !underflowtrap) 364 goto underflow; 365 } 366 /* Now narrow it down to the nibble */ 367 if(Dbl_iszero_hiddenhigh3mantissa(resultp1)) 368 { 369 /* The lower nibble contains the normalizing one */ 370 Dbl_leftshiftby4(resultp1,resultp2); 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 = Dbl_hiddenhigh3mantissa(resultp1)) > 7) 377 { 378 /* Already normalized */ 379 if(result_exponent <= 0) goto underflow; 380 Dbl_set_sign(resultp1,/*using*/sign_save); 381 Dbl_set_exponent(resultp1,/*using*/result_exponent); 382 Dbl_copytoptr(resultp1,resultp2,dstptr); 383 return(NOEXCEPTION); 384 } 385 Dbl_sethigh4bits(resultp1,/*using*/sign_save); 386 switch(jumpsize) 387 { 388 case 1: 389 { 390 Dbl_leftshiftby3(resultp1,resultp2); 391 result_exponent -= 3; 392 break; 393 } 394 case 2: 395 case 3: 396 { 397 Dbl_leftshiftby2(resultp1,resultp2); 398 result_exponent -= 2; 399 break; 400 } 401 case 4: 402 case 5: 403 case 6: 404 case 7: 405 { 406 Dbl_leftshiftby1(resultp1,resultp2); 407 result_exponent -= 1; 408 break; 409 } 410 } 411 if(result_exponent > 0) 412 { 413 Dbl_set_exponent(resultp1,/*using*/result_exponent); 414 Dbl_copytoptr(resultp1,resultp2,dstptr); 415 return(NOEXCEPTION); /* Sign bit is already set */ 416 } 417 /* Fixup potential underflows */ 418 underflow: 419 if(Is_underflowtrap_enabled()) 420 { 421 Dbl_set_sign(resultp1,sign_save); 422 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl); 423 Dbl_copytoptr(resultp1,resultp2,dstptr); 424 /* inexact = FALSE */ 425 return(UNDERFLOWEXCEPTION); 426 } 427 /* 428 * Since we cannot get an inexact denormalized result, 429 * we can now return. 430 */ 431 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent); 432 Dbl_clear_signexponent(resultp1); 433 Dbl_set_sign(resultp1,sign_save); 434 Dbl_copytoptr(resultp1,resultp2,dstptr); 435 return(NOEXCEPTION); 436 } /* end if(hidden...)... */ 437 /* Fall through and round */ 438 } /* end if(save < 0)... */ 439 else 440 { 441 /* Add magnitudes */ 442 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2); 443 if(Dbl_isone_hiddenoverflow(resultp1)) 444 { 445 /* Prenormalization required. */ 446 Dbl_rightshiftby1_withextent(resultp2,extent,extent); 447 Dbl_arithrightshiftby1(resultp1,resultp2); 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 Dbl_isone_lowmantissap2(resultp2)) 468 { 469 /* either exactly half way and odd or more than 1/2ulp */ 470 Dbl_increment(resultp1,resultp2); 471 } 472 } 473 break; 474 475 case ROUNDPLUS: 476 if(Dbl_iszero_sign(resultp1)) 477 { 478 /* Round up positive results */ 479 Dbl_increment(resultp1,resultp2); 480 } 481 break; 482 483 case ROUNDMINUS: 484 if(Dbl_isone_sign(resultp1)) 485 { 486 /* Round down negative results */ 487 Dbl_increment(resultp1,resultp2); 488 } 489 490 case ROUNDZERO:; 491 /* truncate is simple */ 492 } /* end switch... */ 493 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++; 494 } 495 if(result_exponent == DBL_INFINITY_EXPONENT) 496 { 497 /* Overflow */ 498 if(Is_overflowtrap_enabled()) 499 { 500 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl); 501 Dbl_copytoptr(resultp1,resultp2,dstptr); 502 if (inexact) { 503 if (Is_inexacttrap_enabled()) 504 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 505 else 506 Set_inexactflag(); 507 } 508 return(OVERFLOWEXCEPTION); 509 } 510 else 511 { 512 inexact = TRUE; 513 Set_overflowflag(); 514 Dbl_setoverflow(resultp1,resultp2); 515 } 516 } 517 else Dbl_set_exponent(resultp1,result_exponent); 518 Dbl_copytoptr(resultp1,resultp2,dstptr); 519 if(inexact) { 520 if(Is_inexacttrap_enabled()) 521 return(INEXACTEXCEPTION); 522 else 523 Set_inexactflag(); 524 } 525 return(NOEXCEPTION); 526 } 527