Annotation of sys/arch/hppa/spmath/sfmpy.c, Revision 1.1
1.1 ! nbrk 1: /* $OpenBSD: sfmpy.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: /* @(#)sfmpy.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:07 */
! 16:
! 17: #include "float.h"
! 18: #include "sgl_float.h"
! 19:
! 20: /*
! 21: * Single Precision Floating-point Multiply
! 22: */
! 23: int
! 24: sgl_fmpy(srcptr1,srcptr2,dstptr,status)
! 25: sgl_floating_point *srcptr1, *srcptr2, *dstptr;
! 26: unsigned int *status;
! 27: {
! 28: register unsigned int opnd1, opnd2, opnd3, result;
! 29: register int dest_exponent, count;
! 30: register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
! 31: int is_tiny;
! 32:
! 33: opnd1 = *srcptr1;
! 34: opnd2 = *srcptr2;
! 35: /*
! 36: * set sign bit of result
! 37: */
! 38: if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
! 39: else Sgl_setzero(result);
! 40: /*
! 41: * check first operand for NaN's or infinity
! 42: */
! 43: if (Sgl_isinfinity_exponent(opnd1)) {
! 44: if (Sgl_iszero_mantissa(opnd1)) {
! 45: if (Sgl_isnotnan(opnd2)) {
! 46: if (Sgl_iszero_exponentmantissa(opnd2)) {
! 47: /*
! 48: * invalid since operands are infinity
! 49: * and zero
! 50: */
! 51: if (Is_invalidtrap_enabled())
! 52: return(INVALIDEXCEPTION);
! 53: Set_invalidflag();
! 54: Sgl_makequietnan(result);
! 55: *dstptr = result;
! 56: return(NOEXCEPTION);
! 57: }
! 58: /*
! 59: * return infinity
! 60: */
! 61: Sgl_setinfinity_exponentmantissa(result);
! 62: *dstptr = result;
! 63: return(NOEXCEPTION);
! 64: }
! 65: }
! 66: else {
! 67: /*
! 68: * is NaN; signaling or quiet?
! 69: */
! 70: if (Sgl_isone_signaling(opnd1)) {
! 71: /* trap if INVALIDTRAP enabled */
! 72: if (Is_invalidtrap_enabled())
! 73: return(INVALIDEXCEPTION);
! 74: /* make NaN quiet */
! 75: Set_invalidflag();
! 76: Sgl_set_quiet(opnd1);
! 77: }
! 78: /*
! 79: * is second operand a signaling NaN?
! 80: */
! 81: else if (Sgl_is_signalingnan(opnd2)) {
! 82: /* trap if INVALIDTRAP enabled */
! 83: if (Is_invalidtrap_enabled())
! 84: return(INVALIDEXCEPTION);
! 85: /* make NaN quiet */
! 86: Set_invalidflag();
! 87: Sgl_set_quiet(opnd2);
! 88: *dstptr = opnd2;
! 89: return(NOEXCEPTION);
! 90: }
! 91: /*
! 92: * return quiet NaN
! 93: */
! 94: *dstptr = opnd1;
! 95: return(NOEXCEPTION);
! 96: }
! 97: }
! 98: /*
! 99: * check second operand for NaN's or infinity
! 100: */
! 101: if (Sgl_isinfinity_exponent(opnd2)) {
! 102: if (Sgl_iszero_mantissa(opnd2)) {
! 103: if (Sgl_iszero_exponentmantissa(opnd1)) {
! 104: /* invalid since operands are zero & infinity */
! 105: if (Is_invalidtrap_enabled())
! 106: return(INVALIDEXCEPTION);
! 107: Set_invalidflag();
! 108: Sgl_makequietnan(opnd2);
! 109: *dstptr = opnd2;
! 110: return(NOEXCEPTION);
! 111: }
! 112: /*
! 113: * return infinity
! 114: */
! 115: Sgl_setinfinity_exponentmantissa(result);
! 116: *dstptr = result;
! 117: return(NOEXCEPTION);
! 118: }
! 119: /*
! 120: * is NaN; signaling or quiet?
! 121: */
! 122: if (Sgl_isone_signaling(opnd2)) {
! 123: /* trap if INVALIDTRAP enabled */
! 124: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 125:
! 126: /* make NaN quiet */
! 127: Set_invalidflag();
! 128: Sgl_set_quiet(opnd2);
! 129: }
! 130: /*
! 131: * return quiet NaN
! 132: */
! 133: *dstptr = opnd2;
! 134: return(NOEXCEPTION);
! 135: }
! 136: /*
! 137: * Generate exponent
! 138: */
! 139: dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
! 140:
! 141: /*
! 142: * Generate mantissa
! 143: */
! 144: if (Sgl_isnotzero_exponent(opnd1)) {
! 145: /* set hidden bit */
! 146: Sgl_clear_signexponent_set_hidden(opnd1);
! 147: }
! 148: else {
! 149: /* check for zero */
! 150: if (Sgl_iszero_mantissa(opnd1)) {
! 151: Sgl_setzero_exponentmantissa(result);
! 152: *dstptr = result;
! 153: return(NOEXCEPTION);
! 154: }
! 155: /* is denormalized, adjust exponent */
! 156: Sgl_clear_signexponent(opnd1);
! 157: Sgl_leftshiftby1(opnd1);
! 158: Sgl_normalize(opnd1,dest_exponent);
! 159: }
! 160: /* opnd2 needs to have hidden bit set with msb in hidden bit */
! 161: if (Sgl_isnotzero_exponent(opnd2)) {
! 162: Sgl_clear_signexponent_set_hidden(opnd2);
! 163: }
! 164: else {
! 165: /* check for zero */
! 166: if (Sgl_iszero_mantissa(opnd2)) {
! 167: Sgl_setzero_exponentmantissa(result);
! 168: *dstptr = result;
! 169: return(NOEXCEPTION);
! 170: }
! 171: /* is denormalized; want to normalize */
! 172: Sgl_clear_signexponent(opnd2);
! 173: Sgl_leftshiftby1(opnd2);
! 174: Sgl_normalize(opnd2,dest_exponent);
! 175: }
! 176:
! 177: /* Multiply two source mantissas together */
! 178:
! 179: Sgl_leftshiftby4(opnd2); /* make room for guard bits */
! 180: Sgl_setzero(opnd3);
! 181: /*
! 182: * Four bits at a time are inspected in each loop, and a
! 183: * simple shift and add multiply algorithm is used.
! 184: */
! 185: for (count=1;count<SGL_P;count+=4) {
! 186: stickybit |= Slow4(opnd3);
! 187: Sgl_rightshiftby4(opnd3);
! 188: if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
! 189: if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
! 190: if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
! 191: if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
! 192: Sgl_rightshiftby4(opnd1);
! 193: }
! 194: /* make sure result is left-justified */
! 195: if (Sgl_iszero_sign(opnd3)) {
! 196: Sgl_leftshiftby1(opnd3);
! 197: }
! 198: else {
! 199: /* result mantissa >= 2. */
! 200: dest_exponent++;
! 201: }
! 202: /* check for denormalized result */
! 203: while (Sgl_iszero_sign(opnd3)) {
! 204: Sgl_leftshiftby1(opnd3);
! 205: dest_exponent--;
! 206: }
! 207: /*
! 208: * check for guard, sticky and inexact bits
! 209: */
! 210: stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
! 211: guardbit = Sbit24(opnd3);
! 212: inexact = guardbit | stickybit;
! 213:
! 214: /* re-align mantissa */
! 215: Sgl_rightshiftby8(opnd3);
! 216:
! 217: /*
! 218: * round result
! 219: */
! 220: if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
! 221: Sgl_clear_signexponent(opnd3);
! 222: switch (Rounding_mode()) {
! 223: case ROUNDPLUS:
! 224: if (Sgl_iszero_sign(result))
! 225: Sgl_increment(opnd3);
! 226: break;
! 227: case ROUNDMINUS:
! 228: if (Sgl_isone_sign(result))
! 229: Sgl_increment(opnd3);
! 230: break;
! 231: case ROUNDNEAREST:
! 232: if (guardbit &&
! 233: (stickybit || Sgl_isone_lowmantissa(opnd3)))
! 234: Sgl_increment(opnd3);
! 235: break;
! 236: }
! 237: if (Sgl_isone_hidden(opnd3)) dest_exponent++;
! 238: }
! 239: Sgl_set_mantissa(result,opnd3);
! 240:
! 241: /*
! 242: * Test for overflow
! 243: */
! 244: if (dest_exponent >= SGL_INFINITY_EXPONENT) {
! 245: /* trap if OVERFLOWTRAP enabled */
! 246: if (Is_overflowtrap_enabled()) {
! 247: /*
! 248: * Adjust bias of result
! 249: */
! 250: Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
! 251: *dstptr = result;
! 252: if (inexact) {
! 253: if (Is_inexacttrap_enabled())
! 254: return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
! 255: else Set_inexactflag();
! 256: }
! 257: return(OVERFLOWEXCEPTION);
! 258: }
! 259: inexact = TRUE;
! 260: Set_overflowflag();
! 261: /* set result to infinity or largest number */
! 262: Sgl_setoverflow(result);
! 263: }
! 264: /*
! 265: * Test for underflow
! 266: */
! 267: else if (dest_exponent <= 0) {
! 268: /* trap if UNDERFLOWTRAP enabled */
! 269: if (Is_underflowtrap_enabled()) {
! 270: /*
! 271: * Adjust bias of result
! 272: */
! 273: Sgl_setwrapped_exponent(result,dest_exponent,unfl);
! 274: *dstptr = result;
! 275: if (inexact) {
! 276: if (Is_inexacttrap_enabled())
! 277: return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
! 278: else Set_inexactflag();
! 279: }
! 280: return(UNDERFLOWEXCEPTION);
! 281: }
! 282:
! 283: /* Determine if should set underflow flag */
! 284: is_tiny = TRUE;
! 285: if (dest_exponent == 0 && inexact) {
! 286: switch (Rounding_mode()) {
! 287: case ROUNDPLUS:
! 288: if (Sgl_iszero_sign(result)) {
! 289: Sgl_increment(opnd3);
! 290: if (Sgl_isone_hiddenoverflow(opnd3))
! 291: is_tiny = FALSE;
! 292: Sgl_decrement(opnd3);
! 293: }
! 294: break;
! 295: case ROUNDMINUS:
! 296: if (Sgl_isone_sign(result)) {
! 297: Sgl_increment(opnd3);
! 298: if (Sgl_isone_hiddenoverflow(opnd3))
! 299: is_tiny = FALSE;
! 300: Sgl_decrement(opnd3);
! 301: }
! 302: break;
! 303: case ROUNDNEAREST:
! 304: if (guardbit && (stickybit ||
! 305: Sgl_isone_lowmantissa(opnd3))) {
! 306: Sgl_increment(opnd3);
! 307: if (Sgl_isone_hiddenoverflow(opnd3))
! 308: is_tiny = FALSE;
! 309: Sgl_decrement(opnd3);
! 310: }
! 311: break;
! 312: }
! 313: }
! 314:
! 315: /*
! 316: * denormalize result or set to signed zero
! 317: */
! 318: stickybit = inexact;
! 319: Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
! 320:
! 321: /* return zero or smallest number */
! 322: if (inexact) {
! 323: switch (Rounding_mode()) {
! 324: case ROUNDPLUS:
! 325: if (Sgl_iszero_sign(result))
! 326: Sgl_increment(opnd3);
! 327: break;
! 328: case ROUNDMINUS:
! 329: if (Sgl_isone_sign(result))
! 330: Sgl_increment(opnd3);
! 331: break;
! 332: case ROUNDNEAREST:
! 333: if (guardbit && (stickybit ||
! 334: Sgl_isone_lowmantissa(opnd3)))
! 335: Sgl_increment(opnd3);
! 336: break;
! 337: }
! 338: if (is_tiny) Set_underflowflag();
! 339: }
! 340: Sgl_set_exponentmantissa(result,opnd3);
! 341: }
! 342: else Sgl_set_exponent(result,dest_exponent);
! 343: *dstptr = result;
! 344:
! 345: /* check for inexact */
! 346: if (inexact) {
! 347: if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
! 348: else Set_inexactflag();
! 349: }
! 350: return(NOEXCEPTION);
! 351: }
CVSweb