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