Home | History | Annotate | Line # | Download | only in spmath
sfmpy.c revision 1.4 
      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