1 /* $OpenBSD: dfdiv.c,v 1.4 2001/03/29 03:58:17 mickey Exp $ */ 2 3 /* 4 * Copyright 1996 1995 by Open Software Foundation, Inc. 5 * All Rights Reserved 6 * 7 * Permission to use, copy, modify, and distribute this software and 8 * its documentation for any purpose and without fee is hereby granted, 9 * provided that the above copyright notice appears in all copies and 10 * that both the copyright notice and this permission notice appear in 11 * supporting documentation. 12 * 13 * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE 14 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 15 * FOR A PARTICULAR PURPOSE. 16 * 17 * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR 18 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM 19 * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT, 20 * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION 21 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 22 * 23 */ 24 /* 25 * pmk1.1 26 */ 27 /* 28 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 29 * 30 * To anyone who acknowledges that this file is provided "AS IS" 31 * without any express or implied warranty: 32 * permission to use, copy, modify, and distribute this file 33 * for any purpose is hereby granted without fee, provided that 34 * the above copyright notice and this notice appears in all 35 * copies, and that the name of Hewlett-Packard Company not be 36 * used in advertising or publicity pertaining to distribution 37 * of the software without specific, written prior permission. 38 * Hewlett-Packard Company makes no representations about the 39 * suitability of this software for any purpose. 40 */ 41 42 #include "../spmath/float.h" 43 #include "../spmath/dbl_float.h" 44 45 /* 46 * Double Precision Floating-point Divide 47 */ 48 49 int 50 dbl_fdiv(srcptr1,srcptr2,dstptr,status) 51 52 dbl_floating_point *srcptr1, *srcptr2, *dstptr; 53 unsigned int *status; 54 { 55 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 56 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2; 57 register int dest_exponent, count; 58 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 59 int is_tiny; 60 61 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 62 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 63 /* 64 * set sign bit of result 65 */ 66 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 67 Dbl_setnegativezerop1(resultp1); 68 else Dbl_setzerop1(resultp1); 69 /* 70 * check first operand for NaN's or infinity 71 */ 72 if (Dbl_isinfinity_exponent(opnd1p1)) { 73 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 74 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 75 if (Dbl_isinfinity(opnd2p1,opnd2p2)) { 76 /* 77 * invalid since both operands 78 * are infinity 79 */ 80 if (Is_invalidtrap_enabled()) 81 return(INVALIDEXCEPTION); 82 Set_invalidflag(); 83 Dbl_makequietnan(resultp1,resultp2); 84 Dbl_copytoptr(resultp1,resultp2,dstptr); 85 return(NOEXCEPTION); 86 } 87 /* 88 * return infinity 89 */ 90 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 91 Dbl_copytoptr(resultp1,resultp2,dstptr); 92 return(NOEXCEPTION); 93 } 94 } 95 else { 96 /* 97 * is NaN; signaling or quiet? 98 */ 99 if (Dbl_isone_signaling(opnd1p1)) { 100 /* trap if INVALIDTRAP enabled */ 101 if (Is_invalidtrap_enabled()) 102 return(INVALIDEXCEPTION); 103 /* make NaN quiet */ 104 Set_invalidflag(); 105 Dbl_set_quiet(opnd1p1); 106 } 107 /* 108 * is second operand a signaling NaN? 109 */ 110 else if (Dbl_is_signalingnan(opnd2p1)) { 111 /* trap if INVALIDTRAP enabled */ 112 if (Is_invalidtrap_enabled()) 113 return(INVALIDEXCEPTION); 114 /* make NaN quiet */ 115 Set_invalidflag(); 116 Dbl_set_quiet(opnd2p1); 117 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 118 return(NOEXCEPTION); 119 } 120 /* 121 * return quiet NaN 122 */ 123 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 124 return(NOEXCEPTION); 125 } 126 } 127 /* 128 * check second operand for NaN's or infinity 129 */ 130 if (Dbl_isinfinity_exponent(opnd2p1)) { 131 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 132 /* 133 * return zero 134 */ 135 Dbl_setzero_exponentmantissa(resultp1,resultp2); 136 Dbl_copytoptr(resultp1,resultp2,dstptr); 137 return(NOEXCEPTION); 138 } 139 /* 140 * is NaN; signaling or quiet? 141 */ 142 if (Dbl_isone_signaling(opnd2p1)) { 143 /* trap if INVALIDTRAP enabled */ 144 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 145 /* make NaN quiet */ 146 Set_invalidflag(); 147 Dbl_set_quiet(opnd2p1); 148 } 149 /* 150 * return quiet NaN 151 */ 152 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 153 return(NOEXCEPTION); 154 } 155 /* 156 * check for division by zero 157 */ 158 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 159 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) { 160 /* invalid since both operands are zero */ 161 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 162 Set_invalidflag(); 163 Dbl_makequietnan(resultp1,resultp2); 164 Dbl_copytoptr(resultp1,resultp2,dstptr); 165 return(NOEXCEPTION); 166 } 167 if (Is_divisionbyzerotrap_enabled()) 168 return(DIVISIONBYZEROEXCEPTION); 169 Set_divisionbyzeroflag(); 170 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 171 Dbl_copytoptr(resultp1,resultp2,dstptr); 172 return(NOEXCEPTION); 173 } 174 /* 175 * Generate exponent 176 */ 177 dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS; 178 179 /* 180 * Generate mantissa 181 */ 182 if (Dbl_isnotzero_exponent(opnd1p1)) { 183 /* set hidden bit */ 184 Dbl_clear_signexponent_set_hidden(opnd1p1); 185 } 186 else { 187 /* check for zero */ 188 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 189 Dbl_setzero_exponentmantissa(resultp1,resultp2); 190 Dbl_copytoptr(resultp1,resultp2,dstptr); 191 return(NOEXCEPTION); 192 } 193 /* is denormalized, want to normalize */ 194 Dbl_clear_signexponent(opnd1p1); 195 Dbl_leftshiftby1(opnd1p1,opnd1p2); 196 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent); 197 } 198 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 199 if (Dbl_isnotzero_exponent(opnd2p1)) { 200 Dbl_clear_signexponent_set_hidden(opnd2p1); 201 } 202 else { 203 /* is denormalized; want to normalize */ 204 Dbl_clear_signexponent(opnd2p1); 205 Dbl_leftshiftby1(opnd2p1,opnd2p2); 206 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) { 207 dest_exponent+=8; 208 Dbl_leftshiftby8(opnd2p1,opnd2p2); 209 } 210 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) { 211 dest_exponent+=4; 212 Dbl_leftshiftby4(opnd2p1,opnd2p2); 213 } 214 while (Dbl_iszero_hidden(opnd2p1)) { 215 dest_exponent++; 216 Dbl_leftshiftby1(opnd2p1,opnd2p2); 217 } 218 } 219 220 /* Divide the source mantissas */ 221 222 /* 223 * A non-restoring divide algorithm is used. 224 */ 225 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 226 Dbl_setzero(opnd3p1,opnd3p2); 227 for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) { 228 Dbl_leftshiftby1(opnd1p1,opnd1p2); 229 Dbl_leftshiftby1(opnd3p1,opnd3p2); 230 if (Dbl_iszero_sign(opnd1p1)) { 231 Dbl_setone_lowmantissap2(opnd3p2); 232 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 233 } 234 else { 235 Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2); 236 } 237 } 238 if (count <= DBL_P) { 239 Dbl_leftshiftby1(opnd3p1,opnd3p2); 240 Dbl_setone_lowmantissap2(opnd3p2); 241 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count)); 242 if (Dbl_iszero_hidden(opnd3p1)) { 243 Dbl_leftshiftby1(opnd3p1,opnd3p2); 244 dest_exponent--; 245 } 246 } 247 else { 248 if (Dbl_iszero_hidden(opnd3p1)) { 249 /* need to get one more bit of result */ 250 Dbl_leftshiftby1(opnd1p1,opnd1p2); 251 Dbl_leftshiftby1(opnd3p1,opnd3p2); 252 if (Dbl_iszero_sign(opnd1p1)) { 253 Dbl_setone_lowmantissap2(opnd3p2); 254 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 255 } 256 else { 257 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 258 } 259 dest_exponent--; 260 } 261 if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE; 262 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2); 263 } 264 inexact = guardbit | stickybit; 265 266 /* 267 * round result 268 */ 269 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) { 270 Dbl_clear_signexponent(opnd3p1); 271 switch (Rounding_mode()) { 272 case ROUNDPLUS: 273 if (Dbl_iszero_sign(resultp1)) 274 Dbl_increment(opnd3p1,opnd3p2); 275 break; 276 case ROUNDMINUS: 277 if (Dbl_isone_sign(resultp1)) 278 Dbl_increment(opnd3p1,opnd3p2); 279 break; 280 case ROUNDNEAREST: 281 if (guardbit && (stickybit || 282 Dbl_isone_lowmantissap2(opnd3p2))) { 283 Dbl_increment(opnd3p1,opnd3p2); 284 } 285 } 286 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++; 287 } 288 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2); 289 290 /* 291 * Test for overflow 292 */ 293 if (dest_exponent >= DBL_INFINITY_EXPONENT) { 294 /* trap if OVERFLOWTRAP enabled */ 295 if (Is_overflowtrap_enabled()) { 296 /* 297 * Adjust bias of result 298 */ 299 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl); 300 Dbl_copytoptr(resultp1,resultp2,dstptr); 301 if (inexact) { 302 if (Is_inexacttrap_enabled()) 303 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 304 else 305 Set_inexactflag(); 306 } 307 return(OVERFLOWEXCEPTION); 308 } 309 Set_overflowflag(); 310 /* set result to infinity or largest number */ 311 Dbl_setoverflow(resultp1,resultp2); 312 inexact = TRUE; 313 } 314 /* 315 * Test for underflow 316 */ 317 else if (dest_exponent <= 0) { 318 /* trap if UNDERFLOWTRAP enabled */ 319 if (Is_underflowtrap_enabled()) { 320 /* 321 * Adjust bias of result 322 */ 323 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 324 Dbl_copytoptr(resultp1,resultp2,dstptr); 325 if (inexact) { 326 if (Is_inexacttrap_enabled()) 327 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 328 else 329 Set_inexactflag(); 330 } 331 return(UNDERFLOWEXCEPTION); 332 } 333 334 /* Determine if should set underflow flag */ 335 is_tiny = TRUE; 336 if (dest_exponent == 0 && inexact) { 337 switch (Rounding_mode()) { 338 case ROUNDPLUS: 339 if (Dbl_iszero_sign(resultp1)) { 340 Dbl_increment(opnd3p1,opnd3p2); 341 if (Dbl_isone_hiddenoverflow(opnd3p1)) 342 is_tiny = FALSE; 343 Dbl_decrement(opnd3p1,opnd3p2); 344 } 345 break; 346 case ROUNDMINUS: 347 if (Dbl_isone_sign(resultp1)) { 348 Dbl_increment(opnd3p1,opnd3p2); 349 if (Dbl_isone_hiddenoverflow(opnd3p1)) 350 is_tiny = FALSE; 351 Dbl_decrement(opnd3p1,opnd3p2); 352 } 353 break; 354 case ROUNDNEAREST: 355 if (guardbit && (stickybit || 356 Dbl_isone_lowmantissap2(opnd3p2))) { 357 Dbl_increment(opnd3p1,opnd3p2); 358 if (Dbl_isone_hiddenoverflow(opnd3p1)) 359 is_tiny = FALSE; 360 Dbl_decrement(opnd3p1,opnd3p2); 361 } 362 break; 363 } 364 } 365 366 /* 367 * denormalize result or set to signed zero 368 */ 369 stickybit = inexact; 370 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit, 371 stickybit,inexact); 372 373 /* return rounded number */ 374 if (inexact) { 375 switch (Rounding_mode()) { 376 case ROUNDPLUS: 377 if (Dbl_iszero_sign(resultp1)) { 378 Dbl_increment(opnd3p1,opnd3p2); 379 } 380 break; 381 case ROUNDMINUS: 382 if (Dbl_isone_sign(resultp1)) { 383 Dbl_increment(opnd3p1,opnd3p2); 384 } 385 break; 386 case ROUNDNEAREST: 387 if (guardbit && (stickybit || 388 Dbl_isone_lowmantissap2(opnd3p2))) { 389 Dbl_increment(opnd3p1,opnd3p2); 390 } 391 break; 392 } 393 if (is_tiny) 394 Set_underflowflag(); 395 } 396 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2); 397 } 398 else Dbl_set_exponent(resultp1,dest_exponent); 399 Dbl_copytoptr(resultp1,resultp2,dstptr); 400 401 /* check for inexact */ 402 if (inexact) { 403 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 404 else Set_inexactflag(); 405 } 406 return(NOEXCEPTION); 407 } 408