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