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 } 272