1 /* $OpenBSD: sfdiv.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 /* @(#)sfdiv.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:05 */ 16 17 #include "float.h" 18 #include "sgl_float.h" 19 20 /* 21 * Single Precision Floating-point Divide 22 */ 23 int 24 sgl_fdiv(srcptr1,srcptr2,dstptr,status) 25 sgl_floating_point *srcptr1, *srcptr2, *dstptr; 26 unsigned int *status; 27 { 28 register unsigned int opnd1, opnd2, opnd3, result; 29 register int dest_exponent, count; 30 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 31 int is_tiny; 32 33 opnd1 = *srcptr1; 34 opnd2 = *srcptr2; 35 /* 36 * set sign bit of result 37 */ 38 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result); 39 else Sgl_setzero(result); 40 /* 41 * check first operand for NaN's or infinity 42 */ 43 if (Sgl_isinfinity_exponent(opnd1)) { 44 if (Sgl_iszero_mantissa(opnd1)) { 45 if (Sgl_isnotnan(opnd2)) { 46 if (Sgl_isinfinity(opnd2)) { 47 /* 48 * invalid since both operands 49 * are infinity 50 */ 51 if (Is_invalidtrap_enabled()) 52 return(INVALIDEXCEPTION); 53 Set_invalidflag(); 54 Sgl_makequietnan(result); 55 *dstptr = result; 56 return(NOEXCEPTION); 57 } 58 /* 59 * return infinity 60 */ 61 Sgl_setinfinity_exponentmantissa(result); 62 *dstptr = result; 63 return(NOEXCEPTION); 64 } 65 } 66 else { 67 /* 68 * is NaN; signaling or quiet? 69 */ 70 if (Sgl_isone_signaling(opnd1)) { 71 /* trap if INVALIDTRAP enabled */ 72 if (Is_invalidtrap_enabled()) 73 return(INVALIDEXCEPTION); 74 /* make NaN quiet */ 75 Set_invalidflag(); 76 Sgl_set_quiet(opnd1); 77 } 78 /* 79 * is second operand a signaling NaN? 80 */ 81 else if (Sgl_is_signalingnan(opnd2)) { 82 /* trap if INVALIDTRAP enabled */ 83 if (Is_invalidtrap_enabled()) 84 return(INVALIDEXCEPTION); 85 /* make NaN quiet */ 86 Set_invalidflag(); 87 Sgl_set_quiet(opnd2); 88 *dstptr = opnd2; 89 return(NOEXCEPTION); 90 } 91 /* 92 * return quiet NaN 93 */ 94 *dstptr = opnd1; 95 return(NOEXCEPTION); 96 } 97 } 98 /* 99 * check second operand for NaN's or infinity 100 */ 101 if (Sgl_isinfinity_exponent(opnd2)) { 102 if (Sgl_iszero_mantissa(opnd2)) { 103 /* 104 * return zero 105 */ 106 Sgl_setzero_exponentmantissa(result); 107 *dstptr = result; 108 return(NOEXCEPTION); 109 } 110 /* 111 * is NaN; signaling or quiet? 112 */ 113 if (Sgl_isone_signaling(opnd2)) { 114 /* trap if INVALIDTRAP enabled */ 115 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 116 /* make NaN quiet */ 117 Set_invalidflag(); 118 Sgl_set_quiet(opnd2); 119 } 120 /* 121 * return quiet NaN 122 */ 123 *dstptr = opnd2; 124 return(NOEXCEPTION); 125 } 126 /* 127 * check for division by zero 128 */ 129 if (Sgl_iszero_exponentmantissa(opnd2)) { 130 if (Sgl_iszero_exponentmantissa(opnd1)) { 131 /* invalid since both operands are zero */ 132 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 133 Set_invalidflag(); 134 Sgl_makequietnan(result); 135 *dstptr = result; 136 return(NOEXCEPTION); 137 } 138 if (Is_divisionbyzerotrap_enabled()) 139 return(DIVISIONBYZEROEXCEPTION); 140 Set_divisionbyzeroflag(); 141 Sgl_setinfinity_exponentmantissa(result); 142 *dstptr = result; 143 return(NOEXCEPTION); 144 } 145 /* 146 * Generate exponent 147 */ 148 dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS; 149 150 /* 151 * Generate mantissa 152 */ 153 if (Sgl_isnotzero_exponent(opnd1)) { 154 /* set hidden bit */ 155 Sgl_clear_signexponent_set_hidden(opnd1); 156 } 157 else { 158 /* check for zero */ 159 if (Sgl_iszero_mantissa(opnd1)) { 160 Sgl_setzero_exponentmantissa(result); 161 *dstptr = result; 162 return(NOEXCEPTION); 163 } 164 /* is denormalized; want to normalize */ 165 Sgl_clear_signexponent(opnd1); 166 Sgl_leftshiftby1(opnd1); 167 Sgl_normalize(opnd1,dest_exponent); 168 } 169 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 170 if (Sgl_isnotzero_exponent(opnd2)) { 171 Sgl_clear_signexponent_set_hidden(opnd2); 172 } 173 else { 174 /* is denormalized; want to normalize */ 175 Sgl_clear_signexponent(opnd2); 176 Sgl_leftshiftby1(opnd2); 177 while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) { 178 Sgl_leftshiftby8(opnd2); 179 dest_exponent += 8; 180 } 181 if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) { 182 Sgl_leftshiftby4(opnd2); 183 dest_exponent += 4; 184 } 185 while(Sgl_iszero_hidden(opnd2)) { 186 Sgl_leftshiftby1(opnd2); 187 dest_exponent += 1; 188 } 189 } 190 191 /* Divide the source mantissas */ 192 193 /* 194 * A non_restoring divide algorithm is used. 195 */ 196 Sgl_subtract(opnd1,opnd2,opnd1); 197 Sgl_setzero(opnd3); 198 for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) { 199 Sgl_leftshiftby1(opnd1); 200 Sgl_leftshiftby1(opnd3); 201 if (Sgl_iszero_sign(opnd1)) { 202 Sgl_setone_lowmantissa(opnd3); 203 Sgl_subtract(opnd1,opnd2,opnd1); 204 } 205 else Sgl_addition(opnd1,opnd2,opnd1); 206 } 207 if (count <= SGL_P) { 208 Sgl_leftshiftby1(opnd3); 209 Sgl_setone_lowmantissa(opnd3); 210 Sgl_leftshift(opnd3,SGL_P-count); 211 if (Sgl_iszero_hidden(opnd3)) { 212 Sgl_leftshiftby1(opnd3); 213 dest_exponent--; 214 } 215 } 216 else { 217 if (Sgl_iszero_hidden(opnd3)) { 218 /* need to get one more bit of result */ 219 Sgl_leftshiftby1(opnd1); 220 Sgl_leftshiftby1(opnd3); 221 if (Sgl_iszero_sign(opnd1)) { 222 Sgl_setone_lowmantissa(opnd3); 223 Sgl_subtract(opnd1,opnd2,opnd1); 224 } 225 else Sgl_addition(opnd1,opnd2,opnd1); 226 dest_exponent--; 227 } 228 if (Sgl_iszero_sign(opnd1)) guardbit = TRUE; 229 stickybit = Sgl_all(opnd1); 230 } 231 inexact = guardbit | stickybit; 232 233 /* 234 * round result 235 */ 236 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) { 237 Sgl_clear_signexponent(opnd3); 238 switch (Rounding_mode()) { 239 case ROUNDPLUS: 240 if (Sgl_iszero_sign(result)) 241 Sgl_increment_mantissa(opnd3); 242 break; 243 case ROUNDMINUS: 244 if (Sgl_isone_sign(result)) 245 Sgl_increment_mantissa(opnd3); 246 break; 247 case ROUNDNEAREST: 248 if (guardbit && 249 (stickybit || Sgl_isone_lowmantissa(opnd3))) 250 Sgl_increment_mantissa(opnd3); 251 } 252 if (Sgl_isone_hidden(opnd3)) dest_exponent++; 253 } 254 Sgl_set_mantissa(result,opnd3); 255 256 /* 257 * Test for overflow 258 */ 259 if (dest_exponent >= SGL_INFINITY_EXPONENT) { 260 /* trap if OVERFLOWTRAP enabled */ 261 if (Is_overflowtrap_enabled()) { 262 /* 263 * Adjust bias of result 264 */ 265 Sgl_setwrapped_exponent(result,dest_exponent,ovfl); 266 *dstptr = result; 267 if (inexact) { 268 if (Is_inexacttrap_enabled()) 269 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 270 else Set_inexactflag(); 271 } 272 return(OVERFLOWEXCEPTION); 273 } 274 Set_overflowflag(); 275 /* set result to infinity or largest number */ 276 Sgl_setoverflow(result); 277 inexact = TRUE; 278 } 279 /* 280 * Test for underflow 281 */ 282 else if (dest_exponent <= 0) { 283 /* trap if UNDERFLOWTRAP enabled */ 284 if (Is_underflowtrap_enabled()) { 285 /* 286 * Adjust bias of result 287 */ 288 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 289 *dstptr = result; 290 if (inexact) { 291 if (Is_inexacttrap_enabled()) 292 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 293 else Set_inexactflag(); 294 } 295 return(UNDERFLOWEXCEPTION); 296 } 297 298 /* Determine if should set underflow flag */ 299 is_tiny = TRUE; 300 if (dest_exponent == 0 && inexact) { 301 switch (Rounding_mode()) { 302 case ROUNDPLUS: 303 if (Sgl_iszero_sign(result)) { 304 Sgl_increment(opnd3); 305 if (Sgl_isone_hiddenoverflow(opnd3)) 306 is_tiny = FALSE; 307 Sgl_decrement(opnd3); 308 } 309 break; 310 case ROUNDMINUS: 311 if (Sgl_isone_sign(result)) { 312 Sgl_increment(opnd3); 313 if (Sgl_isone_hiddenoverflow(opnd3)) 314 is_tiny = FALSE; 315 Sgl_decrement(opnd3); 316 } 317 break; 318 case ROUNDNEAREST: 319 if (guardbit && (stickybit || 320 Sgl_isone_lowmantissa(opnd3))) { 321 Sgl_increment(opnd3); 322 if (Sgl_isone_hiddenoverflow(opnd3)) 323 is_tiny = FALSE; 324 Sgl_decrement(opnd3); 325 } 326 break; 327 } 328 } 329 330 /* 331 * denormalize result or set to signed zero 332 */ 333 stickybit = inexact; 334 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact); 335 336 /* return rounded number */ 337 if (inexact) { 338 switch (Rounding_mode()) { 339 case ROUNDPLUS: 340 if (Sgl_iszero_sign(result)) { 341 Sgl_increment(opnd3); 342 } 343 break; 344 case ROUNDMINUS: 345 if (Sgl_isone_sign(result)) { 346 Sgl_increment(opnd3); 347 } 348 break; 349 case ROUNDNEAREST: 350 if (guardbit && (stickybit || 351 Sgl_isone_lowmantissa(opnd3))) { 352 Sgl_increment(opnd3); 353 } 354 break; 355 } 356 if (is_tiny) Set_underflowflag(); 357 } 358 Sgl_set_exponentmantissa(result,opnd3); 359 } 360 else Sgl_set_exponent(result,dest_exponent); 361 *dstptr = result; 362 /* check for inexact */ 363 if (inexact) { 364 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 365 else Set_inexactflag(); 366 } 367 return(NOEXCEPTION); 368 } 369