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