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