![[BACK]](/icons/back.gif){kind=icon}
Return to dfmpy.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath|
Return to dfmpy.c CVS log | Up to [local] / sys / arch / hppa / spmath |
1.1 ! nbrk 1: /* $OpenBSD: dfmpy.c,v 1.5 2002/05/07 22:19:30 mickey 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: /* @(#)dfmpy.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:41 */
! 16:
! 17: #include "float.h"
! 18: #include "dbl_float.h"
! 19:
! 20: /*
! 21: * Double Precision Floating-point Multiply
! 22: */
! 23:
! 24: int
! 25: dbl_fmpy(srcptr1,srcptr2,dstptr,status)
! 26: dbl_floating_point *srcptr1, *srcptr2, *dstptr;
! 27: unsigned int *status;
! 28: {
! 29: register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
! 30: register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
! 31: register int dest_exponent, count;
! 32: register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
! 33: int is_tiny;
! 34:
! 35: Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
! 36: Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
! 37:
! 38: /*
! 39: * set sign bit of result
! 40: */
! 41: if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
! 42: Dbl_setnegativezerop1(resultp1);
! 43: else Dbl_setzerop1(resultp1);
! 44: /*
! 45: * check first operand for NaN's or infinity
! 46: */
! 47: if (Dbl_isinfinity_exponent(opnd1p1)) {
! 48: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
! 49: if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
! 50: if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
! 51: /*
! 52: * invalid since operands are infinity
! 53: * and zero
! 54: */
! 55: if (Is_invalidtrap_enabled())
! 56: return(INVALIDEXCEPTION);
! 57: Set_invalidflag();
! 58: Dbl_makequietnan(resultp1,resultp2);
! 59: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 60: return(NOEXCEPTION);
! 61: }
! 62: /*
! 63: * return infinity
! 64: */
! 65: Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
! 66: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 67: return(NOEXCEPTION);
! 68: }
! 69: }
! 70: else {
! 71: /*
! 72: * is NaN; signaling or quiet?
! 73: */
! 74: if (Dbl_isone_signaling(opnd1p1)) {
! 75: /* trap if INVALIDTRAP enabled */
! 76: if (Is_invalidtrap_enabled())
! 77: return(INVALIDEXCEPTION);
! 78: /* make NaN quiet */
! 79: Set_invalidflag();
! 80: Dbl_set_quiet(opnd1p1);
! 81: }
! 82: /*
! 83: * is second operand a signaling NaN?
! 84: */
! 85: else if (Dbl_is_signalingnan(opnd2p1)) {
! 86: /* trap if INVALIDTRAP enabled */
! 87: if (Is_invalidtrap_enabled())
! 88: return(INVALIDEXCEPTION);
! 89: /* make NaN quiet */
! 90: Set_invalidflag();
! 91: Dbl_set_quiet(opnd2p1);
! 92: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
! 93: return(NOEXCEPTION);
! 94: }
! 95: /*
! 96: * return quiet NaN
! 97: */
! 98: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
! 99: return(NOEXCEPTION);
! 100: }
! 101: }
! 102: /*
! 103: * check second operand for NaN's or infinity
! 104: */
! 105: if (Dbl_isinfinity_exponent(opnd2p1)) {
! 106: if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
! 107: if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
! 108: /* invalid since operands are zero & infinity */
! 109: if (Is_invalidtrap_enabled())
! 110: return(INVALIDEXCEPTION);
! 111: Set_invalidflag();
! 112: Dbl_makequietnan(opnd2p1,opnd2p2);
! 113: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
! 114: return(NOEXCEPTION);
! 115: }
! 116: /*
! 117: * return infinity
! 118: */
! 119: Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
! 120: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 121: return(NOEXCEPTION);
! 122: }
! 123: /*
! 124: * is NaN; signaling or quiet?
! 125: */
! 126: if (Dbl_isone_signaling(opnd2p1)) {
! 127: /* trap if INVALIDTRAP enabled */
! 128: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 129: /* make NaN quiet */
! 130: Set_invalidflag();
! 131: Dbl_set_quiet(opnd2p1);
! 132: }
! 133: /*
! 134: * return quiet NaN
! 135: */
! 136: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
! 137: return(NOEXCEPTION);
! 138: }
! 139: /*
! 140: * Generate exponent
! 141: */
! 142: dest_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) -DBL_BIAS;
! 143:
! 144: /*
! 145: * Generate mantissa
! 146: */
! 147: if (Dbl_isnotzero_exponent(opnd1p1)) {
! 148: /* set hidden bit */
! 149: Dbl_clear_signexponent_set_hidden(opnd1p1);
! 150: }
! 151: else {
! 152: /* check for zero */
! 153: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
! 154: Dbl_setzero_exponentmantissa(resultp1,resultp2);
! 155: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 156: return(NOEXCEPTION);
! 157: }
! 158: /* is denormalized, adjust exponent */
! 159: Dbl_clear_signexponent(opnd1p1);
! 160: Dbl_leftshiftby1(opnd1p1,opnd1p2);
! 161: Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
! 162: }
! 163: /* opnd2 needs to have hidden bit set with msb in hidden bit */
! 164: if (Dbl_isnotzero_exponent(opnd2p1)) {
! 165: Dbl_clear_signexponent_set_hidden(opnd2p1);
! 166: }
! 167: else {
! 168: /* check for zero */
! 169: if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
! 170: Dbl_setzero_exponentmantissa(resultp1,resultp2);
! 171: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 172: return(NOEXCEPTION);
! 173: }
! 174: /* is denormalized; want to normalize */
! 175: Dbl_clear_signexponent(opnd2p1);
! 176: Dbl_leftshiftby1(opnd2p1,opnd2p2);
! 177: Dbl_normalize(opnd2p1,opnd2p2,dest_exponent);
! 178: }
! 179:
! 180: /* Multiply two source mantissas together */
! 181:
! 182: /* make room for guard bits */
! 183: Dbl_leftshiftby7(opnd2p1,opnd2p2);
! 184: Dbl_setzero(opnd3p1,opnd3p2);
! 185: /*
! 186: * Four bits at a time are inspected in each loop, and a
! 187: * simple shift and add multiply algorithm is used.
! 188: */
! 189: for (count=1;count<=DBL_P;count+=4) {
! 190: stickybit |= Dlow4p2(opnd3p2);
! 191: Dbl_rightshiftby4(opnd3p1,opnd3p2);
! 192: if (Dbit28p2(opnd1p2)) {
! 193: /* Twoword_add should be an ADDC followed by an ADD. */
! 194: Twoword_add(opnd3p1, opnd3p2, opnd2p1<<3 | opnd2p2>>29,
! 195: opnd2p2<<3);
! 196: }
! 197: if (Dbit29p2(opnd1p2)) {
! 198: Twoword_add(opnd3p1, opnd3p2, opnd2p1<<2 | opnd2p2>>30,
! 199: opnd2p2<<2);
! 200: }
! 201: if (Dbit30p2(opnd1p2)) {
! 202: Twoword_add(opnd3p1, opnd3p2, opnd2p1<<1 | opnd2p2>>31,
! 203: opnd2p2<<1);
! 204: }
! 205: if (Dbit31p2(opnd1p2)) {
! 206: Twoword_add(opnd3p1, opnd3p2, opnd2p1, opnd2p2);
! 207: }
! 208: Dbl_rightshiftby4(opnd1p1,opnd1p2);
! 209: }
! 210: if (Dbit3p1(opnd3p1)==0) {
! 211: Dbl_leftshiftby1(opnd3p1,opnd3p2);
! 212: }
! 213: else {
! 214: /* result mantissa >= 2. */
! 215: dest_exponent++;
! 216: }
! 217: /* check for denormalized result */
! 218: while (Dbit3p1(opnd3p1)==0) {
! 219: Dbl_leftshiftby1(opnd3p1,opnd3p2);
! 220: dest_exponent--;
! 221: }
! 222: /*
! 223: * check for guard, sticky and inexact bits
! 224: */
! 225: stickybit |= Dallp2(opnd3p2) << 25;
! 226: guardbit = (Dallp2(opnd3p2) << 24) >> 31;
! 227: inexact = guardbit | stickybit;
! 228:
! 229: /* align result mantissa */
! 230: Dbl_rightshiftby8(opnd3p1,opnd3p2);
! 231:
! 232: /*
! 233: * round result
! 234: */
! 235: if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
! 236: Dbl_clear_signexponent(opnd3p1);
! 237: switch (Rounding_mode()) {
! 238: case ROUNDPLUS:
! 239: if (Dbl_iszero_sign(resultp1))
! 240: Dbl_increment(opnd3p1,opnd3p2);
! 241: break;
! 242: case ROUNDMINUS:
! 243: if (Dbl_isone_sign(resultp1))
! 244: Dbl_increment(opnd3p1,opnd3p2);
! 245: break;
! 246: case ROUNDNEAREST:
! 247: if (guardbit &&
! 248: (stickybit || Dbl_isone_lowmantissap2(opnd3p2)))
! 249: Dbl_increment(opnd3p1,opnd3p2);
! 250: break;
! 251: }
! 252: if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
! 253: }
! 254: Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
! 255:
! 256: /*
! 257: * Test for overflow
! 258: */
! 259: if (dest_exponent >= DBL_INFINITY_EXPONENT) {
! 260: /* trap if OVERFLOWTRAP enabled */
! 261: if (Is_overflowtrap_enabled()) {
! 262: /*
! 263: * Adjust bias of result
! 264: */
! 265: Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
! 266: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 267: if (inexact) {
! 268: if (Is_inexacttrap_enabled())
! 269: return (OVERFLOWEXCEPTION | INEXACTEXCEPTION);
! 270: else
! 271: Set_inexactflag();
! 272: }
! 273: return (OVERFLOWEXCEPTION);
! 274: }
! 275: inexact = TRUE;
! 276: Set_overflowflag();
! 277: /* set result to infinity or largest number */
! 278: Dbl_setoverflow(resultp1,resultp2);
! 279: }
! 280: /*
! 281: * Test for underflow
! 282: */
! 283: else if (dest_exponent <= 0) {
! 284: /* trap if UNDERFLOWTRAP enabled */
! 285: if (Is_underflowtrap_enabled()) {
! 286: /*
! 287: * Adjust bias of result
! 288: */
! 289: Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
! 290: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 291: if (inexact) {
! 292: if (Is_inexacttrap_enabled())
! 293: return (UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
! 294: else
! 295: Set_inexactflag();
! 296: }
! 297: return (UNDERFLOWEXCEPTION);
! 298: }
! 299:
! 300: /* Determine if should set underflow flag */
! 301: is_tiny = TRUE;
! 302: if (dest_exponent == 0 && inexact) {
! 303: switch (Rounding_mode()) {
! 304: case ROUNDPLUS:
! 305: if (Dbl_iszero_sign(resultp1)) {
! 306: Dbl_increment(opnd3p1,opnd3p2);
! 307: if (Dbl_isone_hiddenoverflow(opnd3p1))
! 308: is_tiny = FALSE;
! 309: Dbl_decrement(opnd3p1,opnd3p2);
! 310: }
! 311: break;
! 312: case ROUNDMINUS:
! 313: if (Dbl_isone_sign(resultp1)) {
! 314: Dbl_increment(opnd3p1,opnd3p2);
! 315: if (Dbl_isone_hiddenoverflow(opnd3p1))
! 316: is_tiny = FALSE;
! 317: Dbl_decrement(opnd3p1,opnd3p2);
! 318: }
! 319: break;
! 320: case ROUNDNEAREST:
! 321: if (guardbit && (stickybit ||
! 322: Dbl_isone_lowmantissap2(opnd3p2))) {
! 323: Dbl_increment(opnd3p1,opnd3p2);
! 324: if (Dbl_isone_hiddenoverflow(opnd3p1))
! 325: is_tiny = FALSE;
! 326: Dbl_decrement(opnd3p1,opnd3p2);
! 327: }
! 328: break;
! 329: }
! 330: }
! 331:
! 332: /*
! 333: * denormalize result or set to signed zero
! 334: */
! 335: stickybit = inexact;
! 336: Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
! 337: stickybit,inexact);
! 338:
! 339: /* return zero or smallest number */
! 340: if (inexact) {
! 341: switch (Rounding_mode()) {
! 342: case ROUNDPLUS:
! 343: if (Dbl_iszero_sign(resultp1)) {
! 344: Dbl_increment(opnd3p1,opnd3p2);
! 345: }
! 346: break;
! 347: case ROUNDMINUS:
! 348: if (Dbl_isone_sign(resultp1)) {
! 349: Dbl_increment(opnd3p1,opnd3p2);
! 350: }
! 351: break;
! 352: case ROUNDNEAREST:
! 353: if (guardbit && (stickybit ||
! 354: Dbl_isone_lowmantissap2(opnd3p2))) {
! 355: Dbl_increment(opnd3p1,opnd3p2);
! 356: }
! 357: break;
! 358: }
! 359: if (is_tiny) Set_underflowflag();
! 360: }
! 361: Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
! 362: }
! 363: else Dbl_set_exponent(resultp1,dest_exponent);
! 364: /* check for inexact */
! 365: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 366: if (inexact) {
! 367: if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
! 368: else Set_inexactflag();
! 369: }
! 370: return(NOEXCEPTION);
! 371: }