1 /* $NetBSD: sfsqrt.c,v 1.4 2007/02/22 05:46:31 thorpej 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.4 2007/02/22 05:46:31 thorpej 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