Home | History | Annotate | Download | only in core
      1 /*
      2  * Copyright (C) 2006-2008 The Android Open Source Project
      3  *
      4  * Licensed under the Apache License, Version 2.0 (the "License");
      5  * you may not use this file except in compliance with the License.
      6  * You may obtain a copy of the License at
      7  *
      8  *      http://www.apache.org/licenses/LICENSE-2.0
      9  *
     10  * Unless required by applicable law or agreed to in writing, software
     11  * distributed under the License is distributed on an "AS IS" BASIS,
     12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     13  * See the License for the specific language governing permissions and
     14  * limitations under the License.
     15  */
     16 
     17 #include "SkMath.h"
     18 #include "SkCordic.h"
     19 #include "SkFloatBits.h"
     20 #include "SkFloatingPoint.h"
     21 #include "Sk64.h"
     22 #include "SkScalar.h"
     23 
     24 #ifdef SK_SCALAR_IS_FLOAT
     25     const uint32_t gIEEENotANumber = 0x7FFFFFFF;
     26     const uint32_t gIEEEInfinity = 0x7F800000;
     27 #endif
     28 
     29 #define sub_shift(zeros, x, n)  \
     30     zeros -= n;                 \
     31     x >>= n
     32 
     33 int SkCLZ_portable(uint32_t x) {
     34     if (x == 0) {
     35         return 32;
     36     }
     37 
     38 #ifdef SK_CPU_HAS_CONDITIONAL_INSTR
     39     int zeros = 31;
     40     if (x & 0xFFFF0000) {
     41         sub_shift(zeros, x, 16);
     42     }
     43     if (x & 0xFF00) {
     44         sub_shift(zeros, x, 8);
     45     }
     46     if (x & 0xF0) {
     47         sub_shift(zeros, x, 4);
     48     }
     49     if (x & 0xC) {
     50         sub_shift(zeros, x, 2);
     51     }
     52     if (x & 0x2) {
     53         sub_shift(zeros, x, 1);
     54     }
     55 #else
     56     int zeros = ((x >> 16) - 1) >> 31 << 4;
     57     x <<= zeros;
     58 
     59     int nonzero = ((x >> 24) - 1) >> 31 << 3;
     60     zeros += nonzero;
     61     x <<= nonzero;
     62 
     63     nonzero = ((x >> 28) - 1) >> 31 << 2;
     64     zeros += nonzero;
     65     x <<= nonzero;
     66 
     67     nonzero = ((x >> 30) - 1) >> 31 << 1;
     68     zeros += nonzero;
     69     x <<= nonzero;
     70 
     71     zeros += (~x) >> 31;
     72 #endif
     73 
     74     return zeros;
     75 }
     76 
     77 int32_t SkMulDiv(int32_t numer1, int32_t numer2, int32_t denom) {
     78     SkASSERT(denom);
     79 
     80     Sk64 tmp;
     81     tmp.setMul(numer1, numer2);
     82     tmp.div(denom, Sk64::kTrunc_DivOption);
     83     return tmp.get32();
     84 }
     85 
     86 int32_t SkMulShift(int32_t a, int32_t b, unsigned shift) {
     87     int sign = SkExtractSign(a ^ b);
     88 
     89     if (shift > 63) {
     90         return sign;
     91     }
     92 
     93     a = SkAbs32(a);
     94     b = SkAbs32(b);
     95 
     96     uint32_t ah = a >> 16;
     97     uint32_t al = a & 0xFFFF;
     98     uint32_t bh = b >> 16;
     99     uint32_t bl = b & 0xFFFF;
    100 
    101     uint32_t A = ah * bh;
    102     uint32_t B = ah * bl + al * bh;
    103     uint32_t C = al * bl;
    104 
    105     /*  [  A  ]
    106            [  B  ]
    107               [  C  ]
    108     */
    109     uint32_t lo = C + (B << 16);
    110     int32_t  hi = A + (B >> 16) + (lo < C);
    111 
    112     if (sign < 0) {
    113         hi = -hi - Sk32ToBool(lo);
    114         lo = 0 - lo;
    115     }
    116 
    117     if (shift == 0) {
    118 #ifdef SK_DEBUGx
    119         SkASSERT(((int32_t)lo >> 31) == hi);
    120 #endif
    121         return lo;
    122     } else if (shift >= 32) {
    123         return hi >> (shift - 32);
    124     } else {
    125 #ifdef SK_DEBUGx
    126         int32_t tmp = hi >> shift;
    127         SkASSERT(tmp == 0 || tmp == -1);
    128 #endif
    129         // we want (hi << (32 - shift)) | (lo >> shift) but rounded
    130         int roundBit = (lo >> (shift - 1)) & 1;
    131         return ((hi << (32 - shift)) | (lo >> shift)) + roundBit;
    132     }
    133 }
    134 
    135 SkFixed SkFixedMul_portable(SkFixed a, SkFixed b) {
    136 #if 0
    137     Sk64    tmp;
    138 
    139     tmp.setMul(a, b);
    140     tmp.shiftRight(16);
    141     return tmp.fLo;
    142 #elif defined(SkLONGLONG)
    143     return static_cast<SkFixed>((SkLONGLONG)a * b >> 16);
    144 #else
    145     int sa = SkExtractSign(a);
    146     int sb = SkExtractSign(b);
    147     // now make them positive
    148     a = SkApplySign(a, sa);
    149     b = SkApplySign(b, sb);
    150 
    151     uint32_t    ah = a >> 16;
    152     uint32_t    al = a & 0xFFFF;
    153     uint32_t bh = b >> 16;
    154     uint32_t bl = b & 0xFFFF;
    155 
    156     uint32_t R = ah * b + al * bh + (al * bl >> 16);
    157 
    158     return SkApplySign(R, sa ^ sb);
    159 #endif
    160 }
    161 
    162 SkFract SkFractMul_portable(SkFract a, SkFract b) {
    163 #if 0
    164     Sk64 tmp;
    165     tmp.setMul(a, b);
    166     return tmp.getFract();
    167 #elif defined(SkLONGLONG)
    168     return static_cast<SkFract>((SkLONGLONG)a * b >> 30);
    169 #else
    170     int sa = SkExtractSign(a);
    171     int sb = SkExtractSign(b);
    172     // now make them positive
    173     a = SkApplySign(a, sa);
    174     b = SkApplySign(b, sb);
    175 
    176     uint32_t ah = a >> 16;
    177     uint32_t al = a & 0xFFFF;
    178     uint32_t bh = b >> 16;
    179     uint32_t bl = b & 0xFFFF;
    180 
    181     uint32_t A = ah * bh;
    182     uint32_t B = ah * bl + al * bh;
    183     uint32_t C = al * bl;
    184 
    185     /*  [  A  ]
    186            [  B  ]
    187               [  C  ]
    188     */
    189     uint32_t Lo = C + (B << 16);
    190     uint32_t Hi = A + (B >>16) + (Lo < C);
    191 
    192     SkASSERT((Hi >> 29) == 0);  // else overflow
    193 
    194     int32_t R = (Hi << 2) + (Lo >> 30);
    195 
    196     return SkApplySign(R, sa ^ sb);
    197 #endif
    198 }
    199 
    200 int SkFixedMulCommon(SkFixed a, int b, int bias) {
    201     // this function only works if b is 16bits
    202     SkASSERT(b == (int16_t)b);
    203     SkASSERT(b >= 0);
    204 
    205     int sa = SkExtractSign(a);
    206     a = SkApplySign(a, sa);
    207     uint32_t ah = a >> 16;
    208     uint32_t al = a & 0xFFFF;
    209     uint32_t R = ah * b + ((al * b + bias) >> 16);
    210     return SkApplySign(R, sa);
    211 }
    212 
    213 #ifdef SK_DEBUGx
    214     #define TEST_FASTINVERT
    215 #endif
    216 
    217 SkFixed SkFixedFastInvert(SkFixed x) {
    218 /*  Adapted (stolen) from gglRecip()
    219 */
    220 
    221     if (x == SK_Fixed1) {
    222         return SK_Fixed1;
    223     }
    224 
    225     int      sign = SkExtractSign(x);
    226     uint32_t a = SkApplySign(x, sign);
    227 
    228     if (a <= 2) {
    229         return SkApplySign(SK_MaxS32, sign);
    230     }
    231 
    232 #ifdef TEST_FASTINVERT
    233     SkFixed orig = a;
    234     uint32_t slow = SkFixedDiv(SK_Fixed1, a);
    235 #endif
    236 
    237     // normalize a
    238     int lz = SkCLZ(a);
    239     a = a << lz >> 16;
    240 
    241     // compute 1/a approximation (0.5 <= a < 1.0)
    242     uint32_t r = 0x17400 - a;      // (2.90625 (~2.914) - 2*a) >> 1
    243 
    244     // Newton-Raphson iteration:
    245     // x = r*(2 - a*r) = ((r/2)*(1 - a*r/2))*4
    246     r = ( (0x10000 - ((a*r)>>16)) * r ) >> 15;
    247     r = ( (0x10000 - ((a*r)>>16)) * r ) >> (30 - lz);
    248 
    249 #ifdef TEST_FASTINVERT
    250     SkDebugf("SkFixedFastInvert(%x %g) = %x %g Slow[%x %g]\n",
    251                 orig, orig/65536.,
    252                 r, r/65536.,
    253                 slow, slow/65536.);
    254 #endif
    255 
    256     return SkApplySign(r, sign);
    257 }
    258 
    259 ///////////////////////////////////////////////////////////////////////////////
    260 
    261 #define DIVBITS_ITER(n)                                 \
    262     case n:                                             \
    263         if ((numer = (numer << 1) - denom) >= 0)        \
    264             result |= 1 << (n - 1); else numer += denom
    265 
    266 int32_t SkDivBits(int32_t numer, int32_t denom, int shift_bias) {
    267     SkASSERT(denom != 0);
    268     if (numer == 0) {
    269         return 0;
    270     }
    271 
    272     // make numer and denom positive, and sign hold the resulting sign
    273     int32_t sign = SkExtractSign(numer ^ denom);
    274     numer = SkAbs32(numer);
    275     denom = SkAbs32(denom);
    276 
    277     int nbits = SkCLZ(numer) - 1;
    278     int dbits = SkCLZ(denom) - 1;
    279     int bits = shift_bias - nbits + dbits;
    280 
    281     if (bits < 0) {  // answer will underflow
    282         return 0;
    283     }
    284     if (bits > 31) {  // answer will overflow
    285         return SkApplySign(SK_MaxS32, sign);
    286     }
    287 
    288     denom <<= dbits;
    289     numer <<= nbits;
    290 
    291     SkFixed result = 0;
    292 
    293     // do the first one
    294     if ((numer -= denom) >= 0) {
    295         result = 1;
    296     } else {
    297         numer += denom;
    298     }
    299 
    300     // Now fall into our switch statement if there are more bits to compute
    301     if (bits > 0) {
    302         // make room for the rest of the answer bits
    303         result <<= bits;
    304         switch (bits) {
    305             DIVBITS_ITER(31); DIVBITS_ITER(30); DIVBITS_ITER(29);
    306             DIVBITS_ITER(28); DIVBITS_ITER(27); DIVBITS_ITER(26);
    307             DIVBITS_ITER(25); DIVBITS_ITER(24); DIVBITS_ITER(23);
    308             DIVBITS_ITER(22); DIVBITS_ITER(21); DIVBITS_ITER(20);
    309             DIVBITS_ITER(19); DIVBITS_ITER(18); DIVBITS_ITER(17);
    310             DIVBITS_ITER(16); DIVBITS_ITER(15); DIVBITS_ITER(14);
    311             DIVBITS_ITER(13); DIVBITS_ITER(12); DIVBITS_ITER(11);
    312             DIVBITS_ITER(10); DIVBITS_ITER( 9); DIVBITS_ITER( 8);
    313             DIVBITS_ITER( 7); DIVBITS_ITER( 6); DIVBITS_ITER( 5);
    314             DIVBITS_ITER( 4); DIVBITS_ITER( 3); DIVBITS_ITER( 2);
    315             // we merge these last two together, makes GCC make better ARM
    316             default:
    317             DIVBITS_ITER( 1);
    318         }
    319     }
    320 
    321     if (result < 0) {
    322         result = SK_MaxS32;
    323     }
    324     return SkApplySign(result, sign);
    325 }
    326 
    327 /*  mod(float numer, float denom) seems to always return the sign
    328     of the numer, so that's what we do too
    329 */
    330 SkFixed SkFixedMod(SkFixed numer, SkFixed denom) {
    331     int sn = SkExtractSign(numer);
    332     int sd = SkExtractSign(denom);
    333 
    334     numer = SkApplySign(numer, sn);
    335     denom = SkApplySign(denom, sd);
    336 
    337     if (numer < denom) {
    338         return SkApplySign(numer, sn);
    339     } else if (numer == denom) {
    340         return 0;
    341     } else {
    342         SkFixed div = SkFixedDiv(numer, denom);
    343         return SkApplySign(SkFixedMul(denom, div & 0xFFFF), sn);
    344     }
    345 }
    346 
    347 /* www.worldserver.com/turk/computergraphics/FixedSqrt.pdf
    348 */
    349 int32_t SkSqrtBits(int32_t x, int count) {
    350     SkASSERT(x >= 0 && count > 0 && (unsigned)count <= 30);
    351 
    352     uint32_t    root = 0;
    353     uint32_t    remHi = 0;
    354     uint32_t    remLo = x;
    355 
    356     do {
    357         root <<= 1;
    358 
    359         remHi = (remHi<<2) | (remLo>>30);
    360         remLo <<= 2;
    361 
    362         uint32_t testDiv = (root << 1) + 1;
    363         if (remHi >= testDiv) {
    364             remHi -= testDiv;
    365             root++;
    366         }
    367     } while (--count >= 0);
    368 
    369     return root;
    370 }
    371 
    372 int32_t SkCubeRootBits(int32_t value, int bits) {
    373     SkASSERT(bits > 0);
    374 
    375     int sign = SkExtractSign(value);
    376     value = SkApplySign(value, sign);
    377 
    378     uint32_t root = 0;
    379     uint32_t curr = (uint32_t)value >> 30;
    380     value <<= 2;
    381 
    382     do {
    383         root <<= 1;
    384         uint32_t guess = root * root + root;
    385         guess = (guess << 1) + guess;   // guess *= 3
    386         if (guess < curr) {
    387             curr -= guess + 1;
    388             root |= 1;
    389         }
    390         curr = (curr << 3) | ((uint32_t)value >> 29);
    391         value <<= 3;
    392     } while (--bits);
    393 
    394     return SkApplySign(root, sign);
    395 }
    396 
    397 SkFixed SkFixedMean(SkFixed a, SkFixed b) {
    398     Sk64 tmp;
    399 
    400     tmp.setMul(a, b);
    401     return tmp.getSqrt();
    402 }
    403 
    404 ///////////////////////////////////////////////////////////////////////////////
    405 
    406 #ifdef SK_SCALAR_IS_FLOAT
    407 float SkScalarSinCos(float radians, float* cosValue) {
    408     float sinValue = sk_float_sin(radians);
    409 
    410     if (cosValue) {
    411         *cosValue = sk_float_cos(radians);
    412         if (SkScalarNearlyZero(*cosValue)) {
    413             *cosValue = 0;
    414         }
    415     }
    416 
    417     if (SkScalarNearlyZero(sinValue)) {
    418         sinValue = 0;
    419     }
    420     return sinValue;
    421 }
    422 #endif
    423 
    424 #define INTERP_SINTABLE
    425 #define BUILD_TABLE_AT_RUNTIMEx
    426 
    427 #define kTableSize  256
    428 
    429 #ifdef BUILD_TABLE_AT_RUNTIME
    430     static uint16_t gSkSinTable[kTableSize];
    431 
    432     static void build_sintable(uint16_t table[]) {
    433         for (int i = 0; i < kTableSize; i++) {
    434             double  rad = i * 3.141592653589793 / (2*kTableSize);
    435             double  val = sin(rad);
    436             int     ival = (int)(val * SK_Fixed1);
    437             table[i] = SkToU16(ival);
    438         }
    439     }
    440 #else
    441     #include "SkSinTable.h"
    442 #endif
    443 
    444 #define SK_Fract1024SizeOver2PI     0x28BE60    /* floatToFract(1024 / 2PI) */
    445 
    446 #ifdef INTERP_SINTABLE
    447 static SkFixed interp_table(const uint16_t table[], int index, int partial255) {
    448     SkASSERT((unsigned)index < kTableSize);
    449     SkASSERT((unsigned)partial255 <= 255);
    450 
    451     SkFixed lower = table[index];
    452     SkFixed upper = (index == kTableSize - 1) ? SK_Fixed1 : table[index + 1];
    453 
    454     SkASSERT(lower < upper);
    455     SkASSERT(lower >= 0);
    456     SkASSERT(upper <= SK_Fixed1);
    457 
    458     partial255 += (partial255 >> 7);
    459     return lower + ((upper - lower) * partial255 >> 8);
    460 }
    461 #endif
    462 
    463 SkFixed SkFixedSinCos(SkFixed radians, SkFixed* cosValuePtr) {
    464     SkASSERT(SK_ARRAY_COUNT(gSkSinTable) == kTableSize);
    465 
    466 #ifdef BUILD_TABLE_AT_RUNTIME
    467     static bool gFirstTime = true;
    468     if (gFirstTime) {
    469         build_sintable(gSinTable);
    470         gFirstTime = false;
    471     }
    472 #endif
    473 
    474     // make radians positive
    475     SkFixed sinValue, cosValue;
    476     int32_t cosSign = 0;
    477     int32_t sinSign = SkExtractSign(radians);
    478     radians = SkApplySign(radians, sinSign);
    479     // scale it to 0...1023 ...
    480 
    481 #ifdef INTERP_SINTABLE
    482     radians = SkMulDiv(radians, 2 * kTableSize * 256, SK_FixedPI);
    483     int findex = radians & (kTableSize * 256 - 1);
    484     int index = findex >> 8;
    485     int partial = findex & 255;
    486     sinValue = interp_table(gSkSinTable, index, partial);
    487 
    488     findex = kTableSize * 256 - findex - 1;
    489     index = findex >> 8;
    490     partial = findex & 255;
    491     cosValue = interp_table(gSkSinTable, index, partial);
    492 
    493     int quad = ((unsigned)radians / (kTableSize * 256)) & 3;
    494 #else
    495     radians = SkMulDiv(radians, 2 * kTableSize, SK_FixedPI);
    496     int     index = radians & (kTableSize - 1);
    497 
    498     if (index == 0) {
    499         sinValue = 0;
    500         cosValue = SK_Fixed1;
    501     } else {
    502         sinValue = gSkSinTable[index];
    503         cosValue = gSkSinTable[kTableSize - index];
    504     }
    505     int quad = ((unsigned)radians / kTableSize) & 3;
    506 #endif
    507 
    508     if (quad & 1) {
    509         SkTSwap<SkFixed>(sinValue, cosValue);
    510     }
    511     if (quad & 2) {
    512         sinSign = ~sinSign;
    513     }
    514     if (((quad - 1) & 2) == 0) {
    515         cosSign = ~cosSign;
    516     }
    517 
    518     // restore the sign for negative angles
    519     sinValue = SkApplySign(sinValue, sinSign);
    520     cosValue = SkApplySign(cosValue, cosSign);
    521 
    522 #ifdef SK_DEBUG
    523     if (1) {
    524         SkFixed sin2 = SkFixedMul(sinValue, sinValue);
    525         SkFixed cos2 = SkFixedMul(cosValue, cosValue);
    526         int diff = cos2 + sin2 - SK_Fixed1;
    527         SkASSERT(SkAbs32(diff) <= 7);
    528     }
    529 #endif
    530 
    531     if (cosValuePtr) {
    532         *cosValuePtr = cosValue;
    533     }
    534     return sinValue;
    535 }
    536 
    537 ///////////////////////////////////////////////////////////////////////////////
    538 
    539 SkFixed SkFixedTan(SkFixed radians) { return SkCordicTan(radians); }
    540 SkFixed SkFixedASin(SkFixed x) { return SkCordicASin(x); }
    541 SkFixed SkFixedACos(SkFixed x) { return SkCordicACos(x); }
    542 SkFixed SkFixedATan2(SkFixed y, SkFixed x) { return SkCordicATan2(y, x); }
    543 SkFixed SkFixedExp(SkFixed x) { return SkCordicExp(x); }
    544 SkFixed SkFixedLog(SkFixed x) { return SkCordicLog(x); }
    545 
    546