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