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