Home | History | Annotate | Download | only in core 
      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 
     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 
     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 
     78 SkFixed SkCordicTan(SkFixed a)
     79 {
     80     int32_t cos;
     81     int32_t sin = SkCordicSinCos(a, &cos);
     82     return SkFixedDiv(sin, cos);
     83 }
     84 
     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 
    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 
    124 SkFixed SkCordicACos(SkFixed a) {
    125     int32_t z = SkCordicASin(a);
    126     z = (SK_FixedPI>>1) - z;
    127     return z;
    128 }
    129 
    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 
    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 
    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 
    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 
    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