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