1 // SPDX-License-Identifier: GPL-2.0-or-later 2 /* 3 * Linux/PA-RISC Project (http://www.parisc-linux.org/) 4 * 5 * Floating-point emulation code 6 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> 7 */ 8 /* 9 * BEGIN_DESC 10 * 11 * File: 12 * @(#) pa/spmath/dfsqrt.c $Revision: 1.1 $ 13 * 14 * Purpose: 15 * Double Floating-point Square Root 16 * 17 * External Interfaces: 18 * dbl_fsqrt(srcptr,_nullptr,dstptr,status) 19 * 20 * Internal Interfaces: 21 * 22 * Theory: 23 * <<please update with a overview of the operation of this file>> 24 * 25 * END_DESC 26 */ 27 28 29 #include "float.h" 30 #include "dbl_float.h" 31 32 /* 33 * Double Floating-point Square Root 34 */ 35 36 /*ARGSUSED*/ 37 unsigned int 38 dbl_fsqrt( 39 dbl_floating_point *srcptr, 40 unsigned int *_nullptr, 41 dbl_floating_point *dstptr, 42 unsigned int *status) 43 { 44 register unsigned int srcp1, srcp2, resultp1, resultp2; 45 register unsigned int newbitp1, newbitp2, sump1, sump2; 46 register int src_exponent; 47 register boolean guardbit = FALSE, even_exponent; 48 49 Dbl_copyfromptr(srcptr,srcp1,srcp2); 50 /* 51 * check source operand for NaN or infinity 52 */ 53 if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 54 /* 55 * is signaling NaN? 56 */ 57 if (Dbl_isone_signaling(srcp1)) { 58 /* trap if INVALIDTRAP enabled */ 59 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 60 /* make NaN quiet */ 61 Set_invalidflag(); 62 Dbl_set_quiet(srcp1); 63 } 64 /* 65 * Return quiet NaN or positive infinity. 66 * Fall through to negative test if negative infinity. 67 */ 68 if (Dbl_iszero_sign(srcp1) || 69 Dbl_isnotzero_mantissa(srcp1,srcp2)) { 70 Dbl_copytoptr(srcp1,srcp2,dstptr); 71 return(NOEXCEPTION); 72 } 73 } 74 75 /* 76 * check for zero source operand 77 */ 78 if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 79 Dbl_copytoptr(srcp1,srcp2,dstptr); 80 return(NOEXCEPTION); 81 } 82 83 /* 84 * check for negative source operand 85 */ 86 if (Dbl_isone_sign(srcp1)) { 87 /* trap if INVALIDTRAP enabled */ 88 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 89 /* make NaN quiet */ 90 Set_invalidflag(); 91 Dbl_makequietnan(srcp1,srcp2); 92 Dbl_copytoptr(srcp1,srcp2,dstptr); 93 return(NOEXCEPTION); 94 } 95 96 /* 97 * Generate result 98 */ 99 if (src_exponent > 0) { 100 even_exponent = Dbl_hidden(srcp1); 101 Dbl_clear_signexponent_set_hidden(srcp1); 102 } 103 else { 104 /* normalize operand */ 105 Dbl_clear_signexponent(srcp1); 106 src_exponent++; 107 Dbl_normalize(srcp1,srcp2,src_exponent); 108 even_exponent = src_exponent & 1; 109 } 110 if (even_exponent) { 111 /* exponent is even */ 112 /* Add comment here. Explain why odd exponent needs correction */ 113 Dbl_leftshiftby1(srcp1,srcp2); 114 } 115 /* 116 * Add comment here. Explain following algorithm. 117 * 118 * Trust me, it works. 119 * 120 */ 121 Dbl_setzero(resultp1,resultp2); 122 Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 123 Dbl_setzero_mantissap2(newbitp2); 124 while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 125 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 126 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 127 Dbl_leftshiftby1(newbitp1,newbitp2); 128 /* update result */ 129 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 130 resultp1,resultp2); 131 Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 132 Dbl_rightshiftby2(newbitp1,newbitp2); 133 } 134 else { 135 Dbl_rightshiftby1(newbitp1,newbitp2); 136 } 137 Dbl_leftshiftby1(srcp1,srcp2); 138 } 139 /* correct exponent for pre-shift */ 140 if (even_exponent) { 141 Dbl_rightshiftby1(resultp1,resultp2); 142 } 143 144 /* check for inexact */ 145 if (Dbl_isnotzero(srcp1,srcp2)) { 146 if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 147 Dbl_increment(resultp1,resultp2); 148 } 149 guardbit = Dbl_lowmantissap2(resultp2); 150 Dbl_rightshiftby1(resultp1,resultp2); 151 152 /* now round result */ 153 switch (Rounding_mode()) { 154 case ROUNDPLUS: 155 Dbl_increment(resultp1,resultp2); 156 break; 157 case ROUNDNEAREST: 158 /* stickybit is always true, so guardbit 159 * is enough to determine rounding */ 160 if (guardbit) { 161 Dbl_increment(resultp1,resultp2); 162 } 163 break; 164 } 165 /* increment result exponent by 1 if mantissa overflowed */ 166 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 167 168 if (Is_inexacttrap_enabled()) { 169 Dbl_set_exponent(resultp1, 170 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 171 Dbl_copytoptr(resultp1,resultp2,dstptr); 172 return(INEXACTEXCEPTION); 173 } 174 else Set_inexactflag(); 175 } 176 else { 177 Dbl_rightshiftby1(resultp1,resultp2); 178 } 179 Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 180 Dbl_copytoptr(resultp1,resultp2,dstptr); 181 return(NOEXCEPTION); 182 } 183