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