Home | History | Annotate | Line # | Download | only in spmath
sfdiv.c revision 1.2 
      1 /*	$NetBSD: sfdiv.c,v 1.2 2003/07/15 02:29:42 lukem Exp $	*/
      2 
      3 /*	$OpenBSD: sfdiv.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: sfdiv.c,v 1.2 2003/07/15 02:29:42 lukem Exp $");
     46 
     47 #include "../spmath/float.h"
     48 #include "../spmath/sgl_float.h"
     49 
     50 /*
     51  *  Single Precision Floating-point Divide
     52  */
     53 int
     54 sgl_fdiv(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_isinfinity(opnd2)) {
     78 					/*
     79 					 * invalid since both operands
     80 					 * are infinity
     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 			/*
    135 			 * return zero
    136 			 */
    137 			Sgl_setzero_exponentmantissa(result);
    138 			*dstptr = result;
    139 			return(NOEXCEPTION);
    140 		}
    141 		/*
    142 		 * is NaN; signaling or quiet?
    143 		 */
    144 		if (Sgl_isone_signaling(opnd2)) {
    145 			/* trap if INVALIDTRAP enabled */
    146 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
    147 			/* make NaN quiet */
    148 			Set_invalidflag();
    149 			Sgl_set_quiet(opnd2);
    150 		}
    151 		/*
    152 		 * return quiet NaN
    153 		 */
    154 		*dstptr = opnd2;
    155 		return(NOEXCEPTION);
    156 	}
    157 	/*
    158 	 * check for division by zero
    159 	 */
    160 	if (Sgl_iszero_exponentmantissa(opnd2)) {
    161 		if (Sgl_iszero_exponentmantissa(opnd1)) {
    162 			/* invalid since both operands are zero */
    163 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
    164 			Set_invalidflag();
    165 			Sgl_makequietnan(result);
    166 			*dstptr = result;
    167 			return(NOEXCEPTION);
    168 		}
    169 		if (Is_divisionbyzerotrap_enabled())
    170 			return(DIVISIONBYZEROEXCEPTION);
    171 		Set_divisionbyzeroflag();
    172 		Sgl_setinfinity_exponentmantissa(result);
    173 		*dstptr = result;
    174 		return(NOEXCEPTION);
    175 	}
    176 	/*
    177 	 * Generate exponent
    178 	 */
    179 	dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
    180 
    181 	/*
    182 	 * Generate mantissa
    183 	 */
    184 	if (Sgl_isnotzero_exponent(opnd1)) {
    185 		/* set hidden bit */
    186 		Sgl_clear_signexponent_set_hidden(opnd1);
    187 	}
    188 	else {
    189 		/* check for zero */
    190 		if (Sgl_iszero_mantissa(opnd1)) {
    191 			Sgl_setzero_exponentmantissa(result);
    192 			*dstptr = result;
    193 			return(NOEXCEPTION);
    194 		}
    195 		/* is denormalized; want to normalize */
    196 		Sgl_clear_signexponent(opnd1);
    197 		Sgl_leftshiftby1(opnd1);
    198 		Sgl_normalize(opnd1,dest_exponent);
    199 	}
    200 	/* opnd2 needs to have hidden bit set with msb in hidden bit */
    201 	if (Sgl_isnotzero_exponent(opnd2)) {
    202 		Sgl_clear_signexponent_set_hidden(opnd2);
    203 	}
    204 	else {
    205 		/* is denormalized; want to normalize */
    206 		Sgl_clear_signexponent(opnd2);
    207 		Sgl_leftshiftby1(opnd2);
    208 		while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
    209 			Sgl_leftshiftby8(opnd2);
    210 			dest_exponent += 8;
    211 		}
    212 		if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
    213 			Sgl_leftshiftby4(opnd2);
    214 			dest_exponent += 4;
    215 		}
    216 		while(Sgl_iszero_hidden(opnd2)) {
    217 			Sgl_leftshiftby1(opnd2);
    218 			dest_exponent += 1;
    219 		}
    220 	}
    221 
    222 	/* Divide the source mantissas */
    223 
    224 	/*
    225 	 * A non_restoring divide algorithm is used.
    226 	 */
    227 	Sgl_subtract(opnd1,opnd2,opnd1);
    228 	Sgl_setzero(opnd3);
    229 	for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
    230 		Sgl_leftshiftby1(opnd1);
    231 		Sgl_leftshiftby1(opnd3);
    232 		if (Sgl_iszero_sign(opnd1)) {
    233 			Sgl_setone_lowmantissa(opnd3);
    234 			Sgl_subtract(opnd1,opnd2,opnd1);
    235 		}
    236 		else Sgl_addition(opnd1,opnd2,opnd1);
    237 	}
    238 	if (count <= SGL_P) {
    239 		Sgl_leftshiftby1(opnd3);
    240 		Sgl_setone_lowmantissa(opnd3);
    241 		Sgl_leftshift(opnd3,SGL_P-count);
    242 		if (Sgl_iszero_hidden(opnd3)) {
    243 			Sgl_leftshiftby1(opnd3);
    244 			dest_exponent--;
    245 		}
    246 	}
    247 	else {
    248 		if (Sgl_iszero_hidden(opnd3)) {
    249 			/* need to get one more bit of result */
    250 			Sgl_leftshiftby1(opnd1);
    251 			Sgl_leftshiftby1(opnd3);
    252 			if (Sgl_iszero_sign(opnd1)) {
    253 				Sgl_setone_lowmantissa(opnd3);
    254 				Sgl_subtract(opnd1,opnd2,opnd1);
    255 			}
    256 			else Sgl_addition(opnd1,opnd2,opnd1);
    257 			dest_exponent--;
    258 		}
    259 		if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
    260 		stickybit = Sgl_all(opnd1);
    261 	}
    262 	inexact = guardbit | stickybit;
    263 
    264 	/*
    265 	 * round result
    266 	 */
    267 	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
    268 		Sgl_clear_signexponent(opnd3);
    269 		switch (Rounding_mode()) {
    270 			case ROUNDPLUS:
    271 				if (Sgl_iszero_sign(result))
    272 					Sgl_increment_mantissa(opnd3);
    273 				break;
    274 			case ROUNDMINUS:
    275 				if (Sgl_isone_sign(result))
    276 					Sgl_increment_mantissa(opnd3);
    277 				break;
    278 			case ROUNDNEAREST:
    279 				if (guardbit &&
    280 				    (stickybit || Sgl_isone_lowmantissa(opnd3)))
    281 					Sgl_increment_mantissa(opnd3);
    282 		}
    283 		if (Sgl_isone_hidden(opnd3)) dest_exponent++;
    284 	}
    285 	Sgl_set_mantissa(result,opnd3);
    286 
    287 	/*
    288 	 * Test for overflow
    289 	 */
    290 	if (dest_exponent >= SGL_INFINITY_EXPONENT) {
    291 		/* trap if OVERFLOWTRAP enabled */
    292 		if (Is_overflowtrap_enabled()) {
    293 			/*
    294 			 * Adjust bias of result
    295 			 */
    296 			Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
    297 			*dstptr = result;
    298 			if (inexact) {
    299 			    if (Is_inexacttrap_enabled())
    300 				return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
    301 			    else Set_inexactflag();
    302 			}
    303 			return(OVERFLOWEXCEPTION);
    304 		}
    305 		Set_overflowflag();
    306 		/* set result to infinity or largest number */
    307 		Sgl_setoverflow(result);
    308 		inexact = TRUE;
    309 	}
    310 	/*
    311 	 * Test for underflow
    312 	 */
    313 	else if (dest_exponent <= 0) {
    314 		/* trap if UNDERFLOWTRAP enabled */
    315 		if (Is_underflowtrap_enabled()) {
    316 			/*
    317 			 * Adjust bias of result
    318 			 */
    319 			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
    320 			*dstptr = result;
    321 			if (inexact) {
    322 			    if (Is_inexacttrap_enabled())
    323 				return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
    324 			    else Set_inexactflag();
    325 			}
    326 			return(UNDERFLOWEXCEPTION);
    327 		}
    328 
    329 		/* Determine if should set underflow flag */
    330 		is_tiny = TRUE;
    331 		if (dest_exponent == 0 && inexact) {
    332 			switch (Rounding_mode()) {
    333 			case ROUNDPLUS:
    334 				if (Sgl_iszero_sign(result)) {
    335 					Sgl_increment(opnd3);
    336 					if (Sgl_isone_hiddenoverflow(opnd3))
    337 						is_tiny = FALSE;
    338 					Sgl_decrement(opnd3);
    339 				}
    340 				break;
    341 			case ROUNDMINUS:
    342 				if (Sgl_isone_sign(result)) {
    343 					Sgl_increment(opnd3);
    344 					if (Sgl_isone_hiddenoverflow(opnd3))
    345 						is_tiny = FALSE;
    346 					Sgl_decrement(opnd3);
    347 				}
    348 				break;
    349 			case ROUNDNEAREST:
    350 				if (guardbit && (stickybit ||
    351 				    Sgl_isone_lowmantissa(opnd3))) {
    352 					Sgl_increment(opnd3);
    353 					if (Sgl_isone_hiddenoverflow(opnd3))
    354 						is_tiny = FALSE;
    355 					Sgl_decrement(opnd3);
    356 				}
    357 				break;
    358 			}
    359 		}
    360 
    361 		/*
    362 		 * denormalize result or set to signed zero
    363 		 */
    364 		stickybit = inexact;
    365 		Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
    366 
    367 		/* return rounded number */
    368 		if (inexact) {
    369 			switch (Rounding_mode()) {
    370 			case ROUNDPLUS:
    371 				if (Sgl_iszero_sign(result)) {
    372 					Sgl_increment(opnd3);
    373 				}
    374 				break;
    375 			case ROUNDMINUS:
    376 				if (Sgl_isone_sign(result))  {
    377 					Sgl_increment(opnd3);
    378 				}
    379 				break;
    380 			case ROUNDNEAREST:
    381 				if (guardbit && (stickybit ||
    382 				    Sgl_isone_lowmantissa(opnd3))) {
    383 					Sgl_increment(opnd3);
    384 				}
    385 				break;
    386 			}
    387 			if (is_tiny) Set_underflowflag();
    388 		}
    389 		Sgl_set_exponentmantissa(result,opnd3);
    390 	}
    391 	else Sgl_set_exponent(result,dest_exponent);
    392 	*dstptr = result;
    393 	/* check for inexact */
    394 	if (inexact) {
    395 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
    396 		else  Set_inexactflag();
    397 	}
    398 	return(NOEXCEPTION);
    399 }
    400