1 /* $OpenBSD: dfmpy.c,v 1.5 2002/05/07 22:19:30 mickey Exp $ */ 2 /* 3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY 4 To anyone who acknowledges that this file is provided "AS IS" 5 without any express or implied warranty: 6 permission to use, copy, modify, and distribute this file 7 for any purpose is hereby granted without fee, provided that 8 the above copyright notice and this notice appears in all 9 copies, and that the name of Hewlett-Packard Company not be 10 used in advertising or publicity pertaining to distribution 11 of the software without specific, written prior permission. 12 Hewlett-Packard Company makes no representations about the 13 suitability of this software for any purpose. 14 */ 15 /* @(#)dfmpy.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:41 */ 16 17 #include "float.h" 18 #include "dbl_float.h" 19 20 /* 21 * Double Precision Floating-point Multiply 22 */ 23 24 int 25 dbl_fmpy(srcptr1,srcptr2,dstptr,status) 26 dbl_floating_point *srcptr1, *srcptr2, *dstptr; 27 unsigned int *status; 28 { 29 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 30 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2; 31 register int dest_exponent, count; 32 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 33 int is_tiny; 34 35 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 36 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 37 38 /* 39 * set sign bit of result 40 */ 41 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 42 Dbl_setnegativezerop1(resultp1); 43 else Dbl_setzerop1(resultp1); 44 /* 45 * check first operand for NaN's or infinity 46 */ 47 if (Dbl_isinfinity_exponent(opnd1p1)) { 48 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 49 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 50 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 51 /* 52 * invalid since operands are infinity 53 * and zero 54 */ 55 if (Is_invalidtrap_enabled()) 56 return(INVALIDEXCEPTION); 57 Set_invalidflag(); 58 Dbl_makequietnan(resultp1,resultp2); 59 Dbl_copytoptr(resultp1,resultp2,dstptr); 60 return(NOEXCEPTION); 61 } 62 /* 63 * return infinity 64 */ 65 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 66 Dbl_copytoptr(resultp1,resultp2,dstptr); 67 return(NOEXCEPTION); 68 } 69 } 70 else { 71 /* 72 * is NaN; signaling or quiet? 73 */ 74 if (Dbl_isone_signaling(opnd1p1)) { 75 /* trap if INVALIDTRAP enabled */ 76 if (Is_invalidtrap_enabled()) 77 return(INVALIDEXCEPTION); 78 /* make NaN quiet */ 79 Set_invalidflag(); 80 Dbl_set_quiet(opnd1p1); 81 } 82 /* 83 * is second operand a signaling NaN? 84 */ 85 else if (Dbl_is_signalingnan(opnd2p1)) { 86 /* trap if INVALIDTRAP enabled */ 87 if (Is_invalidtrap_enabled()) 88 return(INVALIDEXCEPTION); 89 /* make NaN quiet */ 90 Set_invalidflag(); 91 Dbl_set_quiet(opnd2p1); 92 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 93 return(NOEXCEPTION); 94 } 95 /* 96 * return quiet NaN 97 */ 98 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 99 return(NOEXCEPTION); 100 } 101 } 102 /* 103 * check second operand for NaN's or infinity 104 */ 105 if (Dbl_isinfinity_exponent(opnd2p1)) { 106 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 107 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) { 108 /* invalid since operands are zero & infinity */ 109 if (Is_invalidtrap_enabled()) 110 return(INVALIDEXCEPTION); 111 Set_invalidflag(); 112 Dbl_makequietnan(opnd2p1,opnd2p2); 113 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 114 return(NOEXCEPTION); 115 } 116 /* 117 * return infinity 118 */ 119 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 120 Dbl_copytoptr(resultp1,resultp2,dstptr); 121 return(NOEXCEPTION); 122 } 123 /* 124 * is NaN; signaling or quiet? 125 */ 126 if (Dbl_isone_signaling(opnd2p1)) { 127 /* trap if INVALIDTRAP enabled */ 128 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 129 /* make NaN quiet */ 130 Set_invalidflag(); 131 Dbl_set_quiet(opnd2p1); 132 } 133 /* 134 * return quiet NaN 135 */ 136 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 137 return(NOEXCEPTION); 138 } 139 /* 140 * Generate exponent 141 */ 142 dest_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) -DBL_BIAS; 143 144 /* 145 * Generate mantissa 146 */ 147 if (Dbl_isnotzero_exponent(opnd1p1)) { 148 /* set hidden bit */ 149 Dbl_clear_signexponent_set_hidden(opnd1p1); 150 } 151 else { 152 /* check for zero */ 153 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 154 Dbl_setzero_exponentmantissa(resultp1,resultp2); 155 Dbl_copytoptr(resultp1,resultp2,dstptr); 156 return(NOEXCEPTION); 157 } 158 /* is denormalized, adjust exponent */ 159 Dbl_clear_signexponent(opnd1p1); 160 Dbl_leftshiftby1(opnd1p1,opnd1p2); 161 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent); 162 } 163 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 164 if (Dbl_isnotzero_exponent(opnd2p1)) { 165 Dbl_clear_signexponent_set_hidden(opnd2p1); 166 } 167 else { 168 /* check for zero */ 169 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 170 Dbl_setzero_exponentmantissa(resultp1,resultp2); 171 Dbl_copytoptr(resultp1,resultp2,dstptr); 172 return(NOEXCEPTION); 173 } 174 /* is denormalized; want to normalize */ 175 Dbl_clear_signexponent(opnd2p1); 176 Dbl_leftshiftby1(opnd2p1,opnd2p2); 177 Dbl_normalize(opnd2p1,opnd2p2,dest_exponent); 178 } 179 180 /* Multiply two source mantissas together */ 181 182 /* make room for guard bits */ 183 Dbl_leftshiftby7(opnd2p1,opnd2p2); 184 Dbl_setzero(opnd3p1,opnd3p2); 185 /* 186 * Four bits at a time are inspected in each loop, and a 187 * simple shift and add multiply algorithm is used. 188 */ 189 for (count=1;count<=DBL_P;count+=4) { 190 stickybit |= Dlow4p2(opnd3p2); 191 Dbl_rightshiftby4(opnd3p1,opnd3p2); 192 if (Dbit28p2(opnd1p2)) { 193 /* Twoword_add should be an ADDC followed by an ADD. */ 194 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<3 | opnd2p2>>29, 195 opnd2p2<<3); 196 } 197 if (Dbit29p2(opnd1p2)) { 198 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<2 | opnd2p2>>30, 199 opnd2p2<<2); 200 } 201 if (Dbit30p2(opnd1p2)) { 202 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<1 | opnd2p2>>31, 203 opnd2p2<<1); 204 } 205 if (Dbit31p2(opnd1p2)) { 206 Twoword_add(opnd3p1, opnd3p2, opnd2p1, opnd2p2); 207 } 208 Dbl_rightshiftby4(opnd1p1,opnd1p2); 209 } 210 if (Dbit3p1(opnd3p1)==0) { 211 Dbl_leftshiftby1(opnd3p1,opnd3p2); 212 } 213 else { 214 /* result mantissa >= 2. */ 215 dest_exponent++; 216 } 217 /* check for denormalized result */ 218 while (Dbit3p1(opnd3p1)==0) { 219 Dbl_leftshiftby1(opnd3p1,opnd3p2); 220 dest_exponent--; 221 } 222 /* 223 * check for guard, sticky and inexact bits 224 */ 225 stickybit |= Dallp2(opnd3p2) << 25; 226 guardbit = (Dallp2(opnd3p2) << 24) >> 31; 227 inexact = guardbit | stickybit; 228 229 /* align result mantissa */ 230 Dbl_rightshiftby8(opnd3p1,opnd3p2); 231 232 /* 233 * round result 234 */ 235 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 236 Dbl_clear_signexponent(opnd3p1); 237 switch (Rounding_mode()) { 238 case ROUNDPLUS: 239 if (Dbl_iszero_sign(resultp1)) 240 Dbl_increment(opnd3p1,opnd3p2); 241 break; 242 case ROUNDMINUS: 243 if (Dbl_isone_sign(resultp1)) 244 Dbl_increment(opnd3p1,opnd3p2); 245 break; 246 case ROUNDNEAREST: 247 if (guardbit && 248 (stickybit || Dbl_isone_lowmantissap2(opnd3p2))) 249 Dbl_increment(opnd3p1,opnd3p2); 250 break; 251 } 252 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++; 253 } 254 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2); 255 256 /* 257 * Test for overflow 258 */ 259 if (dest_exponent >= DBL_INFINITY_EXPONENT) { 260 /* trap if OVERFLOWTRAP enabled */ 261 if (Is_overflowtrap_enabled()) { 262 /* 263 * Adjust bias of result 264 */ 265 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl); 266 Dbl_copytoptr(resultp1,resultp2,dstptr); 267 if (inexact) { 268 if (Is_inexacttrap_enabled()) 269 return (OVERFLOWEXCEPTION | INEXACTEXCEPTION); 270 else 271 Set_inexactflag(); 272 } 273 return (OVERFLOWEXCEPTION); 274 } 275 inexact = TRUE; 276 Set_overflowflag(); 277 /* set result to infinity or largest number */ 278 Dbl_setoverflow(resultp1,resultp2); 279 } 280 /* 281 * Test for underflow 282 */ 283 else if (dest_exponent <= 0) { 284 /* trap if UNDERFLOWTRAP enabled */ 285 if (Is_underflowtrap_enabled()) { 286 /* 287 * Adjust bias of result 288 */ 289 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 290 Dbl_copytoptr(resultp1,resultp2,dstptr); 291 if (inexact) { 292 if (Is_inexacttrap_enabled()) 293 return (UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 294 else 295 Set_inexactflag(); 296 } 297 return (UNDERFLOWEXCEPTION); 298 } 299 300 /* Determine if should set underflow flag */ 301 is_tiny = TRUE; 302 if (dest_exponent == 0 && inexact) { 303 switch (Rounding_mode()) { 304 case ROUNDPLUS: 305 if (Dbl_iszero_sign(resultp1)) { 306 Dbl_increment(opnd3p1,opnd3p2); 307 if (Dbl_isone_hiddenoverflow(opnd3p1)) 308 is_tiny = FALSE; 309 Dbl_decrement(opnd3p1,opnd3p2); 310 } 311 break; 312 case ROUNDMINUS: 313 if (Dbl_isone_sign(resultp1)) { 314 Dbl_increment(opnd3p1,opnd3p2); 315 if (Dbl_isone_hiddenoverflow(opnd3p1)) 316 is_tiny = FALSE; 317 Dbl_decrement(opnd3p1,opnd3p2); 318 } 319 break; 320 case ROUNDNEAREST: 321 if (guardbit && (stickybit || 322 Dbl_isone_lowmantissap2(opnd3p2))) { 323 Dbl_increment(opnd3p1,opnd3p2); 324 if (Dbl_isone_hiddenoverflow(opnd3p1)) 325 is_tiny = FALSE; 326 Dbl_decrement(opnd3p1,opnd3p2); 327 } 328 break; 329 } 330 } 331 332 /* 333 * denormalize result or set to signed zero 334 */ 335 stickybit = inexact; 336 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit, 337 stickybit,inexact); 338 339 /* return zero or smallest number */ 340 if (inexact) { 341 switch (Rounding_mode()) { 342 case ROUNDPLUS: 343 if (Dbl_iszero_sign(resultp1)) { 344 Dbl_increment(opnd3p1,opnd3p2); 345 } 346 break; 347 case ROUNDMINUS: 348 if (Dbl_isone_sign(resultp1)) { 349 Dbl_increment(opnd3p1,opnd3p2); 350 } 351 break; 352 case ROUNDNEAREST: 353 if (guardbit && (stickybit || 354 Dbl_isone_lowmantissap2(opnd3p2))) { 355 Dbl_increment(opnd3p1,opnd3p2); 356 } 357 break; 358 } 359 if (is_tiny) Set_underflowflag(); 360 } 361 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2); 362 } 363 else Dbl_set_exponent(resultp1,dest_exponent); 364 /* check for inexact */ 365 Dbl_copytoptr(resultp1,resultp2,dstptr); 366 if (inexact) { 367 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 368 else Set_inexactflag(); 369 } 370 return(NOEXCEPTION); 371 } 372