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