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