Chromium Code Reviews
chromiumcodereview-hr@appspot.gserviceaccount.com (chromiumcodereview-hr) | Please choose your nickname with Settings | Help | Chromium Project | Gerrit Changes | Sign out
(132)

Side by Side Diff: fusl/src/math/sqrt.c

Issue 1714623002: [fusl] clang-format fusl (Closed) Base URL: git@github.com:domokit/mojo.git@master
Patch Set: headers too Created 4 years, 10 months ago
Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.
Jump to:
View unified diff | Download patch
OLDNEW
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
2 /* 2 /*
3 * ==================================================== 3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * 5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this 7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice 8 * software is freely granted, provided that this notice
9 * is preserved. 9 * is preserved.
10 * ==================================================== 10 * ====================================================
(...skipping 62 matching lines...) Expand 10 before | Expand all | Expand 10 after
73 * sqrt(+-0) = +-0 ... exact 73 * sqrt(+-0) = +-0 ... exact
74 * sqrt(inf) = inf 74 * sqrt(inf) = inf
75 * sqrt(-ve) = NaN ... with invalid signal 75 * sqrt(-ve) = NaN ... with invalid signal
76 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN 76 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
77 */ 77 */
78 78
79 #include "libm.h" 79 #include "libm.h"
80 80
81 static const double tiny = 1.0e-300; 81 static const double tiny = 1.0e-300;
82 82
83 double sqrt(double x) 83 double sqrt(double x) {
84 { 84 double z;
85 » double z; 85 int32_t sign = (int)0x80000000;
86 » int32_t sign = (int)0x80000000; 86 int32_t ix0, s0, q, m, t, i;
87 » int32_t ix0,s0,q,m,t,i; 87 uint32_t r, t1, s1, ix1, q1;
88 » uint32_t r,t1,s1,ix1,q1;
89 88
90 » EXTRACT_WORDS(ix0, ix1, x); 89 EXTRACT_WORDS(ix0, ix1, x);
91 90
92 » /* take care of Inf and NaN */ 91 /* take care of Inf and NaN */
93 » if ((ix0&0x7ff00000) == 0x7ff00000) { 92 if ((ix0 & 0x7ff00000) == 0x7ff00000) {
94 » » return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=s NaN */ 93 return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
95 » } 94 }
96 » /* take care of zero */ 95 /* take care of zero */
97 » if (ix0 <= 0) { 96 if (ix0 <= 0) {
98 » » if (((ix0&~sign)|ix1) == 0) 97 if (((ix0 & ~sign) | ix1) == 0)
99 » » » return x; /* sqrt(+-0) = +-0 */ 98 return x; /* sqrt(+-0) = +-0 */
100 » » if (ix0 < 0) 99 if (ix0 < 0)
101 » » » return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 100 return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
102 » } 101 }
103 » /* normalize x */ 102 /* normalize x */
104 » m = ix0>>20; 103 m = ix0 >> 20;
105 » if (m == 0) { /* subnormal x */ 104 if (m == 0) { /* subnormal x */
106 » » while (ix0 == 0) { 105 while (ix0 == 0) {
107 » » » m -= 21; 106 m -= 21;
108 » » » ix0 |= (ix1>>11); 107 ix0 |= (ix1 >> 11);
109 » » » ix1 <<= 21; 108 ix1 <<= 21;
110 » » } 109 }
111 » » for (i=0; (ix0&0x00100000) == 0; i++) 110 for (i = 0; (ix0 & 0x00100000) == 0; i++)
112 » » » ix0<<=1; 111 ix0 <<= 1;
113 » » m -= i - 1; 112 m -= i - 1;
114 » » ix0 |= ix1>>(32-i); 113 ix0 |= ix1 >> (32 - i);
115 » » ix1 <<= i; 114 ix1 <<= i;
116 » } 115 }
117 » m -= 1023; /* unbias exponent */ 116 m -= 1023; /* unbias exponent */
118 » ix0 = (ix0&0x000fffff)|0x00100000; 117 ix0 = (ix0 & 0x000fffff) | 0x00100000;
119 » if (m & 1) { /* odd m, double x to make it even */ 118 if (m & 1) { /* odd m, double x to make it even */
120 » » ix0 += ix0 + ((ix1&sign)>>31); 119 ix0 += ix0 + ((ix1 & sign) >> 31);
121 » » ix1 += ix1; 120 ix1 += ix1;
122 » } 121 }
123 » m >>= 1; /* m = [m/2] */ 122 m >>= 1; /* m = [m/2] */
124 123
125 » /* generate sqrt(x) bit by bit */ 124 /* generate sqrt(x) bit by bit */
126 » ix0 += ix0 + ((ix1&sign)>>31); 125 ix0 += ix0 + ((ix1 & sign) >> 31);
127 » ix1 += ix1; 126 ix1 += ix1;
128 » q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ 127 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
129 » r = 0x00200000; /* r = moving bit from right to left */ 128 r = 0x00200000; /* r = moving bit from right to left */
130 129
131 » while (r != 0) { 130 while (r != 0) {
132 » » t = s0 + r; 131 t = s0 + r;
133 » » if (t <= ix0) { 132 if (t <= ix0) {
134 » » » s0 = t + r; 133 s0 = t + r;
135 » » » ix0 -= t; 134 ix0 -= t;
136 » » » q += r; 135 q += r;
137 » » } 136 }
138 » » ix0 += ix0 + ((ix1&sign)>>31); 137 ix0 += ix0 + ((ix1 & sign) >> 31);
139 » » ix1 += ix1; 138 ix1 += ix1;
140 » » r >>= 1; 139 r >>= 1;
141 » } 140 }
142 141
143 » r = sign; 142 r = sign;
144 » while (r != 0) { 143 while (r != 0) {
145 » » t1 = s1 + r; 144 t1 = s1 + r;
146 » » t = s0; 145 t = s0;
147 » » if (t < ix0 || (t == ix0 && t1 <= ix1)) { 146 if (t < ix0 || (t == ix0 && t1 <= ix1)) {
148 » » » s1 = t1 + r; 147 s1 = t1 + r;
149 » » » if ((t1&sign) == sign && (s1&sign) == 0) 148 if ((t1 & sign) == sign && (s1 & sign) == 0)
150 » » » » s0++; 149 s0++;
151 » » » ix0 -= t; 150 ix0 -= t;
152 » » » if (ix1 < t1) 151 if (ix1 < t1)
153 » » » » ix0--; 152 ix0--;
154 » » » ix1 -= t1; 153 ix1 -= t1;
155 » » » q1 += r; 154 q1 += r;
156 » » } 155 }
157 » » ix0 += ix0 + ((ix1&sign)>>31); 156 ix0 += ix0 + ((ix1 & sign) >> 31);
158 » » ix1 += ix1; 157 ix1 += ix1;
159 » » r >>= 1; 158 r >>= 1;
160 » } 159 }
161 160
162 » /* use floating add to find out rounding direction */ 161 /* use floating add to find out rounding direction */
163 » if ((ix0|ix1) != 0) { 162 if ((ix0 | ix1) != 0) {
164 » » z = 1.0 - tiny; /* raise inexact flag */ 163 z = 1.0 - tiny; /* raise inexact flag */
165 » » if (z >= 1.0) { 164 if (z >= 1.0) {
166 » » » z = 1.0 + tiny; 165 z = 1.0 + tiny;
167 » » » if (q1 == (uint32_t)0xffffffff) { 166 if (q1 == (uint32_t)0xffffffff) {
168 » » » » q1 = 0; 167 q1 = 0;
169 » » » » q++; 168 q++;
170 » » » } else if (z > 1.0) { 169 } else if (z > 1.0) {
171 » » » » if (q1 == (uint32_t)0xfffffffe) 170 if (q1 == (uint32_t)0xfffffffe)
172 » » » » » q++; 171 q++;
173 » » » » q1 += 2; 172 q1 += 2;
174 » » » } else 173 } else
175 » » » » q1 += q1 & 1; 174 q1 += q1 & 1;
176 » » } 175 }
177 » } 176 }
178 » ix0 = (q>>1) + 0x3fe00000; 177 ix0 = (q >> 1) + 0x3fe00000;
179 » ix1 = q1>>1; 178 ix1 = q1 >> 1;
180 » if (q&1) 179 if (q & 1)
181 » » ix1 |= sign; 180 ix1 |= sign;
182 » ix0 += m << 20; 181 ix0 += m << 20;
183 » INSERT_WORDS(z, ix0, ix1); 182 INSERT_WORDS(z, ix0, ix1);
184 » return z; 183 return z;
185 } 184 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698