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