dfsqrt.c revision 1.1.2.2
1 1.1.2.2 jdolecek /* $NetBSD: dfsqrt.c,v 1.1.2.2 2002/06/23 17:37:12 jdolecek Exp $ */ 2 1.1.2.2 jdolecek 3 1.1.2.2 jdolecek /* $OpenBSD: dfsqrt.c,v 1.5 2001/03/29 03:58:17 mickey Exp $ */ 4 1.1.2.2 jdolecek 5 1.1.2.2 jdolecek /* 6 1.1.2.2 jdolecek * Copyright 1996 1995 by Open Software Foundation, Inc. 7 1.1.2.2 jdolecek * All Rights Reserved 8 1.1.2.2 jdolecek * 9 1.1.2.2 jdolecek * Permission to use, copy, modify, and distribute this software and 10 1.1.2.2 jdolecek * its documentation for any purpose and without fee is hereby granted, 11 1.1.2.2 jdolecek * provided that the above copyright notice appears in all copies and 12 1.1.2.2 jdolecek * that both the copyright notice and this permission notice appear in 13 1.1.2.2 jdolecek * supporting documentation. 14 1.1.2.2 jdolecek * 15 1.1.2.2 jdolecek * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE 16 1.1.2.2 jdolecek * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 17 1.1.2.2 jdolecek * FOR A PARTICULAR PURPOSE. 18 1.1.2.2 jdolecek * 19 1.1.2.2 jdolecek * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR 20 1.1.2.2 jdolecek * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM 21 1.1.2.2 jdolecek * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT, 22 1.1.2.2 jdolecek * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION 23 1.1.2.2 jdolecek * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 24 1.1.2.2 jdolecek * 25 1.1.2.2 jdolecek */ 26 1.1.2.2 jdolecek /* 27 1.1.2.2 jdolecek * pmk1.1 28 1.1.2.2 jdolecek */ 29 1.1.2.2 jdolecek /* 30 1.1.2.2 jdolecek * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 31 1.1.2.2 jdolecek * 32 1.1.2.2 jdolecek * To anyone who acknowledges that this file is provided "AS IS" 33 1.1.2.2 jdolecek * without any express or implied warranty: 34 1.1.2.2 jdolecek * permission to use, copy, modify, and distribute this file 35 1.1.2.2 jdolecek * for any purpose is hereby granted without fee, provided that 36 1.1.2.2 jdolecek * the above copyright notice and this notice appears in all 37 1.1.2.2 jdolecek * copies, and that the name of Hewlett-Packard Company not be 38 1.1.2.2 jdolecek * used in advertising or publicity pertaining to distribution 39 1.1.2.2 jdolecek * of the software without specific, written prior permission. 40 1.1.2.2 jdolecek * Hewlett-Packard Company makes no representations about the 41 1.1.2.2 jdolecek * suitability of this software for any purpose. 42 1.1.2.2 jdolecek */ 43 1.1.2.2 jdolecek 44 1.1.2.2 jdolecek #include "../spmath/float.h" 45 1.1.2.2 jdolecek #include "../spmath/dbl_float.h" 46 1.1.2.2 jdolecek 47 1.1.2.2 jdolecek /* 48 1.1.2.2 jdolecek * Double Floating-point Square Root 49 1.1.2.2 jdolecek */ 50 1.1.2.2 jdolecek 51 1.1.2.2 jdolecek /*ARGSUSED*/ 52 1.1.2.2 jdolecek int 53 1.1.2.2 jdolecek dbl_fsqrt(srcptr,dstptr,status) 54 1.1.2.2 jdolecek 55 1.1.2.2 jdolecek dbl_floating_point *srcptr, *dstptr; 56 1.1.2.2 jdolecek unsigned int *status; 57 1.1.2.2 jdolecek { 58 1.1.2.2 jdolecek register unsigned int srcp1, srcp2, resultp1, resultp2; 59 1.1.2.2 jdolecek register unsigned int newbitp1, newbitp2, sump1, sump2; 60 1.1.2.2 jdolecek register int src_exponent; 61 1.1.2.2 jdolecek register int guardbit = FALSE, even_exponent; 62 1.1.2.2 jdolecek 63 1.1.2.2 jdolecek Dbl_copyfromptr(srcptr,srcp1,srcp2); 64 1.1.2.2 jdolecek /* 65 1.1.2.2 jdolecek * check source operand for NaN or infinity 66 1.1.2.2 jdolecek */ 67 1.1.2.2 jdolecek if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 68 1.1.2.2 jdolecek /* 69 1.1.2.2 jdolecek * is signaling NaN? 70 1.1.2.2 jdolecek */ 71 1.1.2.2 jdolecek if (Dbl_isone_signaling(srcp1)) { 72 1.1.2.2 jdolecek /* trap if INVALIDTRAP enabled */ 73 1.1.2.2 jdolecek if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 74 1.1.2.2 jdolecek /* make NaN quiet */ 75 1.1.2.2 jdolecek Set_invalidflag(); 76 1.1.2.2 jdolecek Dbl_set_quiet(srcp1); 77 1.1.2.2 jdolecek } 78 1.1.2.2 jdolecek /* 79 1.1.2.2 jdolecek * Return quiet NaN or positive infinity. 80 1.1.2.2 jdolecek * Fall thru to negative test if negative infinity. 81 1.1.2.2 jdolecek */ 82 1.1.2.2 jdolecek if (Dbl_iszero_sign(srcp1) || 83 1.1.2.2 jdolecek Dbl_isnotzero_mantissa(srcp1,srcp2)) { 84 1.1.2.2 jdolecek Dbl_copytoptr(srcp1,srcp2,dstptr); 85 1.1.2.2 jdolecek return(NOEXCEPTION); 86 1.1.2.2 jdolecek } 87 1.1.2.2 jdolecek } 88 1.1.2.2 jdolecek 89 1.1.2.2 jdolecek /* 90 1.1.2.2 jdolecek * check for zero source operand 91 1.1.2.2 jdolecek */ 92 1.1.2.2 jdolecek if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 93 1.1.2.2 jdolecek Dbl_copytoptr(srcp1,srcp2,dstptr); 94 1.1.2.2 jdolecek return(NOEXCEPTION); 95 1.1.2.2 jdolecek } 96 1.1.2.2 jdolecek 97 1.1.2.2 jdolecek /* 98 1.1.2.2 jdolecek * check for negative source operand 99 1.1.2.2 jdolecek */ 100 1.1.2.2 jdolecek if (Dbl_isone_sign(srcp1)) { 101 1.1.2.2 jdolecek /* trap if INVALIDTRAP enabled */ 102 1.1.2.2 jdolecek if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 103 1.1.2.2 jdolecek /* make NaN quiet */ 104 1.1.2.2 jdolecek Set_invalidflag(); 105 1.1.2.2 jdolecek Dbl_makequietnan(srcp1,srcp2); 106 1.1.2.2 jdolecek Dbl_copytoptr(srcp1,srcp2,dstptr); 107 1.1.2.2 jdolecek return(NOEXCEPTION); 108 1.1.2.2 jdolecek } 109 1.1.2.2 jdolecek 110 1.1.2.2 jdolecek /* 111 1.1.2.2 jdolecek * Generate result 112 1.1.2.2 jdolecek */ 113 1.1.2.2 jdolecek if (src_exponent > 0) { 114 1.1.2.2 jdolecek even_exponent = Dbl_hidden(srcp1); 115 1.1.2.2 jdolecek Dbl_clear_signexponent_set_hidden(srcp1); 116 1.1.2.2 jdolecek } 117 1.1.2.2 jdolecek else { 118 1.1.2.2 jdolecek /* normalize operand */ 119 1.1.2.2 jdolecek Dbl_clear_signexponent(srcp1); 120 1.1.2.2 jdolecek src_exponent++; 121 1.1.2.2 jdolecek Dbl_normalize(srcp1,srcp2,src_exponent); 122 1.1.2.2 jdolecek even_exponent = src_exponent & 1; 123 1.1.2.2 jdolecek } 124 1.1.2.2 jdolecek if (even_exponent) { 125 1.1.2.2 jdolecek /* exponent is even */ 126 1.1.2.2 jdolecek /* Add comment here. Explain why odd exponent needs correction */ 127 1.1.2.2 jdolecek Dbl_leftshiftby1(srcp1,srcp2); 128 1.1.2.2 jdolecek } 129 1.1.2.2 jdolecek /* 130 1.1.2.2 jdolecek * Add comment here. Explain following algorithm. 131 1.1.2.2 jdolecek * 132 1.1.2.2 jdolecek * Trust me, it works. 133 1.1.2.2 jdolecek * 134 1.1.2.2 jdolecek */ 135 1.1.2.2 jdolecek Dbl_setzero(resultp1,resultp2); 136 1.1.2.2 jdolecek Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 137 1.1.2.2 jdolecek Dbl_setzero_mantissap2(newbitp2); 138 1.1.2.2 jdolecek while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 139 1.1.2.2 jdolecek Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 140 1.1.2.2 jdolecek if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 141 1.1.2.2 jdolecek Dbl_leftshiftby1(newbitp1,newbitp2); 142 1.1.2.2 jdolecek /* update result */ 143 1.1.2.2 jdolecek Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 144 1.1.2.2 jdolecek resultp1,resultp2); 145 1.1.2.2 jdolecek Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 146 1.1.2.2 jdolecek Dbl_rightshiftby2(newbitp1,newbitp2); 147 1.1.2.2 jdolecek } 148 1.1.2.2 jdolecek else { 149 1.1.2.2 jdolecek Dbl_rightshiftby1(newbitp1,newbitp2); 150 1.1.2.2 jdolecek } 151 1.1.2.2 jdolecek Dbl_leftshiftby1(srcp1,srcp2); 152 1.1.2.2 jdolecek } 153 1.1.2.2 jdolecek /* correct exponent for pre-shift */ 154 1.1.2.2 jdolecek if (even_exponent) { 155 1.1.2.2 jdolecek Dbl_rightshiftby1(resultp1,resultp2); 156 1.1.2.2 jdolecek } 157 1.1.2.2 jdolecek 158 1.1.2.2 jdolecek /* check for inexact */ 159 1.1.2.2 jdolecek if (Dbl_isnotzero(srcp1,srcp2)) { 160 1.1.2.2 jdolecek if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 161 1.1.2.2 jdolecek Dbl_increment(resultp1,resultp2); 162 1.1.2.2 jdolecek } 163 1.1.2.2 jdolecek guardbit = Dbl_lowmantissap2(resultp2); 164 1.1.2.2 jdolecek Dbl_rightshiftby1(resultp1,resultp2); 165 1.1.2.2 jdolecek 166 1.1.2.2 jdolecek /* now round result */ 167 1.1.2.2 jdolecek switch (Rounding_mode()) { 168 1.1.2.2 jdolecek case ROUNDPLUS: 169 1.1.2.2 jdolecek Dbl_increment(resultp1,resultp2); 170 1.1.2.2 jdolecek break; 171 1.1.2.2 jdolecek case ROUNDNEAREST: 172 1.1.2.2 jdolecek /* stickybit is always true, so guardbit 173 1.1.2.2 jdolecek * is enough to determine rounding */ 174 1.1.2.2 jdolecek if (guardbit) { 175 1.1.2.2 jdolecek Dbl_increment(resultp1,resultp2); 176 1.1.2.2 jdolecek } 177 1.1.2.2 jdolecek break; 178 1.1.2.2 jdolecek } 179 1.1.2.2 jdolecek /* increment result exponent by 1 if mantissa overflowed */ 180 1.1.2.2 jdolecek if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 181 1.1.2.2 jdolecek 182 1.1.2.2 jdolecek if (Is_inexacttrap_enabled()) { 183 1.1.2.2 jdolecek Dbl_set_exponent(resultp1, 184 1.1.2.2 jdolecek ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 185 1.1.2.2 jdolecek Dbl_copytoptr(resultp1,resultp2,dstptr); 186 1.1.2.2 jdolecek return(INEXACTEXCEPTION); 187 1.1.2.2 jdolecek } 188 1.1.2.2 jdolecek else Set_inexactflag(); 189 1.1.2.2 jdolecek } 190 1.1.2.2 jdolecek else { 191 1.1.2.2 jdolecek Dbl_rightshiftby1(resultp1,resultp2); 192 1.1.2.2 jdolecek } 193 1.1.2.2 jdolecek Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 194 1.1.2.2 jdolecek Dbl_copytoptr(resultp1,resultp2,dstptr); 195 1.1.2.2 jdolecek return(NOEXCEPTION); 196 1.1.2.2 jdolecek } 197
Indexes created Thu Oct 16 14:10:15 GMT 2025