Changeset ca113cf in mainline for uspace/lib/math/generic/sqrt.c
- Timestamp:
- 2021-06-18T16:06:10Z (3 years ago)
- Children:
- f23dbf4
- Parents:
- 197d59d
- git-author:
- Maurizio Lombardi <mlombard@…> (2018-11-25 16:51:03)
- git-committer:
- Maurizio Lombardi <mlombard@…> (2021-06-18 16:06:10)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/math/generic/sqrt.c
r197d59d rca113cf 92 92 #include "internal.h" 93 93 94 static const double tiny =1.0e-300;94 static const double one = 1.0, tiny=1.0e-300; 95 95 96 96 double sqrt(double x) … … 101 101 uint32_t r,t1,s1,ix1,q1; 102 102 103 EXTRACT_WORDS(ix0, ix1, x); 104 105 /* take care of Inf and NaN */ 106 if ((ix0&0x7ff00000) == 0x7ff00000) { 107 return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ 108 } 109 /* take care of zero */ 110 if (ix0 <= 0) { 111 if (((ix0&~sign)|ix1) == 0) 112 return x; /* sqrt(+-0) = +-0 */ 113 if (ix0 < 0) 114 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 115 } 116 /* normalize x */ 117 m = ix0>>20; 118 if (m == 0) { /* subnormal x */ 119 while (ix0 == 0) { 120 m -= 21; 121 ix0 |= (ix1>>11); 122 ix1 <<= 21; 123 } 124 for (i=0; (ix0&0x00100000) == 0; i++) 125 ix0<<=1; 126 m -= i - 1; 127 ix0 |= ix1>>(32-i); 128 ix1 <<= i; 129 } 130 m -= 1023; /* unbias exponent */ 103 EXTRACT_WORDS(ix0,ix1,x); 104 105 /* take care of Inf and NaN */ 106 if((ix0&0x7ff00000)==0x7ff00000) { 107 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 108 sqrt(-inf)=sNaN */ 109 } 110 /* take care of zero */ 111 if(ix0<=0) { 112 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ 113 else if(ix0<0) 114 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 115 } 116 /* normalize x */ 117 m = (ix0>>20); 118 if(m==0) { /* subnormal x */ 119 while(ix0==0) { 120 m -= 21; 121 ix0 |= (ix1>>11); ix1 <<= 21; 122 } 123 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; 124 m -= i-1; 125 ix0 |= (ix1>>(32-i)); 126 ix1 <<= i; 127 } 128 m -= 1023; /* unbias exponent */ 131 129 ix0 = (ix0&0x000fffff)|0x00100000; 132 if (m & 1) {/* odd m, double x to make it even */133 134 135 } 136 m >>= 1; 137 138 130 if(m&1){ /* odd m, double x to make it even */ 131 ix0 += ix0 + ((ix1&sign)>>31); 132 ix1 += ix1; 133 } 134 m >>= 1; /* m = [m/2] */ 135 136 /* generate sqrt(x) bit by bit */ 139 137 ix0 += ix0 + ((ix1&sign)>>31); 140 138 ix1 += ix1; 141 q = q1 = s0 = s1 = 0; 142 r = 0x00200000; 143 144 while (r !=0) {145 t = s0 + r;146 if (t <= ix0) {147 s0 = t + r;148 ix0 -= t;149 q += r;150 }151 152 153 r >>=1;139 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ 140 r = 0x00200000; /* r = moving bit from right to left */ 141 142 while(r!=0) { 143 t = s0+r; 144 if(t<=ix0) { 145 s0 = t+r; 146 ix0 -= t; 147 q += r; 148 } 149 ix0 += ix0 + ((ix1&sign)>>31); 150 ix1 += ix1; 151 r>>=1; 154 152 } 155 153 156 154 r = sign; 157 while (r != 0) { 158 t1 = s1 + r; 159 t = s0; 160 if (t < ix0 || (t == ix0 && t1 <= ix1)) { 161 s1 = t1 + r; 162 if ((t1&sign) == (uint32_t)sign && (s1&sign) == 0) 163 s0++; 164 ix0 -= t; 165 if (ix1 < t1) 166 ix0--; 167 ix1 -= t1; 168 q1 += r; 169 } 170 ix0 += ix0 + ((ix1&sign)>>31); 171 ix1 += ix1; 172 r >>= 1; 173 } 174 175 /* use floating add to find out rounding direction */ 176 if ((ix0|ix1) != 0) { 177 z = 1.0 - tiny; /* raise inexact flag */ 178 if (z >= 1.0) { 179 z = 1.0 + tiny; 180 if (q1 == (uint32_t)0xffffffff) { 181 q1 = 0; 182 q++; 183 } else if (z > 1.0) { 184 if (q1 == (uint32_t)0xfffffffe) 185 q++; 186 q1 += 2; 187 } else 188 q1 += q1 & 1; 189 } 190 } 191 ix0 = (q>>1) + 0x3fe00000; 192 ix1 = q1>>1; 193 if (q&1) 194 ix1 |= sign; 195 ix0 += m << 20; 196 INSERT_WORDS(z, ix0, ix1); 155 while(r!=0) { 156 t1 = s1+r; 157 t = s0; 158 if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 159 s1 = t1+r; 160 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1; 161 ix0 -= t; 162 if (ix1 < t1) ix0 -= 1; 163 ix1 -= t1; 164 q1 += r; 165 } 166 ix0 += ix0 + ((ix1&sign)>>31); 167 ix1 += ix1; 168 r>>=1; 169 } 170 171 /* use floating add to find out rounding direction */ 172 if((ix0|ix1)!=0) { 173 z = one-tiny; /* trigger inexact flag */ 174 if (z>=one) { 175 z = one+tiny; 176 if (q1==(uint32_t)0xffffffff) { q1=0; q += 1;} 177 else if (z>one) { 178 if (q1==(uint32_t)0xfffffffe) q+=1; 179 q1+=2; 180 } else 181 q1 += (q1&1); 182 } 183 } 184 ix0 = (q>>1)+0x3fe00000; 185 ix1 = q1>>1; 186 if ((q&1)==1) ix1 |= sign; 187 ix0 += (m <<20); 188 INSERT_WORDS(z,ix0,ix1); 197 189 return z; 198 190 }
Note:
See TracChangeset
for help on using the changeset viewer.