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