Annotation of sys/arch/hppa/spmath/dfadd.c, Revision 1.1
1.1 ! nbrk 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: }
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to dfadd.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath