1 /* $OpenBSD: sfsqrt.c,v 1.4 2000/01/11 08:18:43 mickey Exp $ */ 2 3 /* 4 * Copyright 1996 1995 by Open Software Foundation, Inc. 5 * All Rights Reserved 6 * 7 * Permission to use, copy, modify, and distribute this software and 8 * its documentation for any purpose and without fee is hereby granted, 9 * provided that the above copyright notice appears in all copies and 10 * that both the copyright notice and this permission notice appear in 11 * supporting documentation. 12 * 13 * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE 14 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 15 * FOR A PARTICULAR PURPOSE. 16 * 17 * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR 18 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM 19 * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT, 20 * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION 21 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 22 * 23 */ 24 /* 25 * pmk1.1 26 */ 27 /* 28 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 29 * 30 * To anyone who acknowledges that this file is provided "AS IS" 31 * without any express or implied warranty: 32 * permission to use, copy, modify, and distribute this file 33 * for any purpose is hereby granted without fee, provided that 34 * the above copyright notice and this notice appears in all 35 * copies, and that the name of Hewlett-Packard Company not be 36 * used in advertising or publicity pertaining to distribution 37 * of the software without specific, written prior permission. 38 * Hewlett-Packard Company makes no representations about the 39 * suitability of this software for any purpose. 40 */ 41 42 #include "../spmath/float.h" 43 #include "../spmath/sgl_float.h" 44 45 /* 46 * Single Floating-point Square Root 47 */ 48 49 /*ARGSUSED*/ 50 int 51 sgl_fsqrt(srcptr,dstptr,status) 52 53 sgl_floating_point *srcptr, *dstptr; 54 unsigned int *status; 55 { 56 register unsigned int src, result; 57 register int src_exponent, newbit, sum; 58 register boolean guardbit = FALSE, even_exponent; 59 60 src = *srcptr; 61 /* 62 * check source operand for NaN or infinity 63 */ 64 if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) { 65 /* 66 * is signaling NaN? 67 */ 68 if (Sgl_isone_signaling(src)) { 69 /* trap if INVALIDTRAP enabled */ 70 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 71 /* make NaN quiet */ 72 Set_invalidflag(); 73 Sgl_set_quiet(src); 74 } 75 /* 76 * Return quiet NaN or positive infinity. 77 * Fall thru to negative test if negative infinity. 78 */ 79 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) { 80 *dstptr = src; 81 return(NOEXCEPTION); 82 } 83 } 84 85 /* 86 * check for zero source operand 87 */ 88 if (Sgl_iszero_exponentmantissa(src)) { 89 *dstptr = src; 90 return(NOEXCEPTION); 91 } 92 93 /* 94 * check for negative source operand 95 */ 96 if (Sgl_isone_sign(src)) { 97 /* trap if INVALIDTRAP enabled */ 98 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 99 /* make NaN quiet */ 100 Set_invalidflag(); 101 Sgl_makequietnan(src); 102 *dstptr = src; 103 return(NOEXCEPTION); 104 } 105 106 /* 107 * Generate result 108 */ 109 if (src_exponent > 0) { 110 even_exponent = Sgl_hidden(src); 111 Sgl_clear_signexponent_set_hidden(src); 112 } 113 else { 114 /* normalize operand */ 115 Sgl_clear_signexponent(src); 116 src_exponent++; 117 Sgl_normalize(src,src_exponent); 118 even_exponent = src_exponent & 1; 119 } 120 if (even_exponent) { 121 /* exponent is even */ 122 /* Add comment here. Explain why odd exponent needs correction */ 123 Sgl_leftshiftby1(src); 124 } 125 /* 126 * Add comment here. Explain following algorithm. 127 * 128 * Trust me, it works. 129 * 130 */ 131 Sgl_setzero(result); 132 newbit = 1 << SGL_P; 133 while (newbit && Sgl_isnotzero(src)) { 134 Sgl_addition(result,newbit,sum); 135 if(sum <= Sgl_all(src)) { 136 /* update result */ 137 Sgl_addition(result,(newbit<<1),result); 138 Sgl_subtract(src,sum,src); 139 } 140 Sgl_rightshiftby1(newbit); 141 Sgl_leftshiftby1(src); 142 } 143 /* correct exponent for pre-shift */ 144 if (even_exponent) { 145 Sgl_rightshiftby1(result); 146 } 147 148 /* check for inexact */ 149 if (Sgl_isnotzero(src)) { 150 if (!even_exponent & Sgl_islessthan(result,src)) 151 Sgl_increment(result); 152 guardbit = Sgl_lowmantissa(result); 153 Sgl_rightshiftby1(result); 154 155 /* now round result */ 156 switch (Rounding_mode()) { 157 case ROUNDPLUS: 158 Sgl_increment(result); 159 break; 160 case ROUNDNEAREST: 161 /* stickybit is always true, so guardbit 162 * is enough to determine rounding */ 163 if (guardbit) { 164 Sgl_increment(result); 165 } 166 break; 167 } 168 /* increment result exponent by 1 if mantissa overflowed */ 169 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2; 170 171 if (Is_inexacttrap_enabled()) { 172 Sgl_set_exponent(result, 173 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 174 *dstptr = result; 175 return(INEXACTEXCEPTION); 176 } 177 else Set_inexactflag(); 178 } 179 else { 180 Sgl_rightshiftby1(result); 181 } 182 Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 183 *dstptr = result; 184 return(NOEXCEPTION); 185 } 186