1 /* $NetBSD: dfsqrt.c,v 1.1 2002/06/05 01:04:24 fredette Exp $ */ 2 3 /* $OpenBSD: dfsqrt.c,v 1.5 2001/03/29 03:58:17 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/dbl_float.h" 46 47 /* 48 * Double Floating-point Square Root 49 */ 50 51 /*ARGSUSED*/ 52 int 53 dbl_fsqrt(srcptr,dstptr,status) 54 55 dbl_floating_point *srcptr, *dstptr; 56 unsigned int *status; 57 { 58 register unsigned int srcp1, srcp2, resultp1, resultp2; 59 register unsigned int newbitp1, newbitp2, sump1, sump2; 60 register int src_exponent; 61 register int guardbit = FALSE, even_exponent; 62 63 Dbl_copyfromptr(srcptr,srcp1,srcp2); 64 /* 65 * check source operand for NaN or infinity 66 */ 67 if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 68 /* 69 * is signaling NaN? 70 */ 71 if (Dbl_isone_signaling(srcp1)) { 72 /* trap if INVALIDTRAP enabled */ 73 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 74 /* make NaN quiet */ 75 Set_invalidflag(); 76 Dbl_set_quiet(srcp1); 77 } 78 /* 79 * Return quiet NaN or positive infinity. 80 * Fall thru to negative test if negative infinity. 81 */ 82 if (Dbl_iszero_sign(srcp1) || 83 Dbl_isnotzero_mantissa(srcp1,srcp2)) { 84 Dbl_copytoptr(srcp1,srcp2,dstptr); 85 return(NOEXCEPTION); 86 } 87 } 88 89 /* 90 * check for zero source operand 91 */ 92 if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 93 Dbl_copytoptr(srcp1,srcp2,dstptr); 94 return(NOEXCEPTION); 95 } 96 97 /* 98 * check for negative source operand 99 */ 100 if (Dbl_isone_sign(srcp1)) { 101 /* trap if INVALIDTRAP enabled */ 102 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 103 /* make NaN quiet */ 104 Set_invalidflag(); 105 Dbl_makequietnan(srcp1,srcp2); 106 Dbl_copytoptr(srcp1,srcp2,dstptr); 107 return(NOEXCEPTION); 108 } 109 110 /* 111 * Generate result 112 */ 113 if (src_exponent > 0) { 114 even_exponent = Dbl_hidden(srcp1); 115 Dbl_clear_signexponent_set_hidden(srcp1); 116 } 117 else { 118 /* normalize operand */ 119 Dbl_clear_signexponent(srcp1); 120 src_exponent++; 121 Dbl_normalize(srcp1,srcp2,src_exponent); 122 even_exponent = src_exponent & 1; 123 } 124 if (even_exponent) { 125 /* exponent is even */ 126 /* Add comment here. Explain why odd exponent needs correction */ 127 Dbl_leftshiftby1(srcp1,srcp2); 128 } 129 /* 130 * Add comment here. Explain following algorithm. 131 * 132 * Trust me, it works. 133 * 134 */ 135 Dbl_setzero(resultp1,resultp2); 136 Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 137 Dbl_setzero_mantissap2(newbitp2); 138 while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 139 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 140 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 141 Dbl_leftshiftby1(newbitp1,newbitp2); 142 /* update result */ 143 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 144 resultp1,resultp2); 145 Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 146 Dbl_rightshiftby2(newbitp1,newbitp2); 147 } 148 else { 149 Dbl_rightshiftby1(newbitp1,newbitp2); 150 } 151 Dbl_leftshiftby1(srcp1,srcp2); 152 } 153 /* correct exponent for pre-shift */ 154 if (even_exponent) { 155 Dbl_rightshiftby1(resultp1,resultp2); 156 } 157 158 /* check for inexact */ 159 if (Dbl_isnotzero(srcp1,srcp2)) { 160 if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 161 Dbl_increment(resultp1,resultp2); 162 } 163 guardbit = Dbl_lowmantissap2(resultp2); 164 Dbl_rightshiftby1(resultp1,resultp2); 165 166 /* now round result */ 167 switch (Rounding_mode()) { 168 case ROUNDPLUS: 169 Dbl_increment(resultp1,resultp2); 170 break; 171 case ROUNDNEAREST: 172 /* stickybit is always true, so guardbit 173 * is enough to determine rounding */ 174 if (guardbit) { 175 Dbl_increment(resultp1,resultp2); 176 } 177 break; 178 } 179 /* increment result exponent by 1 if mantissa overflowed */ 180 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 181 182 if (Is_inexacttrap_enabled()) { 183 Dbl_set_exponent(resultp1, 184 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 185 Dbl_copytoptr(resultp1,resultp2,dstptr); 186 return(INEXACTEXCEPTION); 187 } 188 else Set_inexactflag(); 189 } 190 else { 191 Dbl_rightshiftby1(resultp1,resultp2); 192 } 193 Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 194 Dbl_copytoptr(resultp1,resultp2,dstptr); 195 return(NOEXCEPTION); 196 } 197