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