![[BACK]](/icons/back.gif){kind=icon}
Return to dfrem.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath|
Return to dfrem.c CVS log | Up to [local] / sys / arch / hppa / spmath |
1.1 ! nbrk 1: /* $OpenBSD: dfrem.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: /* @(#)dfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:05:43 */
! 16:
! 17: #include "float.h"
! 18: #include "dbl_float.h"
! 19:
! 20: /*
! 21: * Double Precision Floating-point Remainder
! 22: */
! 23: int
! 24: dbl_frem(srcptr1,srcptr2,dstptr,status)
! 25: dbl_floating_point *srcptr1, *srcptr2, *dstptr;
! 26: unsigned int *status;
! 27: {
! 28: register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
! 29: register unsigned int resultp1, resultp2;
! 30: register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
! 31: register int roundup = FALSE;
! 32:
! 33: Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
! 34: Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
! 35: /*
! 36: * check first operand for NaN's or infinity
! 37: */
! 38: if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
! 39: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
! 40: if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
! 41: /* invalid since first operand is infinity */
! 42: if (Is_invalidtrap_enabled())
! 43: return(INVALIDEXCEPTION);
! 44: Set_invalidflag();
! 45: Dbl_makequietnan(resultp1,resultp2);
! 46: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 47: return(NOEXCEPTION);
! 48: }
! 49: }
! 50: else {
! 51: /*
! 52: * is NaN; signaling or quiet?
! 53: */
! 54: if (Dbl_isone_signaling(opnd1p1)) {
! 55: /* trap if INVALIDTRAP enabled */
! 56: if (Is_invalidtrap_enabled())
! 57: return(INVALIDEXCEPTION);
! 58: /* make NaN quiet */
! 59: Set_invalidflag();
! 60: Dbl_set_quiet(opnd1p1);
! 61: }
! 62: /*
! 63: * is second operand a signaling NaN?
! 64: */
! 65: else if (Dbl_is_signalingnan(opnd2p1)) {
! 66: /* trap if INVALIDTRAP enabled */
! 67: if (Is_invalidtrap_enabled())
! 68: return(INVALIDEXCEPTION);
! 69: /* make NaN quiet */
! 70: Set_invalidflag();
! 71: Dbl_set_quiet(opnd2p1);
! 72: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
! 73: return(NOEXCEPTION);
! 74: }
! 75: /*
! 76: * return quiet NaN
! 77: */
! 78: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
! 79: return(NOEXCEPTION);
! 80: }
! 81: }
! 82: /*
! 83: * check second operand for NaN's or infinity
! 84: */
! 85: if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
! 86: if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
! 87: /*
! 88: * return first operand
! 89: */
! 90: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
! 91: return(NOEXCEPTION);
! 92: }
! 93: /*
! 94: * is NaN; signaling or quiet?
! 95: */
! 96: if (Dbl_isone_signaling(opnd2p1)) {
! 97: /* trap if INVALIDTRAP enabled */
! 98: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 99: /* make NaN quiet */
! 100: Set_invalidflag();
! 101: Dbl_set_quiet(opnd2p1);
! 102: }
! 103: /*
! 104: * return quiet NaN
! 105: */
! 106: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
! 107: return(NOEXCEPTION);
! 108: }
! 109: /*
! 110: * check second operand for zero
! 111: */
! 112: if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
! 113: /* invalid since second operand is zero */
! 114: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
! 115: Set_invalidflag();
! 116: Dbl_makequietnan(resultp1,resultp2);
! 117: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 118: return(NOEXCEPTION);
! 119: }
! 120:
! 121: /*
! 122: * get sign of result
! 123: */
! 124: resultp1 = opnd1p1;
! 125:
! 126: /*
! 127: * check for denormalized operands
! 128: */
! 129: if (opnd1_exponent == 0) {
! 130: /* check for zero */
! 131: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
! 132: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
! 133: return(NOEXCEPTION);
! 134: }
! 135: /* normalize, then continue */
! 136: opnd1_exponent = 1;
! 137: Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
! 138: }
! 139: else {
! 140: Dbl_clear_signexponent_set_hidden(opnd1p1);
! 141: }
! 142: if (opnd2_exponent == 0) {
! 143: /* normalize, then continue */
! 144: opnd2_exponent = 1;
! 145: Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
! 146: }
! 147: else {
! 148: Dbl_clear_signexponent_set_hidden(opnd2p1);
! 149: }
! 150:
! 151: /* find result exponent and divide step loop count */
! 152: dest_exponent = opnd2_exponent - 1;
! 153: stepcount = opnd1_exponent - opnd2_exponent;
! 154:
! 155: /*
! 156: * check for opnd1/opnd2 < 1
! 157: */
! 158: if (stepcount < 0) {
! 159: /*
! 160: * check for opnd1/opnd2 > 1/2
! 161: *
! 162: * In this case n will round to 1, so
! 163: * r = opnd1 - opnd2
! 164: */
! 165: if (stepcount == -1 &&
! 166: Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
! 167: /* set sign */
! 168: Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
! 169: /* align opnd2 with opnd1 */
! 170: Dbl_leftshiftby1(opnd2p1,opnd2p2);
! 171: Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
! 172: opnd2p1,opnd2p2);
! 173: /* now normalize */
! 174: while (Dbl_iszero_hidden(opnd2p1)) {
! 175: Dbl_leftshiftby1(opnd2p1,opnd2p2);
! 176: dest_exponent--;
! 177: }
! 178: Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
! 179: goto testforunderflow;
! 180: }
! 181: /*
! 182: * opnd1/opnd2 <= 1/2
! 183: *
! 184: * In this case n will round to zero, so
! 185: * r = opnd1
! 186: */
! 187: Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
! 188: dest_exponent = opnd1_exponent;
! 189: goto testforunderflow;
! 190: }
! 191:
! 192: /*
! 193: * Generate result
! 194: *
! 195: * Do iterative subtract until remainder is less than operand 2.
! 196: */
! 197: while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
! 198: if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
! 199: Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
! 200: }
! 201: Dbl_leftshiftby1(opnd1p1,opnd1p2);
! 202: }
! 203: /*
! 204: * Do last subtract, then determine which way to round if remainder
! 205: * is exactly 1/2 of opnd2
! 206: */
! 207: if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
! 208: Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
! 209: roundup = TRUE;
! 210: }
! 211: if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
! 212: /* division is exact, remainder is zero */
! 213: Dbl_setzero_exponentmantissa(resultp1,resultp2);
! 214: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 215: return(NOEXCEPTION);
! 216: }
! 217:
! 218: /*
! 219: * Check for cases where opnd1/opnd2 < n
! 220: *
! 221: * In this case the result's sign will be opposite that of
! 222: * opnd1. The mantissa also needs some correction.
! 223: */
! 224: Dbl_leftshiftby1(opnd1p1,opnd1p2);
! 225: if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
! 226: Dbl_invert_sign(resultp1);
! 227: Dbl_leftshiftby1(opnd2p1,opnd2p2);
! 228: Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
! 229: }
! 230: /* check for remainder being exactly 1/2 of opnd2 */
! 231: else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
! 232: Dbl_invert_sign(resultp1);
! 233: }
! 234:
! 235: /* normalize result's mantissa */
! 236: while (Dbl_iszero_hidden(opnd1p1)) {
! 237: dest_exponent--;
! 238: Dbl_leftshiftby1(opnd1p1,opnd1p2);
! 239: }
! 240: Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
! 241:
! 242: /*
! 243: * Test for underflow
! 244: */
! 245: testforunderflow:
! 246: if (dest_exponent <= 0) {
! 247: /* trap if UNDERFLOWTRAP enabled */
! 248: if (Is_underflowtrap_enabled()) {
! 249: /*
! 250: * Adjust bias of result
! 251: */
! 252: Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
! 253: /* frem is always exact */
! 254: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 255: return(UNDERFLOWEXCEPTION);
! 256: }
! 257: /*
! 258: * denormalize result or set to signed zero
! 259: */
! 260: if (dest_exponent >= (1 - DBL_P)) {
! 261: Dbl_rightshift_exponentmantissa(resultp1,resultp2,
! 262: 1-dest_exponent);
! 263: }
! 264: else {
! 265: Dbl_setzero_exponentmantissa(resultp1,resultp2);
! 266: }
! 267: }
! 268: else Dbl_set_exponent(resultp1,dest_exponent);
! 269: Dbl_copytoptr(resultp1,resultp2,dstptr);
! 270: return(NOEXCEPTION);
! 271: }