1 /* $NetBSD: sfmpy.c,v 1.1 2002/06/05 01:04:26 fredette Exp $ */ 2 3 /* $OpenBSD: sfmpy.c,v 1.4 2001/03/29 03:58:19 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 "../spmath/float.h" 45 #include "../spmath/sgl_float.h" 46 47 /* 48 * Single Precision Floating-point Multiply 49 */ 50 int 51 sgl_fmpy(srcptr1,srcptr2,dstptr,status) 52 53 sgl_floating_point *srcptr1, *srcptr2, *dstptr; 54 unsigned int *status; 55 { 56 register unsigned int opnd1, opnd2, opnd3, result; 57 register int dest_exponent, count; 58 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 59 int is_tiny; 60 61 opnd1 = *srcptr1; 62 opnd2 = *srcptr2; 63 /* 64 * set sign bit of result 65 */ 66 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result); 67 else Sgl_setzero(result); 68 /* 69 * check first operand for NaN's or infinity 70 */ 71 if (Sgl_isinfinity_exponent(opnd1)) { 72 if (Sgl_iszero_mantissa(opnd1)) { 73 if (Sgl_isnotnan(opnd2)) { 74 if (Sgl_iszero_exponentmantissa(opnd2)) { 75 /* 76 * invalid since operands are infinity 77 * and zero 78 */ 79 if (Is_invalidtrap_enabled()) 80 return(INVALIDEXCEPTION); 81 Set_invalidflag(); 82 Sgl_makequietnan(result); 83 *dstptr = result; 84 return(NOEXCEPTION); 85 } 86 /* 87 * return infinity 88 */ 89 Sgl_setinfinity_exponentmantissa(result); 90 *dstptr = result; 91 return(NOEXCEPTION); 92 } 93 } 94 else { 95 /* 96 * is NaN; signaling or quiet? 97 */ 98 if (Sgl_isone_signaling(opnd1)) { 99 /* trap if INVALIDTRAP enabled */ 100 if (Is_invalidtrap_enabled()) 101 return(INVALIDEXCEPTION); 102 /* make NaN quiet */ 103 Set_invalidflag(); 104 Sgl_set_quiet(opnd1); 105 } 106 /* 107 * is second operand a signaling NaN? 108 */ 109 else if (Sgl_is_signalingnan(opnd2)) { 110 /* trap if INVALIDTRAP enabled */ 111 if (Is_invalidtrap_enabled()) 112 return(INVALIDEXCEPTION); 113 /* make NaN quiet */ 114 Set_invalidflag(); 115 Sgl_set_quiet(opnd2); 116 *dstptr = opnd2; 117 return(NOEXCEPTION); 118 } 119 /* 120 * return quiet NaN 121 */ 122 *dstptr = opnd1; 123 return(NOEXCEPTION); 124 } 125 } 126 /* 127 * check second operand for NaN's or infinity 128 */ 129 if (Sgl_isinfinity_exponent(opnd2)) { 130 if (Sgl_iszero_mantissa(opnd2)) { 131 if (Sgl_iszero_exponentmantissa(opnd1)) { 132 /* invalid since operands are zero & infinity */ 133 if (Is_invalidtrap_enabled()) 134 return(INVALIDEXCEPTION); 135 Set_invalidflag(); 136 Sgl_makequietnan(opnd2); 137 *dstptr = opnd2; 138 return(NOEXCEPTION); 139 } 140 /* 141 * return infinity 142 */ 143 Sgl_setinfinity_exponentmantissa(result); 144 *dstptr = result; 145 return(NOEXCEPTION); 146 } 147 /* 148 * is NaN; signaling or quiet? 149 */ 150 if (Sgl_isone_signaling(opnd2)) { 151 /* trap if INVALIDTRAP enabled */ 152 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 153 154 /* make NaN quiet */ 155 Set_invalidflag(); 156 Sgl_set_quiet(opnd2); 157 } 158 /* 159 * return quiet NaN 160 */ 161 *dstptr = opnd2; 162 return(NOEXCEPTION); 163 } 164 /* 165 * Generate exponent 166 */ 167 dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS; 168 169 /* 170 * Generate mantissa 171 */ 172 if (Sgl_isnotzero_exponent(opnd1)) { 173 /* set hidden bit */ 174 Sgl_clear_signexponent_set_hidden(opnd1); 175 } 176 else { 177 /* check for zero */ 178 if (Sgl_iszero_mantissa(opnd1)) { 179 Sgl_setzero_exponentmantissa(result); 180 *dstptr = result; 181 return(NOEXCEPTION); 182 } 183 /* is denormalized, adjust exponent */ 184 Sgl_clear_signexponent(opnd1); 185 Sgl_leftshiftby1(opnd1); 186 Sgl_normalize(opnd1,dest_exponent); 187 } 188 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 189 if (Sgl_isnotzero_exponent(opnd2)) { 190 Sgl_clear_signexponent_set_hidden(opnd2); 191 } 192 else { 193 /* check for zero */ 194 if (Sgl_iszero_mantissa(opnd2)) { 195 Sgl_setzero_exponentmantissa(result); 196 *dstptr = result; 197 return(NOEXCEPTION); 198 } 199 /* is denormalized; want to normalize */ 200 Sgl_clear_signexponent(opnd2); 201 Sgl_leftshiftby1(opnd2); 202 Sgl_normalize(opnd2,dest_exponent); 203 } 204 205 /* Multiply two source mantissas together */ 206 207 Sgl_leftshiftby4(opnd2); /* make room for guard bits */ 208 Sgl_setzero(opnd3); 209 /* 210 * Four bits at a time are inspected in each loop, and a 211 * simple shift and add multiply algorithm is used. 212 */ 213 for (count=1;count<SGL_P;count+=4) { 214 stickybit |= Slow4(opnd3); 215 Sgl_rightshiftby4(opnd3); 216 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3); 217 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2); 218 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1); 219 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2); 220 Sgl_rightshiftby4(opnd1); 221 } 222 /* make sure result is left-justified */ 223 if (Sgl_iszero_sign(opnd3)) { 224 Sgl_leftshiftby1(opnd3); 225 } 226 else { 227 /* result mantissa >= 2. */ 228 dest_exponent++; 229 } 230 /* check for denormalized result */ 231 while (Sgl_iszero_sign(opnd3)) { 232 Sgl_leftshiftby1(opnd3); 233 dest_exponent--; 234 } 235 /* 236 * check for guard, sticky and inexact bits 237 */ 238 stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1); 239 guardbit = Sbit24(opnd3); 240 inexact = guardbit | stickybit; 241 242 /* re-align mantissa */ 243 Sgl_rightshiftby8(opnd3); 244 245 /* 246 * round result 247 */ 248 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 249 Sgl_clear_signexponent(opnd3); 250 switch (Rounding_mode()) { 251 case ROUNDPLUS: 252 if (Sgl_iszero_sign(result)) 253 Sgl_increment(opnd3); 254 break; 255 case ROUNDMINUS: 256 if (Sgl_isone_sign(result)) 257 Sgl_increment(opnd3); 258 break; 259 case ROUNDNEAREST: 260 if (guardbit && 261 (stickybit || Sgl_isone_lowmantissa(opnd3))) 262 Sgl_increment(opnd3); 263 break; 264 } 265 if (Sgl_isone_hidden(opnd3)) dest_exponent++; 266 } 267 Sgl_set_mantissa(result,opnd3); 268 269 /* 270 * Test for overflow 271 */ 272 if (dest_exponent >= SGL_INFINITY_EXPONENT) { 273 /* trap if OVERFLOWTRAP enabled */ 274 if (Is_overflowtrap_enabled()) { 275 /* 276 * Adjust bias of result 277 */ 278 Sgl_setwrapped_exponent(result,dest_exponent,ovfl); 279 *dstptr = result; 280 if (inexact) { 281 if (Is_inexacttrap_enabled()) 282 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 283 else Set_inexactflag(); 284 } 285 return(OVERFLOWEXCEPTION); 286 } 287 inexact = TRUE; 288 Set_overflowflag(); 289 /* set result to infinity or largest number */ 290 Sgl_setoverflow(result); 291 } 292 /* 293 * Test for underflow 294 */ 295 else if (dest_exponent <= 0) { 296 /* trap if UNDERFLOWTRAP enabled */ 297 if (Is_underflowtrap_enabled()) { 298 /* 299 * Adjust bias of result 300 */ 301 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 302 *dstptr = result; 303 if (inexact) { 304 if (Is_inexacttrap_enabled()) 305 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 306 else Set_inexactflag(); 307 } 308 return(UNDERFLOWEXCEPTION); 309 } 310 311 /* Determine if should set underflow flag */ 312 is_tiny = TRUE; 313 if (dest_exponent == 0 && inexact) { 314 switch (Rounding_mode()) { 315 case ROUNDPLUS: 316 if (Sgl_iszero_sign(result)) { 317 Sgl_increment(opnd3); 318 if (Sgl_isone_hiddenoverflow(opnd3)) 319 is_tiny = FALSE; 320 Sgl_decrement(opnd3); 321 } 322 break; 323 case ROUNDMINUS: 324 if (Sgl_isone_sign(result)) { 325 Sgl_increment(opnd3); 326 if (Sgl_isone_hiddenoverflow(opnd3)) 327 is_tiny = FALSE; 328 Sgl_decrement(opnd3); 329 } 330 break; 331 case ROUNDNEAREST: 332 if (guardbit && (stickybit || 333 Sgl_isone_lowmantissa(opnd3))) { 334 Sgl_increment(opnd3); 335 if (Sgl_isone_hiddenoverflow(opnd3)) 336 is_tiny = FALSE; 337 Sgl_decrement(opnd3); 338 } 339 break; 340 } 341 } 342 343 /* 344 * denormalize result or set to signed zero 345 */ 346 stickybit = inexact; 347 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact); 348 349 /* return zero or smallest number */ 350 if (inexact) { 351 switch (Rounding_mode()) { 352 case ROUNDPLUS: 353 if (Sgl_iszero_sign(result)) 354 Sgl_increment(opnd3); 355 break; 356 case ROUNDMINUS: 357 if (Sgl_isone_sign(result)) 358 Sgl_increment(opnd3); 359 break; 360 case ROUNDNEAREST: 361 if (guardbit && (stickybit || 362 Sgl_isone_lowmantissa(opnd3))) 363 Sgl_increment(opnd3); 364 break; 365 } 366 if (is_tiny) Set_underflowflag(); 367 } 368 Sgl_set_exponentmantissa(result,opnd3); 369 } 370 else Sgl_set_exponent(result,dest_exponent); 371 *dstptr = result; 372 373 /* check for inexact */ 374 if (inexact) { 375 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 376 else Set_inexactflag(); 377 } 378 return(NOEXCEPTION); 379 } 380