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
13const 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
19const int32_t kFixedInvGain1 = 0x18bde0bb;  // 0.607252935
20
21static 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
47SkFixed 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
68SkFixed SkCordicTan(SkFixed a)
69{
70    int32_t cos;
71    int32_t sin = SkCordicSinCos(a, &cos);
72    return SkFixedDiv(sin, cos);
73}
74
75static 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
101SkFixed 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
114SkFixed SkCordicACos(SkFixed a) {
115    int32_t z = SkCordicASin(a);
116    z = (SK_FixedPI>>1) - z;
117    return z;
118}
119
120SkFixed 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
136const 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
143const int32_t kFixedInvGain2 = 0x31330AAA;  // 1.207534495
144
145static 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
177SkFixed 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
189SkFixed 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
200void 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