1 /* $OpenBSD: dfrem.c,v 1.4 2001/03/29 03:58:17 mickey Exp $ */ 2 3 /* 4 * Copyright 1996 1995 by Open Software Foundation, Inc. 5 * All Rights Reserved 6 * 7 * Permission to use, copy, modify, and distribute this software and 8 * its documentation for any purpose and without fee is hereby granted, 9 * provided that the above copyright notice appears in all copies and 10 * that both the copyright notice and this permission notice appear in 11 * supporting documentation. 12 * 13 * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE 14 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 15 * FOR A PARTICULAR PURPOSE. 16 * 17 * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR 18 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM 19 * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT, 20 * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION 21 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 22 * 23 */ 24 /* 25 * pmk1.1 26 */ 27 /* 28 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 29 * 30 * To anyone who acknowledges that this file is provided "AS IS" 31 * without any express or implied warranty: 32 * permission to use, copy, modify, and distribute this file 33 * for any purpose is hereby granted without fee, provided that 34 * the above copyright notice and this notice appears in all 35 * copies, and that the name of Hewlett-Packard Company not be 36 * used in advertising or publicity pertaining to distribution 37 * of the software without specific, written prior permission. 38 * Hewlett-Packard Company makes no representations about the 39 * suitability of this software for any purpose. 40 */ 41 42 43 #include "../spmath/float.h" 44 #include "../spmath/dbl_float.h" 45 46 /* 47 * Double Precision Floating-point Remainder 48 */ 49 int 50 dbl_frem(srcptr1,srcptr2,dstptr,status) 51 52 dbl_floating_point *srcptr1, *srcptr2, *dstptr; 53 unsigned int *status; 54 { 55 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 56 register unsigned int resultp1, resultp2; 57 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount; 58 register int roundup = FALSE; 59 60 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 61 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 62 /* 63 * check first operand for NaN's or infinity 64 */ 65 if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) { 66 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 67 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 68 /* invalid since first operand is infinity */ 69 if (Is_invalidtrap_enabled()) 70 return(INVALIDEXCEPTION); 71 Set_invalidflag(); 72 Dbl_makequietnan(resultp1,resultp2); 73 Dbl_copytoptr(resultp1,resultp2,dstptr); 74 return(NOEXCEPTION); 75 } 76 } 77 else { 78 /* 79 * is NaN; signaling or quiet? 80 */ 81 if (Dbl_isone_signaling(opnd1p1)) { 82 /* trap if INVALIDTRAP enabled */ 83 if (Is_invalidtrap_enabled()) 84 return(INVALIDEXCEPTION); 85 /* make NaN quiet */ 86 Set_invalidflag(); 87 Dbl_set_quiet(opnd1p1); 88 } 89 /* 90 * is second operand a signaling NaN? 91 */ 92 else if (Dbl_is_signalingnan(opnd2p1)) { 93 /* trap if INVALIDTRAP enabled */ 94 if (Is_invalidtrap_enabled()) 95 return(INVALIDEXCEPTION); 96 /* make NaN quiet */ 97 Set_invalidflag(); 98 Dbl_set_quiet(opnd2p1); 99 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 100 return(NOEXCEPTION); 101 } 102 /* 103 * return quiet NaN 104 */ 105 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 106 return(NOEXCEPTION); 107 } 108 } 109 /* 110 * check second operand for NaN's or infinity 111 */ 112 if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) { 113 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 114 /* 115 * return first operand 116 */ 117 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 118 return(NOEXCEPTION); 119 } 120 /* 121 * is NaN; signaling or quiet? 122 */ 123 if (Dbl_isone_signaling(opnd2p1)) { 124 /* trap if INVALIDTRAP enabled */ 125 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 126 /* make NaN quiet */ 127 Set_invalidflag(); 128 Dbl_set_quiet(opnd2p1); 129 } 130 /* 131 * return quiet NaN 132 */ 133 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 134 return(NOEXCEPTION); 135 } 136 /* 137 * check second operand for zero 138 */ 139 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 140 /* invalid since second operand is zero */ 141 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 142 Set_invalidflag(); 143 Dbl_makequietnan(resultp1,resultp2); 144 Dbl_copytoptr(resultp1,resultp2,dstptr); 145 return(NOEXCEPTION); 146 } 147 148 /* 149 * get sign of result 150 */ 151 resultp1 = opnd1p1; 152 153 /* 154 * check for denormalized operands 155 */ 156 if (opnd1_exponent == 0) { 157 /* check for zero */ 158 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 159 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 160 return(NOEXCEPTION); 161 } 162 /* normalize, then continue */ 163 opnd1_exponent = 1; 164 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent); 165 } 166 else { 167 Dbl_clear_signexponent_set_hidden(opnd1p1); 168 } 169 if (opnd2_exponent == 0) { 170 /* normalize, then continue */ 171 opnd2_exponent = 1; 172 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent); 173 } 174 else { 175 Dbl_clear_signexponent_set_hidden(opnd2p1); 176 } 177 178 /* find result exponent and divide step loop count */ 179 dest_exponent = opnd2_exponent - 1; 180 stepcount = opnd1_exponent - opnd2_exponent; 181 182 /* 183 * check for opnd1/opnd2 < 1 184 */ 185 if (stepcount < 0) { 186 /* 187 * check for opnd1/opnd2 > 1/2 188 * 189 * In this case n will round to 1, so 190 * r = opnd1 - opnd2 191 */ 192 if (stepcount == -1 && 193 Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 194 /* set sign */ 195 Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1); 196 /* align opnd2 with opnd1 */ 197 Dbl_leftshiftby1(opnd2p1,opnd2p2); 198 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2, 199 opnd2p1,opnd2p2); 200 /* now normalize */ 201 while (Dbl_iszero_hidden(opnd2p1)) { 202 Dbl_leftshiftby1(opnd2p1,opnd2p2); 203 dest_exponent--; 204 } 205 Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2); 206 goto testforunderflow; 207 } 208 /* 209 * opnd1/opnd2 <= 1/2 210 * 211 * In this case n will round to zero, so 212 * r = opnd1 213 */ 214 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 215 dest_exponent = opnd1_exponent; 216 goto testforunderflow; 217 } 218 219 /* 220 * Generate result 221 * 222 * Do iterative subtract until remainder is less than operand 2. 223 */ 224 while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) { 225 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 226 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 227 } 228 Dbl_leftshiftby1(opnd1p1,opnd1p2); 229 } 230 /* 231 * Do last subtract, then determine which way to round if remainder 232 * is exactly 1/2 of opnd2 233 */ 234 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 235 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 236 roundup = TRUE; 237 } 238 if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) { 239 /* division is exact, remainder is zero */ 240 Dbl_setzero_exponentmantissa(resultp1,resultp2); 241 Dbl_copytoptr(resultp1,resultp2,dstptr); 242 return(NOEXCEPTION); 243 } 244 245 /* 246 * Check for cases where opnd1/opnd2 < n 247 * 248 * In this case the result's sign will be opposite that of 249 * opnd1. The mantissa also needs some correction. 250 */ 251 Dbl_leftshiftby1(opnd1p1,opnd1p2); 252 if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 253 Dbl_invert_sign(resultp1); 254 Dbl_leftshiftby1(opnd2p1,opnd2p2); 255 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2); 256 } 257 /* check for remainder being exactly 1/2 of opnd2 */ 258 else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) { 259 Dbl_invert_sign(resultp1); 260 } 261 262 /* normalize result's mantissa */ 263 while (Dbl_iszero_hidden(opnd1p1)) { 264 dest_exponent--; 265 Dbl_leftshiftby1(opnd1p1,opnd1p2); 266 } 267 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 268 269 /* 270 * Test for underflow 271 */ 272 testforunderflow: 273 if (dest_exponent <= 0) { 274 /* trap if UNDERFLOWTRAP enabled */ 275 if (Is_underflowtrap_enabled()) { 276 /* 277 * Adjust bias of result 278 */ 279 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 280 /* frem is always exact */ 281 Dbl_copytoptr(resultp1,resultp2,dstptr); 282 return(UNDERFLOWEXCEPTION); 283 } 284 /* 285 * denormalize result or set to signed zero 286 */ 287 if (dest_exponent >= (1 - DBL_P)) { 288 Dbl_rightshift_exponentmantissa(resultp1,resultp2, 289 1-dest_exponent); 290 } 291 else { 292 Dbl_setzero_exponentmantissa(resultp1,resultp2); 293 } 294 } 295 else Dbl_set_exponent(resultp1,dest_exponent); 296 Dbl_copytoptr(resultp1,resultp2,dstptr); 297 return(NOEXCEPTION); 298 } 299