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