1 /* $NetBSD: sfrem.c,v 1.1 2002/06/05 01:04:26 fredette Exp $ */ 2 3 /* $OpenBSD: sfrem.c,v 1.4 2001/03/29 03:58:19 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/sgl_float.h" 47 48 /* 49 * Single Precision Floating-point Remainder 50 */ 51 int 52 sgl_frem(srcptr1,srcptr2,dstptr,status) 53 54 sgl_floating_point *srcptr1, *srcptr2, *dstptr; 55 unsigned int *status; 56 { 57 register unsigned int opnd1, opnd2, result; 58 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount; 59 register int roundup = FALSE; 60 61 opnd1 = *srcptr1; 62 opnd2 = *srcptr2; 63 /* 64 * check first operand for NaN's or infinity 65 */ 66 if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) { 67 if (Sgl_iszero_mantissa(opnd1)) { 68 if (Sgl_isnotnan(opnd2)) { 69 /* invalid since first operand is infinity */ 70 if (Is_invalidtrap_enabled()) 71 return(INVALIDEXCEPTION); 72 Set_invalidflag(); 73 Sgl_makequietnan(result); 74 *dstptr = result; 75 return(NOEXCEPTION); 76 } 77 } 78 else { 79 /* 80 * is NaN; signaling or quiet? 81 */ 82 if (Sgl_isone_signaling(opnd1)) { 83 /* trap if INVALIDTRAP enabled */ 84 if (Is_invalidtrap_enabled()) 85 return(INVALIDEXCEPTION); 86 /* make NaN quiet */ 87 Set_invalidflag(); 88 Sgl_set_quiet(opnd1); 89 } 90 /* 91 * is second operand a signaling NaN? 92 */ 93 else if (Sgl_is_signalingnan(opnd2)) { 94 /* trap if INVALIDTRAP enabled */ 95 if (Is_invalidtrap_enabled()) 96 return(INVALIDEXCEPTION); 97 /* make NaN quiet */ 98 Set_invalidflag(); 99 Sgl_set_quiet(opnd2); 100 *dstptr = opnd2; 101 return(NOEXCEPTION); 102 } 103 /* 104 * return quiet NaN 105 */ 106 *dstptr = opnd1; 107 return(NOEXCEPTION); 108 } 109 } 110 /* 111 * check second operand for NaN's or infinity 112 */ 113 if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) { 114 if (Sgl_iszero_mantissa(opnd2)) { 115 /* 116 * return first operand 117 */ 118 *dstptr = opnd1; 119 return(NOEXCEPTION); 120 } 121 /* 122 * is NaN; signaling or quiet? 123 */ 124 if (Sgl_isone_signaling(opnd2)) { 125 /* trap if INVALIDTRAP enabled */ 126 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 127 /* make NaN quiet */ 128 Set_invalidflag(); 129 Sgl_set_quiet(opnd2); 130 } 131 /* 132 * return quiet NaN 133 */ 134 *dstptr = opnd2; 135 return(NOEXCEPTION); 136 } 137 /* 138 * check second operand for zero 139 */ 140 if (Sgl_iszero_exponentmantissa(opnd2)) { 141 /* invalid since second operand is zero */ 142 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 143 Set_invalidflag(); 144 Sgl_makequietnan(result); 145 *dstptr = result; 146 return(NOEXCEPTION); 147 } 148 149 /* 150 * get sign of result 151 */ 152 result = opnd1; 153 154 /* 155 * check for denormalized operands 156 */ 157 if (opnd1_exponent == 0) { 158 /* check for zero */ 159 if (Sgl_iszero_mantissa(opnd1)) { 160 *dstptr = opnd1; 161 return(NOEXCEPTION); 162 } 163 /* normalize, then continue */ 164 opnd1_exponent = 1; 165 Sgl_normalize(opnd1,opnd1_exponent); 166 } 167 else { 168 Sgl_clear_signexponent_set_hidden(opnd1); 169 } 170 if (opnd2_exponent == 0) { 171 /* normalize, then continue */ 172 opnd2_exponent = 1; 173 Sgl_normalize(opnd2,opnd2_exponent); 174 } 175 else { 176 Sgl_clear_signexponent_set_hidden(opnd2); 177 } 178 179 /* find result exponent and divide step loop count */ 180 dest_exponent = opnd2_exponent - 1; 181 stepcount = opnd1_exponent - opnd2_exponent; 182 183 /* 184 * check for opnd1/opnd2 < 1 185 */ 186 if (stepcount < 0) { 187 /* 188 * check for opnd1/opnd2 > 1/2 189 * 190 * In this case n will round to 1, so 191 * r = opnd1 - opnd2 192 */ 193 if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) { 194 Sgl_all(result) = ~Sgl_all(result); /* set sign */ 195 /* align opnd2 with opnd1 */ 196 Sgl_leftshiftby1(opnd2); 197 Sgl_subtract(opnd2,opnd1,opnd2); 198 /* now normalize */ 199 while (Sgl_iszero_hidden(opnd2)) { 200 Sgl_leftshiftby1(opnd2); 201 dest_exponent--; 202 } 203 Sgl_set_exponentmantissa(result,opnd2); 204 goto testforunderflow; 205 } 206 /* 207 * opnd1/opnd2 <= 1/2 208 * 209 * In this case n will round to zero, so 210 * r = opnd1 211 */ 212 Sgl_set_exponentmantissa(result,opnd1); 213 dest_exponent = opnd1_exponent; 214 goto testforunderflow; 215 } 216 217 /* 218 * Generate result 219 * 220 * Do iterative subtract until remainder is less than operand 2. 221 */ 222 while (stepcount-- > 0 && Sgl_all(opnd1)) { 223 if (Sgl_isnotlessthan(opnd1,opnd2)) 224 Sgl_subtract(opnd1,opnd2,opnd1); 225 Sgl_leftshiftby1(opnd1); 226 } 227 /* 228 * Do last subtract, then determine which way to round if remainder 229 * is exactly 1/2 of opnd2 230 */ 231 if (Sgl_isnotlessthan(opnd1,opnd2)) { 232 Sgl_subtract(opnd1,opnd2,opnd1); 233 roundup = TRUE; 234 } 235 if (stepcount > 0 || Sgl_iszero(opnd1)) { 236 /* division is exact, remainder is zero */ 237 Sgl_setzero_exponentmantissa(result); 238 *dstptr = result; 239 return(NOEXCEPTION); 240 } 241 242 /* 243 * Check for cases where opnd1/opnd2 < n 244 * 245 * In this case the result's sign will be opposite that of 246 * opnd1. The mantissa also needs some correction. 247 */ 248 Sgl_leftshiftby1(opnd1); 249 if (Sgl_isgreaterthan(opnd1,opnd2)) { 250 Sgl_invert_sign(result); 251 Sgl_subtract((opnd2<<1),opnd1,opnd1); 252 } 253 /* check for remainder being exactly 1/2 of opnd2 */ 254 else if (Sgl_isequal(opnd1,opnd2) && roundup) { 255 Sgl_invert_sign(result); 256 } 257 258 /* normalize result's mantissa */ 259 while (Sgl_iszero_hidden(opnd1)) { 260 dest_exponent--; 261 Sgl_leftshiftby1(opnd1); 262 } 263 Sgl_set_exponentmantissa(result,opnd1); 264 265 /* 266 * Test for underflow 267 */ 268 testforunderflow: 269 if (dest_exponent <= 0) { 270 /* trap if UNDERFLOWTRAP enabled */ 271 if (Is_underflowtrap_enabled()) { 272 /* 273 * Adjust bias of result 274 */ 275 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 276 *dstptr = result; 277 /* frem is always exact */ 278 return(UNDERFLOWEXCEPTION); 279 } 280 /* 281 * denormalize result or set to signed zero 282 */ 283 if (dest_exponent >= (1 - SGL_P)) { 284 Sgl_rightshift_exponentmantissa(result,1-dest_exponent); 285 } else { 286 Sgl_setzero_exponentmantissa(result); 287 } 288 } 289 else Sgl_set_exponent(result,dest_exponent); 290 *dstptr = result; 291 return(NOEXCEPTION); 292 } 293