• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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