1 /* $OpenBSD: sfmpy.c,v 1.4 2001/03/29 03:58:19 mickey Exp $ */ 2 3 /* 4 * Copyright 1996 1995 by Open Software Foundation, Inc. 5 * All Rights Reserved 6 * 7 * Permission to use, copy, modify, and distribute this software and 8 * its documentation for any purpose and without fee is hereby granted, 9 * provided that the above copyright notice appears in all copies and 10 * that both the copyright notice and this permission notice appear in 11 * supporting documentation. 12 * 13 * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE 14 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 15 * FOR A PARTICULAR PURPOSE. 16 * 17 * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR 18 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM 19 * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT, 20 * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION 21 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 22 * 23 */ 24 /* 25 * pmk1.1 26 */ 27 /* 28 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 29 * 30 * To anyone who acknowledges that this file is provided "AS IS" 31 * without any express or implied warranty: 32 * permission to use, copy, modify, and distribute this file 33 * for any purpose is hereby granted without fee, provided that 34 * the above copyright notice and this notice appears in all 35 * copies, and that the name of Hewlett-Packard Company not be 36 * used in advertising or publicity pertaining to distribution 37 * of the software without specific, written prior permission. 38 * Hewlett-Packard Company makes no representations about the 39 * suitability of this software for any purpose. 40 */ 41 42 #include "../spmath/float.h" 43 #include "../spmath/sgl_float.h" 44 45 /* 46 * Single Precision Floating-point Multiply 47 */ 48 int 49 sgl_fmpy(srcptr1,srcptr2,dstptr,status) 50 51 sgl_floating_point *srcptr1, *srcptr2, *dstptr; 52 unsigned int *status; 53 { 54 register unsigned int opnd1, opnd2, opnd3, result; 55 register int dest_exponent, count; 56 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 57 int is_tiny; 58 59 opnd1 = *srcptr1; 60 opnd2 = *srcptr2; 61 /* 62 * set sign bit of result 63 */ 64 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result); 65 else Sgl_setzero(result); 66 /* 67 * check first operand for NaN's or infinity 68 */ 69 if (Sgl_isinfinity_exponent(opnd1)) { 70 if (Sgl_iszero_mantissa(opnd1)) { 71 if (Sgl_isnotnan(opnd2)) { 72 if (Sgl_iszero_exponentmantissa(opnd2)) { 73 /* 74 * invalid since operands are infinity 75 * and zero 76 */ 77 if (Is_invalidtrap_enabled()) 78 return(INVALIDEXCEPTION); 79 Set_invalidflag(); 80 Sgl_makequietnan(result); 81 *dstptr = result; 82 return(NOEXCEPTION); 83 } 84 /* 85 * return infinity 86 */ 87 Sgl_setinfinity_exponentmantissa(result); 88 *dstptr = result; 89 return(NOEXCEPTION); 90 } 91 } 92 else { 93 /* 94 * is NaN; signaling or quiet? 95 */ 96 if (Sgl_isone_signaling(opnd1)) { 97 /* trap if INVALIDTRAP enabled */ 98 if (Is_invalidtrap_enabled()) 99 return(INVALIDEXCEPTION); 100 /* make NaN quiet */ 101 Set_invalidflag(); 102 Sgl_set_quiet(opnd1); 103 } 104 /* 105 * is second operand a signaling NaN? 106 */ 107 else if (Sgl_is_signalingnan(opnd2)) { 108 /* trap if INVALIDTRAP enabled */ 109 if (Is_invalidtrap_enabled()) 110 return(INVALIDEXCEPTION); 111 /* make NaN quiet */ 112 Set_invalidflag(); 113 Sgl_set_quiet(opnd2); 114 *dstptr = opnd2; 115 return(NOEXCEPTION); 116 } 117 /* 118 * return quiet NaN 119 */ 120 *dstptr = opnd1; 121 return(NOEXCEPTION); 122 } 123 } 124 /* 125 * check second operand for NaN's or infinity 126 */ 127 if (Sgl_isinfinity_exponent(opnd2)) { 128 if (Sgl_iszero_mantissa(opnd2)) { 129 if (Sgl_iszero_exponentmantissa(opnd1)) { 130 /* invalid since operands are zero & infinity */ 131 if (Is_invalidtrap_enabled()) 132 return(INVALIDEXCEPTION); 133 Set_invalidflag(); 134 Sgl_makequietnan(opnd2); 135 *dstptr = opnd2; 136 return(NOEXCEPTION); 137 } 138 /* 139 * return infinity 140 */ 141 Sgl_setinfinity_exponentmantissa(result); 142 *dstptr = result; 143 return(NOEXCEPTION); 144 } 145 /* 146 * is NaN; signaling or quiet? 147 */ 148 if (Sgl_isone_signaling(opnd2)) { 149 /* trap if INVALIDTRAP enabled */ 150 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 151 152 /* make NaN quiet */ 153 Set_invalidflag(); 154 Sgl_set_quiet(opnd2); 155 } 156 /* 157 * return quiet NaN 158 */ 159 *dstptr = opnd2; 160 return(NOEXCEPTION); 161 } 162 /* 163 * Generate exponent 164 */ 165 dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS; 166 167 /* 168 * Generate mantissa 169 */ 170 if (Sgl_isnotzero_exponent(opnd1)) { 171 /* set hidden bit */ 172 Sgl_clear_signexponent_set_hidden(opnd1); 173 } 174 else { 175 /* check for zero */ 176 if (Sgl_iszero_mantissa(opnd1)) { 177 Sgl_setzero_exponentmantissa(result); 178 *dstptr = result; 179 return(NOEXCEPTION); 180 } 181 /* is denormalized, adjust exponent */ 182 Sgl_clear_signexponent(opnd1); 183 Sgl_leftshiftby1(opnd1); 184 Sgl_normalize(opnd1,dest_exponent); 185 } 186 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 187 if (Sgl_isnotzero_exponent(opnd2)) { 188 Sgl_clear_signexponent_set_hidden(opnd2); 189 } 190 else { 191 /* check for zero */ 192 if (Sgl_iszero_mantissa(opnd2)) { 193 Sgl_setzero_exponentmantissa(result); 194 *dstptr = result; 195 return(NOEXCEPTION); 196 } 197 /* is denormalized; want to normalize */ 198 Sgl_clear_signexponent(opnd2); 199 Sgl_leftshiftby1(opnd2); 200 Sgl_normalize(opnd2,dest_exponent); 201 } 202 203 /* Multiply two source mantissas together */ 204 205 Sgl_leftshiftby4(opnd2); /* make room for guard bits */ 206 Sgl_setzero(opnd3); 207 /* 208 * Four bits at a time are inspected in each loop, and a 209 * simple shift and add multiply algorithm is used. 210 */ 211 for (count=1;count<SGL_P;count+=4) { 212 stickybit |= Slow4(opnd3); 213 Sgl_rightshiftby4(opnd3); 214 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3); 215 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2); 216 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1); 217 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2); 218 Sgl_rightshiftby4(opnd1); 219 } 220 /* make sure result is left-justified */ 221 if (Sgl_iszero_sign(opnd3)) { 222 Sgl_leftshiftby1(opnd3); 223 } 224 else { 225 /* result mantissa >= 2. */ 226 dest_exponent++; 227 } 228 /* check for denormalized result */ 229 while (Sgl_iszero_sign(opnd3)) { 230 Sgl_leftshiftby1(opnd3); 231 dest_exponent--; 232 } 233 /* 234 * check for guard, sticky and inexact bits 235 */ 236 stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1); 237 guardbit = Sbit24(opnd3); 238 inexact = guardbit | stickybit; 239 240 /* re-align mantissa */ 241 Sgl_rightshiftby8(opnd3); 242 243 /* 244 * round result 245 */ 246 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 247 Sgl_clear_signexponent(opnd3); 248 switch (Rounding_mode()) { 249 case ROUNDPLUS: 250 if (Sgl_iszero_sign(result)) 251 Sgl_increment(opnd3); 252 break; 253 case ROUNDMINUS: 254 if (Sgl_isone_sign(result)) 255 Sgl_increment(opnd3); 256 break; 257 case ROUNDNEAREST: 258 if (guardbit && 259 (stickybit || Sgl_isone_lowmantissa(opnd3))) 260 Sgl_increment(opnd3); 261 break; 262 } 263 if (Sgl_isone_hidden(opnd3)) dest_exponent++; 264 } 265 Sgl_set_mantissa(result,opnd3); 266 267 /* 268 * Test for overflow 269 */ 270 if (dest_exponent >= SGL_INFINITY_EXPONENT) { 271 /* trap if OVERFLOWTRAP enabled */ 272 if (Is_overflowtrap_enabled()) { 273 /* 274 * Adjust bias of result 275 */ 276 Sgl_setwrapped_exponent(result,dest_exponent,ovfl); 277 *dstptr = result; 278 if (inexact) { 279 if (Is_inexacttrap_enabled()) 280 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 281 else Set_inexactflag(); 282 } 283 return(OVERFLOWEXCEPTION); 284 } 285 inexact = TRUE; 286 Set_overflowflag(); 287 /* set result to infinity or largest number */ 288 Sgl_setoverflow(result); 289 } 290 /* 291 * Test for underflow 292 */ 293 else if (dest_exponent <= 0) { 294 /* trap if UNDERFLOWTRAP enabled */ 295 if (Is_underflowtrap_enabled()) { 296 /* 297 * Adjust bias of result 298 */ 299 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 300 *dstptr = result; 301 if (inexact) { 302 if (Is_inexacttrap_enabled()) 303 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 304 else Set_inexactflag(); 305 } 306 return(UNDERFLOWEXCEPTION); 307 } 308 309 /* Determine if should set underflow flag */ 310 is_tiny = TRUE; 311 if (dest_exponent == 0 && inexact) { 312 switch (Rounding_mode()) { 313 case ROUNDPLUS: 314 if (Sgl_iszero_sign(result)) { 315 Sgl_increment(opnd3); 316 if (Sgl_isone_hiddenoverflow(opnd3)) 317 is_tiny = FALSE; 318 Sgl_decrement(opnd3); 319 } 320 break; 321 case ROUNDMINUS: 322 if (Sgl_isone_sign(result)) { 323 Sgl_increment(opnd3); 324 if (Sgl_isone_hiddenoverflow(opnd3)) 325 is_tiny = FALSE; 326 Sgl_decrement(opnd3); 327 } 328 break; 329 case ROUNDNEAREST: 330 if (guardbit && (stickybit || 331 Sgl_isone_lowmantissa(opnd3))) { 332 Sgl_increment(opnd3); 333 if (Sgl_isone_hiddenoverflow(opnd3)) 334 is_tiny = FALSE; 335 Sgl_decrement(opnd3); 336 } 337 break; 338 } 339 } 340 341 /* 342 * denormalize result or set to signed zero 343 */ 344 stickybit = inexact; 345 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact); 346 347 /* return zero or smallest number */ 348 if (inexact) { 349 switch (Rounding_mode()) { 350 case ROUNDPLUS: 351 if (Sgl_iszero_sign(result)) 352 Sgl_increment(opnd3); 353 break; 354 case ROUNDMINUS: 355 if (Sgl_isone_sign(result)) 356 Sgl_increment(opnd3); 357 break; 358 case ROUNDNEAREST: 359 if (guardbit && (stickybit || 360 Sgl_isone_lowmantissa(opnd3))) 361 Sgl_increment(opnd3); 362 break; 363 } 364 if (is_tiny) Set_underflowflag(); 365 } 366 Sgl_set_exponentmantissa(result,opnd3); 367 } 368 else Sgl_set_exponent(result,dest_exponent); 369 *dstptr = result; 370 371 /* check for inexact */ 372 if (inexact) { 373 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 374 else Set_inexactflag(); 375 } 376 return(NOEXCEPTION); 377 } 378