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