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