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