michael@0: /*
michael@0:  * Copyright 2008 The Android Open Source Project
michael@0:  *
michael@0:  * Use of this source code is governed by a BSD-style license that can be
michael@0:  * found in the LICENSE file.
michael@0:  */
michael@0: 
michael@0: #include "SkMathPriv.h"
michael@0: #include "SkFloatBits.h"
michael@0: #include "SkFloatingPoint.h"
michael@0: #include "SkScalar.h"
michael@0: 
michael@0: const uint32_t gIEEENotANumber = 0x7FFFFFFF;
michael@0: const uint32_t gIEEEInfinity = 0x7F800000;
michael@0: const uint32_t gIEEENegativeInfinity = 0xFF800000;
michael@0: 
michael@0: #define sub_shift(zeros, x, n)  \
michael@0:     zeros -= n;                 \
michael@0:     x >>= n
michael@0: 
michael@0: int SkCLZ_portable(uint32_t x) {
michael@0:     if (x == 0) {
michael@0:         return 32;
michael@0:     }
michael@0: 
michael@0:     int zeros = 31;
michael@0:     if (x & 0xFFFF0000) {
michael@0:         sub_shift(zeros, x, 16);
michael@0:     }
michael@0:     if (x & 0xFF00) {
michael@0:         sub_shift(zeros, x, 8);
michael@0:     }
michael@0:     if (x & 0xF0) {
michael@0:         sub_shift(zeros, x, 4);
michael@0:     }
michael@0:     if (x & 0xC) {
michael@0:         sub_shift(zeros, x, 2);
michael@0:     }
michael@0:     if (x & 0x2) {
michael@0:         sub_shift(zeros, x, 1);
michael@0:     }
michael@0: 
michael@0:     return zeros;
michael@0: }
michael@0: 
michael@0: int32_t SkMulDiv(int32_t numer1, int32_t numer2, int32_t denom) {
michael@0:     SkASSERT(denom);
michael@0: 
michael@0:     int64_t tmp = sk_64_mul(numer1, numer2) / denom;
michael@0:     return sk_64_asS32(tmp);
michael@0: }
michael@0: 
michael@0: SkFixed SkFixedMul_portable(SkFixed a, SkFixed b) {
michael@0: #if defined(SkLONGLONG)
michael@0:     return static_cast<SkFixed>((int64_t)a * b >> 16);
michael@0: #else
michael@0:     int sa = SkExtractSign(a);
michael@0:     int sb = SkExtractSign(b);
michael@0:     // now make them positive
michael@0:     a = SkApplySign(a, sa);
michael@0:     b = SkApplySign(b, sb);
michael@0: 
michael@0:     uint32_t    ah = a >> 16;
michael@0:     uint32_t    al = a & 0xFFFF;
michael@0:     uint32_t bh = b >> 16;
michael@0:     uint32_t bl = b & 0xFFFF;
michael@0: 
michael@0:     uint32_t R = ah * b + al * bh + (al * bl >> 16);
michael@0: 
michael@0:     return SkApplySign(R, sa ^ sb);
michael@0: #endif
michael@0: }
michael@0: 
michael@0: ///////////////////////////////////////////////////////////////////////////////
michael@0: 
michael@0: #define DIVBITS_ITER(n)                                 \
michael@0:     case n:                                             \
michael@0:         if ((numer = (numer << 1) - denom) >= 0)        \
michael@0:             result |= 1 << (n - 1); else numer += denom
michael@0: 
michael@0: int32_t SkDivBits(int32_t numer, int32_t denom, int shift_bias) {
michael@0:     SkASSERT(denom != 0);
michael@0:     if (numer == 0) {
michael@0:         return 0;
michael@0:     }
michael@0: 
michael@0:     // make numer and denom positive, and sign hold the resulting sign
michael@0:     int32_t sign = SkExtractSign(numer ^ denom);
michael@0:     numer = SkAbs32(numer);
michael@0:     denom = SkAbs32(denom);
michael@0: 
michael@0:     int nbits = SkCLZ(numer) - 1;
michael@0:     int dbits = SkCLZ(denom) - 1;
michael@0:     int bits = shift_bias - nbits + dbits;
michael@0: 
michael@0:     if (bits < 0) {  // answer will underflow
michael@0:         return 0;
michael@0:     }
michael@0:     if (bits > 31) {  // answer will overflow
michael@0:         return SkApplySign(SK_MaxS32, sign);
michael@0:     }
michael@0: 
michael@0:     denom <<= dbits;
michael@0:     numer <<= nbits;
michael@0: 
michael@0:     SkFixed result = 0;
michael@0: 
michael@0:     // do the first one
michael@0:     if ((numer -= denom) >= 0) {
michael@0:         result = 1;
michael@0:     } else {
michael@0:         numer += denom;
michael@0:     }
michael@0: 
michael@0:     // Now fall into our switch statement if there are more bits to compute
michael@0:     if (bits > 0) {
michael@0:         // make room for the rest of the answer bits
michael@0:         result <<= bits;
michael@0:         switch (bits) {
michael@0:             DIVBITS_ITER(31); DIVBITS_ITER(30); DIVBITS_ITER(29);
michael@0:             DIVBITS_ITER(28); DIVBITS_ITER(27); DIVBITS_ITER(26);
michael@0:             DIVBITS_ITER(25); DIVBITS_ITER(24); DIVBITS_ITER(23);
michael@0:             DIVBITS_ITER(22); DIVBITS_ITER(21); DIVBITS_ITER(20);
michael@0:             DIVBITS_ITER(19); DIVBITS_ITER(18); DIVBITS_ITER(17);
michael@0:             DIVBITS_ITER(16); DIVBITS_ITER(15); DIVBITS_ITER(14);
michael@0:             DIVBITS_ITER(13); DIVBITS_ITER(12); DIVBITS_ITER(11);
michael@0:             DIVBITS_ITER(10); DIVBITS_ITER( 9); DIVBITS_ITER( 8);
michael@0:             DIVBITS_ITER( 7); DIVBITS_ITER( 6); DIVBITS_ITER( 5);
michael@0:             DIVBITS_ITER( 4); DIVBITS_ITER( 3); DIVBITS_ITER( 2);
michael@0:             // we merge these last two together, makes GCC make better ARM
michael@0:             default:
michael@0:             DIVBITS_ITER( 1);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     if (result < 0) {
michael@0:         result = SK_MaxS32;
michael@0:     }
michael@0:     return SkApplySign(result, sign);
michael@0: }
michael@0: 
michael@0: /* www.worldserver.com/turk/computergraphics/FixedSqrt.pdf
michael@0: */
michael@0: int32_t SkSqrtBits(int32_t x, int count) {
michael@0:     SkASSERT(x >= 0 && count > 0 && (unsigned)count <= 30);
michael@0: 
michael@0:     uint32_t    root = 0;
michael@0:     uint32_t    remHi = 0;
michael@0:     uint32_t    remLo = x;
michael@0: 
michael@0:     do {
michael@0:         root <<= 1;
michael@0: 
michael@0:         remHi = (remHi<<2) | (remLo>>30);
michael@0:         remLo <<= 2;
michael@0: 
michael@0:         uint32_t testDiv = (root << 1) + 1;
michael@0:         if (remHi >= testDiv) {
michael@0:             remHi -= testDiv;
michael@0:             root++;
michael@0:         }
michael@0:     } while (--count >= 0);
michael@0: 
michael@0:     return root;
michael@0: }
michael@0: 
michael@0: ///////////////////////////////////////////////////////////////////////////////
michael@0: 
michael@0: float SkScalarSinCos(float radians, float* cosValue) {
michael@0:     float sinValue = sk_float_sin(radians);
michael@0: 
michael@0:     if (cosValue) {
michael@0:         *cosValue = sk_float_cos(radians);
michael@0:         if (SkScalarNearlyZero(*cosValue)) {
michael@0:             *cosValue = 0;
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     if (SkScalarNearlyZero(sinValue)) {
michael@0:         sinValue = 0;
michael@0:     }
michael@0:     return sinValue;
michael@0: }
michael@0: 
michael@0: #define INTERP_SINTABLE
michael@0: #define BUILD_TABLE_AT_RUNTIMEx
michael@0: 
michael@0: #define kTableSize  256
michael@0: 
michael@0: #ifdef BUILD_TABLE_AT_RUNTIME
michael@0:     static uint16_t gSkSinTable[kTableSize];
michael@0: 
michael@0:     static void build_sintable(uint16_t table[]) {
michael@0:         for (int i = 0; i < kTableSize; i++) {
michael@0:             double  rad = i * 3.141592653589793 / (2*kTableSize);
michael@0:             double  val = sin(rad);
michael@0:             int     ival = (int)(val * SK_Fixed1);
michael@0:             table[i] = SkToU16(ival);
michael@0:         }
michael@0:     }
michael@0: #else
michael@0:     #include "SkSinTable.h"
michael@0: #endif
michael@0: 
michael@0: #define SK_Fract1024SizeOver2PI     0x28BE60    /* floatToFract(1024 / 2PI) */
michael@0: 
michael@0: #ifdef INTERP_SINTABLE
michael@0: static SkFixed interp_table(const uint16_t table[], int index, int partial255) {
michael@0:     SkASSERT((unsigned)index < kTableSize);
michael@0:     SkASSERT((unsigned)partial255 <= 255);
michael@0: 
michael@0:     SkFixed lower = table[index];
michael@0:     SkFixed upper = (index == kTableSize - 1) ? SK_Fixed1 : table[index + 1];
michael@0: 
michael@0:     SkASSERT(lower < upper);
michael@0:     SkASSERT(lower >= 0);
michael@0:     SkASSERT(upper <= SK_Fixed1);
michael@0: 
michael@0:     partial255 += (partial255 >> 7);
michael@0:     return lower + ((upper - lower) * partial255 >> 8);
michael@0: }
michael@0: #endif
michael@0: 
michael@0: SkFixed SkFixedSinCos(SkFixed radians, SkFixed* cosValuePtr) {
michael@0:     SkASSERT(SK_ARRAY_COUNT(gSkSinTable) == kTableSize);
michael@0: 
michael@0: #ifdef BUILD_TABLE_AT_RUNTIME
michael@0:     static bool gFirstTime = true;
michael@0:     if (gFirstTime) {
michael@0:         build_sintable(gSinTable);
michael@0:         gFirstTime = false;
michael@0:     }
michael@0: #endif
michael@0: 
michael@0:     // make radians positive
michael@0:     SkFixed sinValue, cosValue;
michael@0:     int32_t cosSign = 0;
michael@0:     int32_t sinSign = SkExtractSign(radians);
michael@0:     radians = SkApplySign(radians, sinSign);
michael@0:     // scale it to 0...1023 ...
michael@0: 
michael@0: #ifdef INTERP_SINTABLE
michael@0:     radians = SkMulDiv(radians, 2 * kTableSize * 256, SK_FixedPI);
michael@0:     int findex = radians & (kTableSize * 256 - 1);
michael@0:     int index = findex >> 8;
michael@0:     int partial = findex & 255;
michael@0:     sinValue = interp_table(gSkSinTable, index, partial);
michael@0: 
michael@0:     findex = kTableSize * 256 - findex - 1;
michael@0:     index = findex >> 8;
michael@0:     partial = findex & 255;
michael@0:     cosValue = interp_table(gSkSinTable, index, partial);
michael@0: 
michael@0:     int quad = ((unsigned)radians / (kTableSize * 256)) & 3;
michael@0: #else
michael@0:     radians = SkMulDiv(radians, 2 * kTableSize, SK_FixedPI);
michael@0:     int     index = radians & (kTableSize - 1);
michael@0: 
michael@0:     if (index == 0) {
michael@0:         sinValue = 0;
michael@0:         cosValue = SK_Fixed1;
michael@0:     } else {
michael@0:         sinValue = gSkSinTable[index];
michael@0:         cosValue = gSkSinTable[kTableSize - index];
michael@0:     }
michael@0:     int quad = ((unsigned)radians / kTableSize) & 3;
michael@0: #endif
michael@0: 
michael@0:     if (quad & 1) {
michael@0:         SkTSwap<SkFixed>(sinValue, cosValue);
michael@0:     }
michael@0:     if (quad & 2) {
michael@0:         sinSign = ~sinSign;
michael@0:     }
michael@0:     if (((quad - 1) & 2) == 0) {
michael@0:         cosSign = ~cosSign;
michael@0:     }
michael@0: 
michael@0:     // restore the sign for negative angles
michael@0:     sinValue = SkApplySign(sinValue, sinSign);
michael@0:     cosValue = SkApplySign(cosValue, cosSign);
michael@0: 
michael@0: #ifdef SK_DEBUG
michael@0:     if (1) {
michael@0:         SkFixed sin2 = SkFixedMul(sinValue, sinValue);
michael@0:         SkFixed cos2 = SkFixedMul(cosValue, cosValue);
michael@0:         int diff = cos2 + sin2 - SK_Fixed1;
michael@0:         SkASSERT(SkAbs32(diff) <= 7);
michael@0:     }
michael@0: #endif
michael@0: 
michael@0:     if (cosValuePtr) {
michael@0:         *cosValuePtr = cosValue;
michael@0:     }
michael@0:     return sinValue;
michael@0: }