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
15const 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
21const int32_t kFixedInvGain1 = 0x18bde0bb;  // 0.607252935
22
23static 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
49SkFixed 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
70SkFixed SkCordicTan(SkFixed a)
71{
72    int32_t cos;
73    int32_t sin = SkCordicSinCos(a, &cos);
74    return SkFixedDiv(sin, cos);
75}
76
77static 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
103SkFixed 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
116SkFixed SkCordicACos(SkFixed a) {
117    int32_t z = SkCordicASin(a);
118    z = (SK_FixedPI>>1) - z;
119    return z;
120}
121
122SkFixed 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
138const 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
145const int32_t kFixedInvGain2 = 0x31330AAA;  // 1.207534495
146
147static 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
179SkFixed 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
191SkFixed 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
204void 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