1 /* $OpenBSD: dfdiv.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 /* @(#)dfdiv.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:39 */ 16 17 #include "float.h" 18 #include "dbl_float.h" 19 20 /* 21 * Double Precision Floating-point Divide 22 */ 23 24 int 25 dbl_fdiv(srcptr1,srcptr2,dstptr,status) 26 dbl_floating_point *srcptr1, *srcptr2, *dstptr; 27 unsigned int *status; 28 { 29 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 30 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2; 31 register int dest_exponent, count; 32 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 33 int is_tiny; 34 35 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 36 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 37 /* 38 * set sign bit of result 39 */ 40 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 41 Dbl_setnegativezerop1(resultp1); 42 else Dbl_setzerop1(resultp1); 43 /* 44 * check first operand for NaN's or infinity 45 */ 46 if (Dbl_isinfinity_exponent(opnd1p1)) { 47 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 48 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 49 if (Dbl_isinfinity(opnd2p1,opnd2p2)) { 50 /* 51 * invalid since both operands 52 * are infinity 53 */ 54 if (Is_invalidtrap_enabled()) 55 return(INVALIDEXCEPTION); 56 Set_invalidflag(); 57 Dbl_makequietnan(resultp1,resultp2); 58 Dbl_copytoptr(resultp1,resultp2,dstptr); 59 return(NOEXCEPTION); 60 } 61 /* 62 * return infinity 63 */ 64 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 65 Dbl_copytoptr(resultp1,resultp2,dstptr); 66 return(NOEXCEPTION); 67 } 68 } 69 else { 70 /* 71 * is NaN; signaling or quiet? 72 */ 73 if (Dbl_isone_signaling(opnd1p1)) { 74 /* trap if INVALIDTRAP enabled */ 75 if (Is_invalidtrap_enabled()) 76 return(INVALIDEXCEPTION); 77 /* make NaN quiet */ 78 Set_invalidflag(); 79 Dbl_set_quiet(opnd1p1); 80 } 81 /* 82 * is second operand a signaling NaN? 83 */ 84 else if (Dbl_is_signalingnan(opnd2p1)) { 85 /* trap if INVALIDTRAP enabled */ 86 if (Is_invalidtrap_enabled()) 87 return(INVALIDEXCEPTION); 88 /* make NaN quiet */ 89 Set_invalidflag(); 90 Dbl_set_quiet(opnd2p1); 91 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 92 return(NOEXCEPTION); 93 } 94 /* 95 * return quiet NaN 96 */ 97 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 98 return(NOEXCEPTION); 99 } 100 } 101 /* 102 * check second operand for NaN's or infinity 103 */ 104 if (Dbl_isinfinity_exponent(opnd2p1)) { 105 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 106 /* 107 * return zero 108 */ 109 Dbl_setzero_exponentmantissa(resultp1,resultp2); 110 Dbl_copytoptr(resultp1,resultp2,dstptr); 111 return(NOEXCEPTION); 112 } 113 /* 114 * is NaN; signaling or quiet? 115 */ 116 if (Dbl_isone_signaling(opnd2p1)) { 117 /* trap if INVALIDTRAP enabled */ 118 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 119 /* make NaN quiet */ 120 Set_invalidflag(); 121 Dbl_set_quiet(opnd2p1); 122 } 123 /* 124 * return quiet NaN 125 */ 126 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 127 return(NOEXCEPTION); 128 } 129 /* 130 * check for division by zero 131 */ 132 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 133 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) { 134 /* invalid since both operands are zero */ 135 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 136 Set_invalidflag(); 137 Dbl_makequietnan(resultp1,resultp2); 138 Dbl_copytoptr(resultp1,resultp2,dstptr); 139 return(NOEXCEPTION); 140 } 141 if (Is_divisionbyzerotrap_enabled()) 142 return(DIVISIONBYZEROEXCEPTION); 143 Set_divisionbyzeroflag(); 144 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 145 Dbl_copytoptr(resultp1,resultp2,dstptr); 146 return(NOEXCEPTION); 147 } 148 /* 149 * Generate exponent 150 */ 151 dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS; 152 153 /* 154 * Generate mantissa 155 */ 156 if (Dbl_isnotzero_exponent(opnd1p1)) { 157 /* set hidden bit */ 158 Dbl_clear_signexponent_set_hidden(opnd1p1); 159 } 160 else { 161 /* check for zero */ 162 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 163 Dbl_setzero_exponentmantissa(resultp1,resultp2); 164 Dbl_copytoptr(resultp1,resultp2,dstptr); 165 return(NOEXCEPTION); 166 } 167 /* is denormalized, want to normalize */ 168 Dbl_clear_signexponent(opnd1p1); 169 Dbl_leftshiftby1(opnd1p1,opnd1p2); 170 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent); 171 } 172 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 173 if (Dbl_isnotzero_exponent(opnd2p1)) { 174 Dbl_clear_signexponent_set_hidden(opnd2p1); 175 } 176 else { 177 /* is denormalized; want to normalize */ 178 Dbl_clear_signexponent(opnd2p1); 179 Dbl_leftshiftby1(opnd2p1,opnd2p2); 180 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) { 181 dest_exponent+=8; 182 Dbl_leftshiftby8(opnd2p1,opnd2p2); 183 } 184 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) { 185 dest_exponent+=4; 186 Dbl_leftshiftby4(opnd2p1,opnd2p2); 187 } 188 while (Dbl_iszero_hidden(opnd2p1)) { 189 dest_exponent++; 190 Dbl_leftshiftby1(opnd2p1,opnd2p2); 191 } 192 } 193 194 /* Divide the source mantissas */ 195 196 /* 197 * A non-restoring divide algorithm is used. 198 */ 199 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 200 Dbl_setzero(opnd3p1,opnd3p2); 201 for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) { 202 Dbl_leftshiftby1(opnd1p1,opnd1p2); 203 Dbl_leftshiftby1(opnd3p1,opnd3p2); 204 if (Dbl_iszero_sign(opnd1p1)) { 205 Dbl_setone_lowmantissap2(opnd3p2); 206 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 207 } 208 else { 209 Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2); 210 } 211 } 212 if (count <= DBL_P) { 213 Dbl_leftshiftby1(opnd3p1,opnd3p2); 214 Dbl_setone_lowmantissap2(opnd3p2); 215 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count)); 216 if (Dbl_iszero_hidden(opnd3p1)) { 217 Dbl_leftshiftby1(opnd3p1,opnd3p2); 218 dest_exponent--; 219 } 220 } 221 else { 222 if (Dbl_iszero_hidden(opnd3p1)) { 223 /* need to get one more bit of result */ 224 Dbl_leftshiftby1(opnd1p1,opnd1p2); 225 Dbl_leftshiftby1(opnd3p1,opnd3p2); 226 if (Dbl_iszero_sign(opnd1p1)) { 227 Dbl_setone_lowmantissap2(opnd3p2); 228 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 229 } 230 else { 231 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 232 } 233 dest_exponent--; 234 } 235 if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE; 236 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2); 237 } 238 inexact = guardbit | stickybit; 239 240 /* 241 * round result 242 */ 243 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) { 244 Dbl_clear_signexponent(opnd3p1); 245 switch (Rounding_mode()) { 246 case ROUNDPLUS: 247 if (Dbl_iszero_sign(resultp1)) 248 Dbl_increment(opnd3p1,opnd3p2); 249 break; 250 case ROUNDMINUS: 251 if (Dbl_isone_sign(resultp1)) 252 Dbl_increment(opnd3p1,opnd3p2); 253 break; 254 case ROUNDNEAREST: 255 if (guardbit && (stickybit || 256 Dbl_isone_lowmantissap2(opnd3p2))) { 257 Dbl_increment(opnd3p1,opnd3p2); 258 } 259 } 260 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++; 261 } 262 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2); 263 264 /* 265 * Test for overflow 266 */ 267 if (dest_exponent >= DBL_INFINITY_EXPONENT) { 268 /* trap if OVERFLOWTRAP enabled */ 269 if (Is_overflowtrap_enabled()) { 270 /* 271 * Adjust bias of result 272 */ 273 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl); 274 Dbl_copytoptr(resultp1,resultp2,dstptr); 275 if (inexact) { 276 if (Is_inexacttrap_enabled()) 277 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 278 else 279 Set_inexactflag(); 280 } 281 return(OVERFLOWEXCEPTION); 282 } 283 Set_overflowflag(); 284 /* set result to infinity or largest number */ 285 Dbl_setoverflow(resultp1,resultp2); 286 inexact = TRUE; 287 } 288 /* 289 * Test for underflow 290 */ 291 else if (dest_exponent <= 0) { 292 /* trap if UNDERFLOWTRAP enabled */ 293 if (Is_underflowtrap_enabled()) { 294 /* 295 * Adjust bias of result 296 */ 297 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 298 Dbl_copytoptr(resultp1,resultp2,dstptr); 299 if (inexact) { 300 if (Is_inexacttrap_enabled()) 301 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 302 else 303 Set_inexactflag(); 304 } 305 return(UNDERFLOWEXCEPTION); 306 } 307 308 /* Determine if should set underflow flag */ 309 is_tiny = TRUE; 310 if (dest_exponent == 0 && inexact) { 311 switch (Rounding_mode()) { 312 case ROUNDPLUS: 313 if (Dbl_iszero_sign(resultp1)) { 314 Dbl_increment(opnd3p1,opnd3p2); 315 if (Dbl_isone_hiddenoverflow(opnd3p1)) 316 is_tiny = FALSE; 317 Dbl_decrement(opnd3p1,opnd3p2); 318 } 319 break; 320 case ROUNDMINUS: 321 if (Dbl_isone_sign(resultp1)) { 322 Dbl_increment(opnd3p1,opnd3p2); 323 if (Dbl_isone_hiddenoverflow(opnd3p1)) 324 is_tiny = FALSE; 325 Dbl_decrement(opnd3p1,opnd3p2); 326 } 327 break; 328 case ROUNDNEAREST: 329 if (guardbit && (stickybit || 330 Dbl_isone_lowmantissap2(opnd3p2))) { 331 Dbl_increment(opnd3p1,opnd3p2); 332 if (Dbl_isone_hiddenoverflow(opnd3p1)) 333 is_tiny = FALSE; 334 Dbl_decrement(opnd3p1,opnd3p2); 335 } 336 break; 337 } 338 } 339 340 /* 341 * denormalize result or set to signed zero 342 */ 343 stickybit = inexact; 344 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit, 345 stickybit,inexact); 346 347 /* return rounded number */ 348 if (inexact) { 349 switch (Rounding_mode()) { 350 case ROUNDPLUS: 351 if (Dbl_iszero_sign(resultp1)) { 352 Dbl_increment(opnd3p1,opnd3p2); 353 } 354 break; 355 case ROUNDMINUS: 356 if (Dbl_isone_sign(resultp1)) { 357 Dbl_increment(opnd3p1,opnd3p2); 358 } 359 break; 360 case ROUNDNEAREST: 361 if (guardbit && (stickybit || 362 Dbl_isone_lowmantissap2(opnd3p2))) { 363 Dbl_increment(opnd3p1,opnd3p2); 364 } 365 break; 366 } 367 if (is_tiny) 368 Set_underflowflag(); 369 } 370 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2); 371 } 372 else Dbl_set_exponent(resultp1,dest_exponent); 373 Dbl_copytoptr(resultp1,resultp2,dstptr); 374 375 /* check for inexact */ 376 if (inexact) { 377 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 378 else Set_inexactflag(); 379 } 380 return(NOEXCEPTION); 381 } 382