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