dfrem.c revision 1.4
1 /* $NetBSD: dfrem.c,v 1.4 2007/02/22 05:46:29 thorpej 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 #include <sys/cdefs.h> 45 __KERNEL_RCSID(0, "$NetBSD: dfrem.c,v 1.4 2007/02/22 05:46:29 thorpej Exp $"); 46 47 #include "../spmath/float.h" 48 #include "../spmath/dbl_float.h" 49 50 /* 51 * Double Precision Floating-point Remainder 52 */ 53 int 54 dbl_frem(srcptr1,srcptr2,dstptr,status) 55 56 dbl_floating_point *srcptr1, *srcptr2, *dstptr; 57 unsigned int *status; 58 { 59 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 60 register unsigned int resultp1, resultp2; 61 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount; 62 register int roundup = false; 63 64 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 65 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 66 /* 67 * check first operand for NaN's or infinity 68 */ 69 if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) { 70 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 71 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 72 /* invalid since first operand is infinity */ 73 if (Is_invalidtrap_enabled()) 74 return(INVALIDEXCEPTION); 75 Set_invalidflag(); 76 Dbl_makequietnan(resultp1,resultp2); 77 Dbl_copytoptr(resultp1,resultp2,dstptr); 78 return(NOEXCEPTION); 79 } 80 } 81 else { 82 /* 83 * is NaN; signaling or quiet? 84 */ 85 if (Dbl_isone_signaling(opnd1p1)) { 86 /* trap if INVALIDTRAP enabled */ 87 if (Is_invalidtrap_enabled()) 88 return(INVALIDEXCEPTION); 89 /* make NaN quiet */ 90 Set_invalidflag(); 91 Dbl_set_quiet(opnd1p1); 92 } 93 /* 94 * is second operand a signaling NaN? 95 */ 96 else if (Dbl_is_signalingnan(opnd2p1)) { 97 /* trap if INVALIDTRAP enabled */ 98 if (Is_invalidtrap_enabled()) 99 return(INVALIDEXCEPTION); 100 /* make NaN quiet */ 101 Set_invalidflag(); 102 Dbl_set_quiet(opnd2p1); 103 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 104 return(NOEXCEPTION); 105 } 106 /* 107 * return quiet NaN 108 */ 109 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 110 return(NOEXCEPTION); 111 } 112 } 113 /* 114 * check second operand for NaN's or infinity 115 */ 116 if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) { 117 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 118 /* 119 * return first operand 120 */ 121 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 122 return(NOEXCEPTION); 123 } 124 /* 125 * is NaN; signaling or quiet? 126 */ 127 if (Dbl_isone_signaling(opnd2p1)) { 128 /* trap if INVALIDTRAP enabled */ 129 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 130 /* make NaN quiet */ 131 Set_invalidflag(); 132 Dbl_set_quiet(opnd2p1); 133 } 134 /* 135 * return quiet NaN 136 */ 137 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 138 return(NOEXCEPTION); 139 } 140 /* 141 * check second operand for zero 142 */ 143 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 144 /* invalid since second operand is zero */ 145 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 146 Set_invalidflag(); 147 Dbl_makequietnan(resultp1,resultp2); 148 Dbl_copytoptr(resultp1,resultp2,dstptr); 149 return(NOEXCEPTION); 150 } 151 152 /* 153 * get sign of result 154 */ 155 resultp1 = opnd1p1; 156 157 /* 158 * check for denormalized operands 159 */ 160 if (opnd1_exponent == 0) { 161 /* check for zero */ 162 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 163 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 164 return(NOEXCEPTION); 165 } 166 /* normalize, then continue */ 167 opnd1_exponent = 1; 168 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent); 169 } 170 else { 171 Dbl_clear_signexponent_set_hidden(opnd1p1); 172 } 173 if (opnd2_exponent == 0) { 174 /* normalize, then continue */ 175 opnd2_exponent = 1; 176 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent); 177 } 178 else { 179 Dbl_clear_signexponent_set_hidden(opnd2p1); 180 } 181 182 /* find result exponent and divide step loop count */ 183 dest_exponent = opnd2_exponent - 1; 184 stepcount = opnd1_exponent - opnd2_exponent; 185 186 /* 187 * check for opnd1/opnd2 < 1 188 */ 189 if (stepcount < 0) { 190 /* 191 * check for opnd1/opnd2 > 1/2 192 * 193 * In this case n will round to 1, so 194 * r = opnd1 - opnd2 195 */ 196 if (stepcount == -1 && 197 Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 198 /* set sign */ 199 Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1); 200 /* align opnd2 with opnd1 */ 201 Dbl_leftshiftby1(opnd2p1,opnd2p2); 202 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2, 203 opnd2p1,opnd2p2); 204 /* now normalize */ 205 while (Dbl_iszero_hidden(opnd2p1)) { 206 Dbl_leftshiftby1(opnd2p1,opnd2p2); 207 dest_exponent--; 208 } 209 Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2); 210 goto testforunderflow; 211 } 212 /* 213 * opnd1/opnd2 <= 1/2 214 * 215 * In this case n will round to zero, so 216 * r = opnd1 217 */ 218 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 219 dest_exponent = opnd1_exponent; 220 goto testforunderflow; 221 } 222 223 /* 224 * Generate result 225 * 226 * Do iterative subtract until remainder is less than operand 2. 227 */ 228 while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) { 229 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 230 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 231 } 232 Dbl_leftshiftby1(opnd1p1,opnd1p2); 233 } 234 /* 235 * Do last subtract, then determine which way to round if remainder 236 * is exactly 1/2 of opnd2 237 */ 238 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 239 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 240 roundup = true; 241 } 242 if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) { 243 /* division is exact, remainder is zero */ 244 Dbl_setzero_exponentmantissa(resultp1,resultp2); 245 Dbl_copytoptr(resultp1,resultp2,dstptr); 246 return(NOEXCEPTION); 247 } 248 249 /* 250 * Check for cases where opnd1/opnd2 < n 251 * 252 * In this case the result's sign will be opposite that of 253 * opnd1. The mantissa also needs some correction. 254 */ 255 Dbl_leftshiftby1(opnd1p1,opnd1p2); 256 if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 257 Dbl_invert_sign(resultp1); 258 Dbl_leftshiftby1(opnd2p1,opnd2p2); 259 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2); 260 } 261 /* check for remainder being exactly 1/2 of opnd2 */ 262 else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) { 263 Dbl_invert_sign(resultp1); 264 } 265 266 /* normalize result's mantissa */ 267 while (Dbl_iszero_hidden(opnd1p1)) { 268 dest_exponent--; 269 Dbl_leftshiftby1(opnd1p1,opnd1p2); 270 } 271 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 272 273 /* 274 * Test for underflow 275 */ 276 testforunderflow: 277 if (dest_exponent <= 0) { 278 /* trap if UNDERFLOWTRAP enabled */ 279 if (Is_underflowtrap_enabled()) { 280 /* 281 * Adjust bias of result 282 */ 283 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 284 /* frem is always exact */ 285 Dbl_copytoptr(resultp1,resultp2,dstptr); 286 return(UNDERFLOWEXCEPTION); 287 } 288 /* 289 * denormalize result or set to signed zero 290 */ 291 if (dest_exponent >= (1 - DBL_P)) { 292 Dbl_rightshift_exponentmantissa(resultp1,resultp2, 293 1-dest_exponent); 294 } 295 else { 296 Dbl_setzero_exponentmantissa(resultp1,resultp2); 297 } 298 } 299 else Dbl_set_exponent(resultp1,dest_exponent); 300 Dbl_copytoptr(resultp1,resultp2,dstptr); 301 return(NOEXCEPTION); 302 } 303
Indexes created Fri Sep 26 20:09:58 GMT 2025