Home | History | Annotate | Line # | Download | only in spmath
sfsqrt.c revision 1.2 
      1 /*	$NetBSD: sfsqrt.c,v 1.2 2003/07/15 02:29:43 lukem Exp $	*/
      2 
      3 /*	$OpenBSD: sfsqrt.c,v 1.5 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: sfsqrt.c,v 1.2 2003/07/15 02:29:43 lukem Exp $");
     46 
     47 #include "../spmath/float.h"
     48 #include "../spmath/sgl_float.h"
     49 
     50 /*
     51  *  Single Floating-point Square Root
     52  */
     53 
     54 /*ARGSUSED*/
     55 int
     56 sgl_fsqrt(srcptr,dstptr,status)
     57 
     58 sgl_floating_point *srcptr, *dstptr;
     59 unsigned int *status;
     60 {
     61 	register unsigned int src, result;
     62 	register int src_exponent, newbit, sum;
     63 	register int guardbit = FALSE, even_exponent;
     64 
     65 	src = *srcptr;
     66 	/*
     67 	 * check source operand for NaN or infinity
     68 	 */
     69 	if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
     70 		/*
     71 		 * is signaling NaN?
     72 		 */
     73 		if (Sgl_isone_signaling(src)) {
     74 			/* trap if INVALIDTRAP enabled */
     75 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
     76 			/* make NaN quiet */
     77 			Set_invalidflag();
     78 			Sgl_set_quiet(src);
     79 		}
     80 		/*
     81 		 * Return quiet NaN or positive infinity.
     82 		 *  Fall thru to negative test if negative infinity.
     83 		 */
     84 		if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
     85 			*dstptr = src;
     86 			return(NOEXCEPTION);
     87 		}
     88 	}
     89 
     90 	/*
     91 	 * check for zero source operand
     92 	 */
     93 	if (Sgl_iszero_exponentmantissa(src)) {
     94 		*dstptr = src;
     95 		return(NOEXCEPTION);
     96 	}
     97 
     98 	/*
     99 	 * check for negative source operand
    100 	 */
    101 	if (Sgl_isone_sign(src)) {
    102 		/* trap if INVALIDTRAP enabled */
    103 		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
    104 		/* make NaN quiet */
    105 		Set_invalidflag();
    106 		Sgl_makequietnan(src);
    107 		*dstptr = src;
    108 		return(NOEXCEPTION);
    109 	}
    110 
    111 	/*
    112 	 * Generate result
    113 	 */
    114 	if (src_exponent > 0) {
    115 		even_exponent = Sgl_hidden(src);
    116 		Sgl_clear_signexponent_set_hidden(src);
    117 	}
    118 	else {
    119 		/* normalize operand */
    120 		Sgl_clear_signexponent(src);
    121 		src_exponent++;
    122 		Sgl_normalize(src,src_exponent);
    123 		even_exponent = src_exponent & 1;
    124 	}
    125 	if (even_exponent) {
    126 		/* exponent is even */
    127 		/* Add comment here.  Explain why odd exponent needs correction */
    128 		Sgl_leftshiftby1(src);
    129 	}
    130 	/*
    131 	 * Add comment here.  Explain following algorithm.
    132 	 *
    133 	 * Trust me, it works.
    134 	 *
    135 	 */
    136 	Sgl_setzero(result);
    137 	newbit = 1 << SGL_P;
    138 	while (newbit && Sgl_isnotzero(src)) {
    139 		Sgl_addition(result,newbit,sum);
    140 		if(sum <= Sgl_all(src)) {
    141 			/* update result */
    142 			Sgl_addition(result,(newbit<<1),result);
    143 			Sgl_subtract(src,sum,src);
    144 		}
    145 		Sgl_rightshiftby1(newbit);
    146 		Sgl_leftshiftby1(src);
    147 	}
    148 	/* correct exponent for pre-shift */
    149 	if (even_exponent) {
    150 		Sgl_rightshiftby1(result);
    151 	}
    152 
    153 	/* check for inexact */
    154 	if (Sgl_isnotzero(src)) {
    155 		if (!even_exponent & Sgl_islessthan(result,src))
    156 			Sgl_increment(result);
    157 		guardbit = Sgl_lowmantissa(result);
    158 		Sgl_rightshiftby1(result);
    159 
    160 		/*  now round result  */
    161 		switch (Rounding_mode()) {
    162 		case ROUNDPLUS:
    163 		     Sgl_increment(result);
    164 		     break;
    165 		case ROUNDNEAREST:
    166 		     /* stickybit is always true, so guardbit
    167 		      * is enough to determine rounding */
    168 		     if (guardbit) {
    169 			Sgl_increment(result);
    170 		     }
    171 		     break;
    172 		}
    173 		/* increment result exponent by 1 if mantissa overflowed */
    174 		if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
    175 
    176 		if (Is_inexacttrap_enabled()) {
    177 			Sgl_set_exponent(result,
    178 			 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
    179 			*dstptr = result;
    180 			return(INEXACTEXCEPTION);
    181 		}
    182 		else Set_inexactflag();
    183 	}
    184 	else {
    185 		Sgl_rightshiftby1(result);
    186 	}
    187 	Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
    188 	*dstptr = result;
    189 	return(NOEXCEPTION);
    190 }
    191