1 /* $OpenBSD: dfsqrt.c,v 1.7 2003/04/10 17:27:58 mickey Exp $ */ 2 /* 3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY 4 To anyone who acknowledges that this file is provided "AS IS" 5 without any express or implied warranty: 6 permission to use, copy, modify, and distribute this file 7 for any purpose is hereby granted without fee, provided that 8 the above copyright notice and this notice appears in all 9 copies, and that the name of Hewlett-Packard Company not be 10 used in advertising or publicity pertaining to distribution 11 of the software without specific, written prior permission. 12 Hewlett-Packard Company makes no representations about the 13 suitability of this software for any purpose. 14 */ 15 /* @(#)dfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:05:46 */ 16 17 #include "float.h" 18 #include "dbl_float.h" 19 20 /* 21 * Double Floating-point Square Root 22 */ 23 24 /*ARGSUSED*/ 25 int 26 dbl_fsqrt(srcptr, null, dstptr, status) 27 dbl_floating_point *srcptr, *null, *dstptr; 28 unsigned int *status; 29 { 30 register unsigned int srcp1, srcp2, resultp1, resultp2; 31 register unsigned int newbitp1, newbitp2, sump1, sump2; 32 register int src_exponent; 33 register int guardbit = FALSE, even_exponent; 34 35 Dbl_copyfromptr(srcptr,srcp1,srcp2); 36 /* 37 * check source operand for NaN or infinity 38 */ 39 if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 40 /* 41 * is signaling NaN? 42 */ 43 if (Dbl_isone_signaling(srcp1)) { 44 /* trap if INVALIDTRAP enabled */ 45 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 46 /* make NaN quiet */ 47 Set_invalidflag(); 48 Dbl_set_quiet(srcp1); 49 } 50 /* 51 * Return quiet NaN or positive infinity. 52 * Fall thru to negative test if negative infinity. 53 */ 54 if (Dbl_iszero_sign(srcp1) || 55 Dbl_isnotzero_mantissa(srcp1,srcp2)) { 56 Dbl_copytoptr(srcp1,srcp2,dstptr); 57 return(NOEXCEPTION); 58 } 59 } 60 61 /* 62 * check for zero source operand 63 */ 64 if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 65 Dbl_copytoptr(srcp1,srcp2,dstptr); 66 return(NOEXCEPTION); 67 } 68 69 /* 70 * check for negative source operand 71 */ 72 if (Dbl_isone_sign(srcp1)) { 73 /* trap if INVALIDTRAP enabled */ 74 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 75 /* make NaN quiet */ 76 Set_invalidflag(); 77 Dbl_makequietnan(srcp1,srcp2); 78 Dbl_copytoptr(srcp1,srcp2,dstptr); 79 return(NOEXCEPTION); 80 } 81 82 /* 83 * Generate result 84 */ 85 if (src_exponent > 0) { 86 even_exponent = Dbl_hidden(srcp1); 87 Dbl_clear_signexponent_set_hidden(srcp1); 88 } 89 else { 90 /* normalize operand */ 91 Dbl_clear_signexponent(srcp1); 92 src_exponent++; 93 Dbl_normalize(srcp1,srcp2,src_exponent); 94 even_exponent = src_exponent & 1; 95 } 96 if (even_exponent) { 97 /* exponent is even */ 98 /* Add comment here. Explain why odd exponent needs correction */ 99 Dbl_leftshiftby1(srcp1,srcp2); 100 } 101 /* 102 * Add comment here. Explain following algorithm. 103 * 104 * Trust me, it works. 105 * 106 */ 107 Dbl_setzero(resultp1,resultp2); 108 Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 109 Dbl_setzero_mantissap2(newbitp2); 110 while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 111 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 112 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 113 Dbl_leftshiftby1(newbitp1,newbitp2); 114 /* update result */ 115 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 116 resultp1,resultp2); 117 Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 118 Dbl_rightshiftby2(newbitp1,newbitp2); 119 } 120 else { 121 Dbl_rightshiftby1(newbitp1,newbitp2); 122 } 123 Dbl_leftshiftby1(srcp1,srcp2); 124 } 125 /* correct exponent for pre-shift */ 126 if (even_exponent) { 127 Dbl_rightshiftby1(resultp1,resultp2); 128 } 129 130 /* check for inexact */ 131 if (Dbl_isnotzero(srcp1,srcp2)) { 132 if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 133 Dbl_increment(resultp1,resultp2); 134 } 135 guardbit = Dbl_lowmantissap2(resultp2); 136 Dbl_rightshiftby1(resultp1,resultp2); 137 138 /* now round result */ 139 switch (Rounding_mode()) { 140 case ROUNDPLUS: 141 Dbl_increment(resultp1,resultp2); 142 break; 143 case ROUNDNEAREST: 144 /* stickybit is always true, so guardbit 145 * is enough to determine rounding */ 146 if (guardbit) { 147 Dbl_increment(resultp1,resultp2); 148 } 149 break; 150 } 151 /* increment result exponent by 1 if mantissa overflowed */ 152 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 153 154 if (Is_inexacttrap_enabled()) { 155 Dbl_set_exponent(resultp1, 156 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 157 Dbl_copytoptr(resultp1,resultp2,dstptr); 158 return(INEXACTEXCEPTION); 159 } 160 else Set_inexactflag(); 161 } 162 else { 163 Dbl_rightshiftby1(resultp1,resultp2); 164 } 165 Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 166 Dbl_copytoptr(resultp1,resultp2,dstptr); 167 return(NOEXCEPTION); 168 } 169