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