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