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