1 /* libs/corecg/SkCordic.cpp
2 **
3 ** Copyright 2006, The Android Open Source Project
4 **
5 ** Licensed under the Apache License, Version 2.0 (the "License");
6 ** you may not use this file except in compliance with the License.
7 ** You may obtain a copy of the License at
8 **
9 ** http://www.apache.org/licenses/LICENSE-2.0
10 **
11 ** Unless required by applicable law or agreed to in writing, software
12 ** distributed under the License is distributed on an "AS IS" BASIS,
13 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 ** See the License for the specific language governing permissions and
15 ** limitations under the License.
16 */
17
18 #include "SkCordic.h"
19 #include "SkMath.h"
20 #include "Sk64.h"
21
22 // 0x20000000 equals pi / 4
23 const int32_t kATanDegrees[] = { 0x20000000,
24 0x12E4051D, 0x9FB385B, 0x51111D4, 0x28B0D43, 0x145D7E1, 0xA2F61E, 0x517C55,
25 0x28BE53, 0x145F2E, 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
26 0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
27 0xA, 0x5, 0x2, 0x1 };
28
29 const int32_t kFixedInvGain1 = 0x18bde0bb; // 0.607252935
30
SkCircularRotation(int32_t * x0,int32_t * y0,int32_t * z0)31 static void SkCircularRotation(int32_t* x0, int32_t* y0, int32_t* z0)
32 {
33 int32_t t = 0;
34 int32_t x = *x0;
35 int32_t y = *y0;
36 int32_t z = *z0;
37 const int32_t* tanPtr = kATanDegrees;
38 do {
39 int32_t x1 = y >> t;
40 int32_t y1 = x >> t;
41 int32_t tan = *tanPtr++;
42 if (z >= 0) {
43 x -= x1;
44 y += y1;
45 z -= tan;
46 } else {
47 x += x1;
48 y -= y1;
49 z += tan;
50 }
51 } while (++t < 16); // 30);
52 *x0 = x;
53 *y0 = y;
54 *z0 = z;
55 }
56
SkCordicSinCos(SkFixed radians,SkFixed * cosp)57 SkFixed SkCordicSinCos(SkFixed radians, SkFixed* cosp)
58 {
59 int32_t scaledRadians = radians * 0x28be; // scale radians to 65536 / PI()
60 int quadrant = scaledRadians >> 30;
61 quadrant += 1;
62 if (quadrant & 2)
63 scaledRadians = -scaledRadians + 0x80000000;
64 /* |a| <= 90 degrees as a 1.31 number */
65 SkFixed sin = 0;
66 SkFixed cos = kFixedInvGain1;
67 SkCircularRotation(&cos, &sin, &scaledRadians);
68 Sk64 scaled;
69 scaled.setMul(sin, 0x6488d);
70 sin = scaled.fHi;
71 scaled.setMul(cos, 0x6488d);
72 if (quadrant & 2)
73 scaled.fHi = - scaled.fHi;
74 *cosp = scaled.fHi;
75 return sin;
76 }
77
SkCordicTan(SkFixed a)78 SkFixed SkCordicTan(SkFixed a)
79 {
80 int32_t cos;
81 int32_t sin = SkCordicSinCos(a, &cos);
82 return SkFixedDiv(sin, cos);
83 }
84
SkCircularVector(int32_t * y0,int32_t * x0,int32_t vecMode)85 static int32_t SkCircularVector(int32_t* y0, int32_t* x0, int32_t vecMode)
86 {
87 int32_t x = *x0;
88 int32_t y = *y0;
89 int32_t z = 0;
90 int32_t t = 0;
91 const int32_t* tanPtr = kATanDegrees;
92 do {
93 int32_t x1 = y >> t;
94 int32_t y1 = x >> t;
95 int32_t tan = *tanPtr++;
96 if (y < vecMode) {
97 x -= x1;
98 y += y1;
99 z -= tan;
100 } else {
101 x += x1;
102 y -= y1;
103 z += tan;
104 }
105 } while (++t < 16); // 30
106 Sk64 scaled;
107 scaled.setMul(z, 0x6488d); // scale back into the SkScalar space (0x100000000/0x28be)
108 return scaled.fHi;
109 }
110
SkCordicASin(SkFixed a)111 SkFixed SkCordicASin(SkFixed a) {
112 int32_t sign = SkExtractSign(a);
113 int32_t z = SkFixedAbs(a);
114 if (z >= SK_Fixed1)
115 return SkApplySign(SK_FixedPI>>1, sign);
116 int32_t x = kFixedInvGain1;
117 int32_t y = 0;
118 z *= 0x28be;
119 z = SkCircularVector(&y, &x, z);
120 z = SkApplySign(z, ~sign);
121 return z;
122 }
123
SkCordicACos(SkFixed a)124 SkFixed SkCordicACos(SkFixed a) {
125 int32_t z = SkCordicASin(a);
126 z = (SK_FixedPI>>1) - z;
127 return z;
128 }
129
SkCordicATan2(SkFixed y,SkFixed x)130 SkFixed SkCordicATan2(SkFixed y, SkFixed x) {
131 if ((x | y) == 0)
132 return 0;
133 int32_t xsign = SkExtractSign(x);
134 x = SkFixedAbs(x);
135 int32_t result = SkCircularVector(&y, &x, 0);
136 if (xsign) {
137 int32_t rsign = SkExtractSign(result);
138 if (y == 0)
139 rsign = 0;
140 SkFixed pi = SkApplySign(SK_FixedPI, rsign);
141 result = pi - result;
142 }
143 return result;
144 }
145
146 const int32_t kATanHDegrees[] = {
147 0x1661788D, 0xA680D61, 0x51EA6FC, 0x28CBFDD, 0x1460E34,
148 0xA2FCE8, 0x517D2E, 0x28BE6E, 0x145F32,
149 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
150 0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
151 0xA, 0x5, 0x2, 0x1 };
152
153 const int32_t kFixedInvGain2 = 0x31330AAA; // 1.207534495
154
SkHyperbolic(int32_t * x0,int32_t * y0,int32_t * z0,int mode)155 static void SkHyperbolic(int32_t* x0, int32_t* y0, int32_t* z0, int mode)
156 {
157 int32_t t = 1;
158 int32_t x = *x0;
159 int32_t y = *y0;
160 int32_t z = *z0;
161 const int32_t* tanPtr = kATanHDegrees;
162 int k = -3;
163 do {
164 int32_t x1 = y >> t;
165 int32_t y1 = x >> t;
166 int32_t tan = *tanPtr++;
167 int count = 2 + (k >> 31);
168 if (++k == 1)
169 k = -2;
170 do {
171 if (((y >> 31) & mode) | ~((z >> 31) | mode)) {
172 x += x1;
173 y += y1;
174 z -= tan;
175 } else {
176 x -= x1;
177 y -= y1;
178 z += tan;
179 }
180 } while (--count);
181 } while (++t < 30);
182 *x0 = x;
183 *y0 = y;
184 *z0 = z;
185 }
186
SkCordicLog(SkFixed a)187 SkFixed SkCordicLog(SkFixed a) {
188 a *= 0x28be;
189 int32_t x = a + 0x28BE60DB; // 1.0
190 int32_t y = a - 0x28BE60DB;
191 int32_t z = 0;
192 SkHyperbolic(&x, &y, &z, -1);
193 Sk64 scaled;
194 scaled.setMul(z, 0x6488d);
195 z = scaled.fHi;
196 return z << 1;
197 }
198
SkCordicExp(SkFixed a)199 SkFixed SkCordicExp(SkFixed a) {
200 int32_t cosh = kFixedInvGain2;
201 int32_t sinh = 0;
202 SkHyperbolic(&cosh, &sinh, &a, 0);
203 return cosh + sinh;
204 }
205
206 #ifdef SK_DEBUG
207
208 #ifdef SK_CAN_USE_FLOAT
209 #include "SkFloatingPoint.h"
210 #endif
211
SkCordic_UnitTest()212 void SkCordic_UnitTest()
213 {
214 #if defined(SK_SUPPORT_UNITTEST) && defined(SK_CAN_USE_FLOAT)
215 float val;
216 for (float angle = -720; angle < 720; angle += 30) {
217 float radian = angle * 3.1415925358f / 180.0f;
218 SkFixed f_angle = (int) (radian * 65536.0f);
219 // sincos
220 float sine = sinf(radian);
221 float cosine = cosf(radian);
222 SkFixed f_cosine;
223 SkFixed f_sine = SkCordicSinCos(f_angle, &f_cosine);
224 float sine2 = (float) f_sine / 65536.0f;
225 float cosine2 = (float) f_cosine / 65536.0f;
226 float error = fabsf(sine - sine2);
227 if (error > 0.001)
228 SkDebugf("sin error : angle = %g ; sin = %g ; cordic = %g\n", angle, sine, sine2);
229 error = fabsf(cosine - cosine2);
230 if (error > 0.001)
231 SkDebugf("cos error : angle = %g ; cos = %g ; cordic = %g\n", angle, cosine, cosine2);
232 // tan
233 float _tan = tanf(radian);
234 SkFixed f_tan = SkCordicTan(f_angle);
235 float tan2 = (float) f_tan / 65536.0f;
236 error = fabsf(_tan - tan2);
237 if (error > 0.05 && fabsf(_tan) < 1e6)
238 SkDebugf("tan error : angle = %g ; tan = %g ; cordic = %g\n", angle, _tan, tan2);
239 }
240 for (val = -1; val <= 1; val += .1f) {
241 SkFixed f_val = (int) (val * 65536.0f);
242 // asin
243 float arcsine = asinf(val);
244 SkFixed f_arcsine = SkCordicASin(f_val);
245 float arcsine2 = (float) f_arcsine / 65536.0f;
246 float error = fabsf(arcsine - arcsine2);
247 if (error > 0.001)
248 SkDebugf("asin error : val = %g ; asin = %g ; cordic = %g\n", val, arcsine, arcsine2);
249 }
250 #if 1
251 for (val = -1; val <= 1; val += .1f) {
252 #else
253 val = .5; {
254 #endif
255 SkFixed f_val = (int) (val * 65536.0f);
256 // acos
257 float arccos = acosf(val);
258 SkFixed f_arccos = SkCordicACos(f_val);
259 float arccos2 = (float) f_arccos / 65536.0f;
260 float error = fabsf(arccos - arccos2);
261 if (error > 0.001)
262 SkDebugf("acos error : val = %g ; acos = %g ; cordic = %g\n", val, arccos, arccos2);
263 }
264 // atan2
265 #if 1
266 for (val = -1000; val <= 1000; val += 500.f) {
267 for (float val2 = -1000; val2 <= 1000; val2 += 500.f) {
268 #else
269 val = 0; {
270 float val2 = -1000; {
271 #endif
272 SkFixed f_val = (int) (val * 65536.0f);
273 SkFixed f_val2 = (int) (val2 * 65536.0f);
274 float arctan = atan2f(val, val2);
275 SkFixed f_arctan = SkCordicATan2(f_val, f_val2);
276 float arctan2 = (float) f_arctan / 65536.0f;
277 float error = fabsf(arctan - arctan2);
278 if (error > 0.001)
279 SkDebugf("atan2 error : val = %g ; val2 = %g ; atan2 = %g ; cordic = %g\n", val, val2, arctan, arctan2);
280 }
281 }
282 // log
283 #if 1
284 for (val = 0.125f; val <= 8.f; val *= 2.0f) {
285 #else
286 val = .5; {
287 #endif
288 SkFixed f_val = (int) (val * 65536.0f);
289 // acos
290 float log = logf(val);
291 SkFixed f_log = SkCordicLog(f_val);
292 float log2 = (float) f_log / 65536.0f;
293 float error = fabsf(log - log2);
294 if (error > 0.001)
295 SkDebugf("log error : val = %g ; log = %g ; cordic = %g\n", val, log, log2);
296 }
297 // exp
298 #endif
299 }
300
301 #endif
302