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