Annotation of sys/arch/hppa/spmath/sfsub.c, Revision 1.1
1.1 ! nbrk 1: /* $OpenBSD: sfsub.c,v 1.6 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: /* @(#)sfsub.c: Revision: 2.7.88.1 Date: 93/12/07 15:07:15 */
! 16:
! 17: #include "float.h"
! 18: #include "sgl_float.h"
! 19:
! 20: /*
! 21: * Single_subtract: subtract two single precision values.
! 22: */
! 23: int
! 24: sgl_fsub(leftptr, rightptr, dstptr, status)
! 25: sgl_floating_point *leftptr, *rightptr, *dstptr;
! 26: unsigned int *status;
! 27: {
! 28: register unsigned int left, right, result, extent;
! 29: register unsigned int signless_upper_left, signless_upper_right, save;
! 30:
! 31: register int result_exponent, right_exponent, diff_exponent;
! 32: register int sign_save, jumpsize;
! 33: register int inexact = FALSE, underflowtrap;
! 34:
! 35: /* Create local copies of the numbers */
! 36: left = *leftptr;
! 37: right = *rightptr;
! 38:
! 39: /* A zero "save" helps discover equal operands (for later), *
! 40: * and is used in swapping operands (if needed). */
! 41: Sgl_xortointp1(left,right,/*to*/save);
! 42:
! 43: /*
! 44: * check first operand for NaN's or infinity
! 45: */
! 46: if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT)
! 47: {
! 48: if (Sgl_iszero_mantissa(left))
! 49: {
! 50: if (Sgl_isnotnan(right))
! 51: {
! 52: if (Sgl_isinfinity(right) && save==0)
! 53: {
! 54: /*
! 55: * invalid since operands are same signed infinity's
! 56: */
! 57: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 58: Set_invalidflag();
! 59: Sgl_makequietnan(result);
! 60: *dstptr = result;
! 61: return(NOEXCEPTION);
! 62: }
! 63: /*
! 64: * return infinity
! 65: */
! 66: *dstptr = left;
! 67: return(NOEXCEPTION);
! 68: }
! 69: }
! 70: else
! 71: {
! 72: /*
! 73: * is NaN; signaling or quiet?
! 74: */
! 75: if (Sgl_isone_signaling(left))
! 76: {
! 77: /* trap if INVALIDTRAP enabled */
! 78: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 79: /* make NaN quiet */
! 80: Set_invalidflag();
! 81: Sgl_set_quiet(left);
! 82: }
! 83: /*
! 84: * is second operand a signaling NaN?
! 85: */
! 86: else if (Sgl_is_signalingnan(right))
! 87: {
! 88: /* trap if INVALIDTRAP enabled */
! 89: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 90: /* make NaN quiet */
! 91: Set_invalidflag();
! 92: Sgl_set_quiet(right);
! 93: *dstptr = right;
! 94: return(NOEXCEPTION);
! 95: }
! 96: /*
! 97: * return quiet NaN
! 98: */
! 99: *dstptr = left;
! 100: return(NOEXCEPTION);
! 101: }
! 102: } /* End left NaN or Infinity processing */
! 103: /*
! 104: * check second operand for NaN's or infinity
! 105: */
! 106: if (Sgl_isinfinity_exponent(right))
! 107: {
! 108: if (Sgl_iszero_mantissa(right))
! 109: {
! 110: /* return infinity */
! 111: Sgl_invert_sign(right);
! 112: *dstptr = right;
! 113: return(NOEXCEPTION);
! 114: }
! 115: /*
! 116: * is NaN; signaling or quiet?
! 117: */
! 118: if (Sgl_isone_signaling(right))
! 119: {
! 120: /* trap if INVALIDTRAP enabled */
! 121: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 122: /* make NaN quiet */
! 123: Set_invalidflag();
! 124: Sgl_set_quiet(right);
! 125: }
! 126: /*
! 127: * return quiet NaN
! 128: */
! 129: *dstptr = right;
! 130: return(NOEXCEPTION);
! 131: } /* End right NaN or Infinity processing */
! 132:
! 133: /* Invariant: Must be dealing with finite numbers */
! 134:
! 135: /* Compare operands by removing the sign */
! 136: Sgl_copytoint_exponentmantissa(left,signless_upper_left);
! 137: Sgl_copytoint_exponentmantissa(right,signless_upper_right);
! 138:
! 139: /* sign difference selects sub or add operation. */
! 140: if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right))
! 141: {
! 142: /* Set the left operand to the larger one by XOR swap *
! 143: * First finish the first word using "save" */
! 144: Sgl_xorfromintp1(save,right,/*to*/right);
! 145: Sgl_xorfromintp1(save,left,/*to*/left);
! 146: result_exponent = Sgl_exponent(left);
! 147: Sgl_invert_sign(left);
! 148: }
! 149: /* Invariant: left is not smaller than right. */
! 150:
! 151: if((right_exponent = Sgl_exponent(right)) == 0)
! 152: {
! 153: /* Denormalized operands. First look for zeroes */
! 154: if(Sgl_iszero_mantissa(right))
! 155: {
! 156: /* right is zero */
! 157: if(Sgl_iszero_exponentmantissa(left))
! 158: {
! 159: /* Both operands are zeros */
! 160: Sgl_invert_sign(right);
! 161: if(Is_rounding_mode(ROUNDMINUS))
! 162: {
! 163: Sgl_or_signs(left,/*with*/right);
! 164: }
! 165: else
! 166: {
! 167: Sgl_and_signs(left,/*with*/right);
! 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 = Sgl_signextendedsign(left);
! 179: Sgl_leftshiftby1(left);
! 180: Sgl_normalize(left,result_exponent);
! 181: Sgl_set_sign(left,/*using*/sign_save);
! 182: Sgl_setwrapped_exponent(left,result_exponent,unfl);
! 183: *dstptr = left;
! 184: /* inexact = FALSE */
! 185: return(UNDERFLOWEXCEPTION);
! 186: }
! 187: }
! 188: *dstptr = left;
! 189: return(NOEXCEPTION);
! 190: }
! 191:
! 192: /* Neither are zeroes */
! 193: Sgl_clear_sign(right); /* 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: Sgl_subtract(left,/*minus*/right,/*into*/result);
! 202: if(Sgl_iszero_mantissa(result))
! 203: {
! 204: if(Is_rounding_mode(ROUNDMINUS))
! 205: {
! 206: Sgl_setone_sign(result);
! 207: }
! 208: else
! 209: {
! 210: Sgl_setzero_sign(result);
! 211: }
! 212: *dstptr = result;
! 213: return(NOEXCEPTION);
! 214: }
! 215: }
! 216: else
! 217: {
! 218: Sgl_addition(left,right,/*into*/result);
! 219: if(Sgl_isone_hidden(result))
! 220: {
! 221: *dstptr = result;
! 222: return(NOEXCEPTION);
! 223: }
! 224: }
! 225: if(Is_underflowtrap_enabled())
! 226: {
! 227: /* need to normalize result */
! 228: sign_save = Sgl_signextendedsign(result);
! 229: Sgl_leftshiftby1(result);
! 230: Sgl_normalize(result,result_exponent);
! 231: Sgl_set_sign(result,/*using*/sign_save);
! 232: Sgl_setwrapped_exponent(result,result_exponent,unfl);
! 233: *dstptr = result;
! 234: /* inexact = FALSE */
! 235: return(UNDERFLOWEXCEPTION);
! 236: }
! 237: *dstptr = result;
! 238: return(NOEXCEPTION);
! 239: }
! 240: right_exponent = 1; /* Set exponent to reflect different bias
! 241: * with denomalized numbers. */
! 242: }
! 243: else
! 244: {
! 245: Sgl_clear_signexponent_set_hidden(right);
! 246: }
! 247: Sgl_clear_exponent_set_hidden(left);
! 248: diff_exponent = result_exponent - right_exponent;
! 249:
! 250: /*
! 251: * Special case alignment of operands that would force alignment
! 252: * beyond the extent of the extension. A further optimization
! 253: * could special case this but only reduces the path length for this
! 254: * infrequent case.
! 255: */
! 256: if(diff_exponent > SGL_THRESHOLD)
! 257: {
! 258: diff_exponent = SGL_THRESHOLD;
! 259: }
! 260:
! 261: /* Align right operand by shifting to right */
! 262: Sgl_right_align(/*operand*/right,/*shifted by*/diff_exponent,
! 263: /*and lower to*/extent);
! 264:
! 265: /* Treat sum and difference of the operands separately. */
! 266: if( (/*signed*/int) save >= 0 )
! 267: {
! 268: /*
! 269: * Difference of the two operands. Their can be no overflow. A
! 270: * borrow can occur out of the hidden bit and force a post
! 271: * normalization phase.
! 272: */
! 273: Sgl_subtract_withextension(left,/*minus*/right,/*with*/extent,/*into*/result);
! 274: if(Sgl_iszero_hidden(result))
! 275: {
! 276: /* Handle normalization */
! 277: /* A straight forward algorithm would now shift the result
! 278: * and extension left until the hidden bit becomes one. Not
! 279: * all of the extension bits need participate in the shift.
! 280: * Only the two most significant bits (round and guard) are
! 281: * needed. If only a single shift is needed then the guard
! 282: * bit becomes a significant low order bit and the extension
! 283: * must participate in the rounding. If more than a single
! 284: * shift is needed, then all bits to the right of the guard
! 285: * bit are zeros, and the guard bit may or may not be zero. */
! 286: sign_save = Sgl_signextendedsign(result);
! 287: Sgl_leftshiftby1_withextent(result,extent,result);
! 288:
! 289: /* Need to check for a zero result. The sign and exponent
! 290: * fields have already been zeroed. The more efficient test
! 291: * of the full object can be used.
! 292: */
! 293: if(Sgl_iszero(result))
! 294: /* Must have been "x-x" or "x+(-x)". */
! 295: {
! 296: if(Is_rounding_mode(ROUNDMINUS)) Sgl_setone_sign(result);
! 297: *dstptr = result;
! 298: return(NOEXCEPTION);
! 299: }
! 300: result_exponent--;
! 301: /* Look to see if normalization is finished. */
! 302: if(Sgl_isone_hidden(result))
! 303: {
! 304: if(result_exponent==0)
! 305: {
! 306: /* Denormalized, exponent should be zero. Left operand *
! 307: * was normalized, so extent (guard, round) was zero */
! 308: goto underflow;
! 309: }
! 310: else
! 311: {
! 312: /* No further normalization is needed. */
! 313: Sgl_set_sign(result,/*using*/sign_save);
! 314: Ext_leftshiftby1(extent);
! 315: goto round;
! 316: }
! 317: }
! 318:
! 319: /* Check for denormalized, exponent should be zero. Left *
! 320: * operand was normalized, so extent (guard, round) was zero */
! 321: if(!(underflowtrap = Is_underflowtrap_enabled()) &&
! 322: result_exponent==0) goto underflow;
! 323:
! 324: /* Shift extension to complete one bit of normalization and
! 325: * update exponent. */
! 326: Ext_leftshiftby1(extent);
! 327:
! 328: /* Discover first one bit to determine shift amount. Use a
! 329: * modified binary search. We have already shifted the result
! 330: * one position right and still not found a one so the remainder
! 331: * of the extension must be zero and simplifies rounding. */
! 332: /* Scan bytes */
! 333: while(Sgl_iszero_hiddenhigh7mantissa(result))
! 334: {
! 335: Sgl_leftshiftby8(result);
! 336: if((result_exponent -= 8) <= 0 && !underflowtrap)
! 337: goto underflow;
! 338: }
! 339: /* Now narrow it down to the nibble */
! 340: if(Sgl_iszero_hiddenhigh3mantissa(result))
! 341: {
! 342: /* The lower nibble contains the normalizing one */
! 343: Sgl_leftshiftby4(result);
! 344: if((result_exponent -= 4) <= 0 && !underflowtrap)
! 345: goto underflow;
! 346: }
! 347: /* Select case were first bit is set (already normalized)
! 348: * otherwise select the proper shift. */
! 349: if((jumpsize = Sgl_hiddenhigh3mantissa(result)) > 7)
! 350: {
! 351: /* Already normalized */
! 352: if(result_exponent <= 0) goto underflow;
! 353: Sgl_set_sign(result,/*using*/sign_save);
! 354: Sgl_set_exponent(result,/*using*/result_exponent);
! 355: *dstptr = result;
! 356: return(NOEXCEPTION);
! 357: }
! 358: Sgl_sethigh4bits(result,/*using*/sign_save);
! 359: switch(jumpsize)
! 360: {
! 361: case 1:
! 362: {
! 363: Sgl_leftshiftby3(result);
! 364: result_exponent -= 3;
! 365: break;
! 366: }
! 367: case 2:
! 368: case 3:
! 369: {
! 370: Sgl_leftshiftby2(result);
! 371: result_exponent -= 2;
! 372: break;
! 373: }
! 374: case 4:
! 375: case 5:
! 376: case 6:
! 377: case 7:
! 378: {
! 379: Sgl_leftshiftby1(result);
! 380: result_exponent -= 1;
! 381: break;
! 382: }
! 383: }
! 384: if(result_exponent > 0)
! 385: {
! 386: Sgl_set_exponent(result,/*using*/result_exponent);
! 387: *dstptr = result; /* Sign bit is already set */
! 388: return(NOEXCEPTION);
! 389: }
! 390: /* Fixup potential underflows */
! 391: underflow:
! 392: if(Is_underflowtrap_enabled())
! 393: {
! 394: Sgl_set_sign(result,sign_save);
! 395: Sgl_setwrapped_exponent(result,result_exponent,unfl);
! 396: *dstptr = result;
! 397: /* inexact = FALSE */
! 398: return(UNDERFLOWEXCEPTION);
! 399: }
! 400: /*
! 401: * Since we cannot get an inexact denormalized result,
! 402: * we can now return.
! 403: */
! 404: Sgl_right_align(result,/*by*/(1-result_exponent),extent);
! 405: Sgl_clear_signexponent(result);
! 406: Sgl_set_sign(result,sign_save);
! 407: *dstptr = result;
! 408: return(NOEXCEPTION);
! 409: } /* end if(hidden...)... */
! 410: /* Fall through and round */
! 411: } /* end if(save >= 0)... */
! 412: else
! 413: {
! 414: /* Add magnitudes */
! 415: Sgl_addition(left,right,/*to*/result);
! 416: if(Sgl_isone_hiddenoverflow(result))
! 417: {
! 418: /* Prenormalization required. */
! 419: Sgl_rightshiftby1_withextent(result,extent,extent);
! 420: Sgl_arithrightshiftby1(result);
! 421: result_exponent++;
! 422: } /* end if hiddenoverflow... */
! 423: } /* end else ...sub magnitudes... */
! 424:
! 425: /* Round the result. If the extension is all zeros,then the result is
! 426: * exact. Otherwise round in the correct direction. No underflow is
! 427: * possible. If a postnormalization is necessary, then the mantissa is
! 428: * all zeros so no shift is needed. */
! 429: round:
! 430: if(Ext_isnotzero(extent))
! 431: {
! 432: inexact = TRUE;
! 433: switch(Rounding_mode())
! 434: {
! 435: case ROUNDNEAREST: /* The default. */
! 436: if(Ext_isone_sign(extent))
! 437: {
! 438: /* at least 1/2 ulp */
! 439: if(Ext_isnotzero_lower(extent) ||
! 440: Sgl_isone_lowmantissa(result))
! 441: {
! 442: /* either exactly half way and odd or more than 1/2ulp */
! 443: Sgl_increment(result);
! 444: }
! 445: }
! 446: break;
! 447:
! 448: case ROUNDPLUS:
! 449: if(Sgl_iszero_sign(result))
! 450: {
! 451: /* Round up positive results */
! 452: Sgl_increment(result);
! 453: }
! 454: break;
! 455:
! 456: case ROUNDMINUS:
! 457: if(Sgl_isone_sign(result))
! 458: {
! 459: /* Round down negative results */
! 460: Sgl_increment(result);
! 461: }
! 462:
! 463: case ROUNDZERO:;
! 464: /* truncate is simple */
! 465: } /* end switch... */
! 466: if(Sgl_isone_hiddenoverflow(result)) result_exponent++;
! 467: }
! 468: if(result_exponent == SGL_INFINITY_EXPONENT)
! 469: {
! 470: /* Overflow */
! 471: if(Is_overflowtrap_enabled())
! 472: {
! 473: Sgl_setwrapped_exponent(result,result_exponent,ovfl);
! 474: *dstptr = result;
! 475: if (inexact) {
! 476: if (Is_inexacttrap_enabled())
! 477: return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
! 478: else Set_inexactflag();
! 479: }
! 480: return(OVERFLOWEXCEPTION);
! 481: }
! 482: else
! 483: {
! 484: Set_overflowflag();
! 485: inexact = TRUE;
! 486: Sgl_setoverflow(result);
! 487: }
! 488: }
! 489: else Sgl_set_exponent(result,result_exponent);
! 490: *dstptr = result;
! 491: if(inexact) {
! 492: if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
! 493: else Set_inexactflag();
! 494: }
! 495: return(NOEXCEPTION);
! 496: }
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to sfsub.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath