Changeset 899


Ignore:
Timestamp:
Sep 4, 2011, 2:28:52 PM (8 years ago)
Author:
sam
Message:

core: implement accelerated lol_sincos() and lol_tan().

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/core.h

    r897 r899  
    2424#   define LOL_FEATURE_CHEAP_BRANCHES
    2525#endif
     26
     27// Optimisation helpers
     28#define __likely(x)   __builtin_expect(!!(x), 1)
     29#define __unlikely(x) __builtin_expect(!!(x), 0)
    2630
    2731// Base types
  • trunk/src/trig.cpp

    r898 r899  
    6666};
    6767
     68/* These coefficients use Sloane’s http://oeis.org/A002430 and
     69 * http://oeis.org/A036279 sequences for the Taylor series of tan().
     70 * Note: the last value should be 443861162*pi^18/1856156927625, ie.
     71 * 2.12485922978838540352881e5, but we tweak it in order to get
     72 * sub 1e-11 precision in a larger range. */
     73static const double TC[] =
     74{
     75    3.28986813369645287294483e0, // pi^2/3
     76    1.29878788045336582981920e1, // 2*pi^4/15
     77    5.18844961612069061254404e1, // 17*pi^6/315
     78    2.07509320280908496804928e2, // 62*pi^8/2835
     79    8.30024701695986756361561e2, // 1382*pi^10/155925
     80    3.32009324029001216460018e3, // 21844*pi^12/6081075
     81    1.32803704909665483598490e4, // 929569*pi^14/638512875
     82    5.31214808666037709352112e4, // 6404582*pi^16/10854718875
     83    2.373e5, // XXX: last value tweaked to improve precision
     84    //2.12485922978838540352881e5, // 443861162*pi^18/1856156927625
     85};
     86
    6887/* Custom intrinsics */
    6988#define INLINEATTR __attribute__((always_inline))
     
    237256#endif
    238257
     258    /* Compute a Tailor series for sin() and combine sign information. */
    239259    sign *= (x >= 0.0) ? PI : NEG_PI;
    240260
     
    314334}
    315335
     336void lol_sincos(double x, double *sinx, double *cosx)
     337{
     338    double absx = lol_fabs(x * INV_PI);
     339
     340#if defined LOL_FEATURE_CHEAP_BRANCHES
     341    if (absx < QUARTER)
     342    {
     343        double x2 = absx * absx;
     344        double x4 = x2 * x2;
     345
     346        double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
     347        double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
     348        double taylors = subs2 * x2 + subs1;
     349        *sinx = x * taylors;
     350
     351        double subc1 = (CC[5] * x4 + CC[3]) * x4 + CC[1];
     352        double subc2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
     353        double taylorc = (subc1 * x2 + subc2) * x2 + ONE;
     354        *cosx = taylorc;
     355
     356        return;
     357    }
     358#endif
     359
     360#if defined __CELLOS_LV2__
     361    double num_cycles = lol_round(absx);
     362    double is_even = lol_trunc(num_cycles * HALF) - (num_cycles * HALF);
     363
     364    double sin_sign = lol_fsel(x, PI, NEG_PI);
     365    sin_sign = lol_fsel(is_even, sin_sign, -sin_sign);
     366    double cos_sign = lol_fsel(is_even, ONE, NEG_ONE);
     367#else
     368    double num_cycles = absx + TWO_EXP_52;
     369    __asm__("" : "+m" (num_cycles)); num_cycles -= TWO_EXP_52;
     370
     371    double is_even = TWO * num_cycles - ONE;
     372    __asm__("" : "+m" (is_even)); is_even += TWO_EXP_54;
     373    __asm__("" : "+m" (is_even)); is_even -= TWO_EXP_54;
     374    __asm__("" : "+m" (is_even));
     375    is_even -= TWO * num_cycles - ONE;
     376    double sin_sign = is_even;
     377    double cos_sign = is_even;
     378#endif
     379    absx -= num_cycles;
     380
     381#if defined LOL_FEATURE_VERY_CHEAP_BRANCHES
     382    if (lol_fabs(absx) > QUARTER)
     383    {
     384        sign = (x * absx >= 0.0) ? sign : -sign;
     385
     386        double x1 = HALF - lol_fabs(absx);
     387        double x2 = x1 * x1;
     388        double x4 = x2 * x2;
     389
     390        double subs1 = ((CC[5] * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
     391        double subs2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
     392        double taylors = subs2 * x2 + subs1;
     393        *sinx = taylors * sign;
     394
     395        double subc1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
     396        double subc2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
     397        double taylorc = subc2 * x2 + subc1;
     398        *cosx = x1 * taylorc * sign * PI;
     399
     400        return;
     401    }
     402#endif
     403
     404    sin_sign *= (x >= 0.0) ? PI : NEG_PI;
     405
     406    double x2 = absx * absx;
     407    double x4 = x2 * x2;
     408#if defined LOL_FEATURE_VERY_CHEAP_BRANCHES
     409    double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
     410    double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
     411    double subc1 = ((CC[5] * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
     412    double subc2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
     413#else
     414    double subs1 = (((SC[7] * x4 + SC[5]) * x4 + SC[3]) * x4 + SC[1]) * x4 + ONE;
     415    double subs2 = ((SC[6] * x4 + SC[4]) * x4 + SC[2]) * x4 + SC[0];
     416    double subc1 = (((CC[7] * x4 + CC[5]) * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
     417    double subc2 = ((CC[6] * x4 + CC[4]) * x4 + CC[2]) * x4 + CC[0];
     418#endif
     419    double taylors = subs2 * x2 + subs1;
     420    *sinx = absx * taylors * sin_sign;
     421
     422    double taylorc = subc2 * x2 + subc1;
     423    *cosx = taylorc * cos_sign;
     424}
     425
     426void lol_sincos(float x, float *sinx, float *cosx)
     427{
     428    double x2 = static_cast<double>(x);
     429    double s2, c2;
     430    lol_sincos(x2, &s2, &c2);
     431    *sinx = static_cast<float>(s2);
     432    *cosx = static_cast<float>(c2);
     433}
     434
     435double lol_tan(double x)
     436{
     437#if defined LOL_FEATURE_CHEAP_BRANCHES
     438    double absx = lol_fabs(x * INV_PI);
     439
     440    /* This value was determined empirically to ensure an error of no
     441     * more than x * 1e-11 in this range. */
     442    if (absx < 0.163)
     443    {
     444        double x2 = absx * absx;
     445        double x4 = x2 * x2;
     446        double sub1 = (((TC[7] * x4 + TC[5]) * x4
     447                           + TC[3]) * x4 + TC[1]) * x4 + ONE;
     448        double sub2 = (((TC[8] * x4 + TC[6]) * x4
     449                           + TC[4]) * x4 + TC[2]) * x4 + TC[0];
     450        double taylor = sub2 * x2 + sub1;
     451        return x * taylor;
     452    }
     453#endif
     454
     455    double sinx, cosx;
     456    lol_sincos(x, &sinx, &cosx);
     457
     458    /* Ensure cosx isn't zero. FIXME: we lose the cosx sign here. */
     459    double absc = lol_fabs(cosx);
     460#if defined __CELLOS_LV2__
     461    double is_cos_not_zero = absc - VERY_SMALL_NUMBER;
     462    cosx = lol_fsel(is_cos_not_zero, cosx, VERY_SMALL_NUMBER);
     463    return __fdiv(sinx, cosx);
     464#else
     465    if (__unlikely(absc < VERY_SMALL_NUMBER))
     466        cosx = VERY_SMALL_NUMBER;
     467    return sinx / cosx;
     468#endif
     469}
     470
    316471} /* namespace lol */
    317472
  • trunk/test/lol-bench.cpp

    r897 r899  
    7575static void bench_trig(int mode)
    7676{
    77     float result[7] = { 0.0f };
     77    float result[12] = { 0.0f };
    7878    Timer timer;
    7979
     
    8181    float *pf = new float[TRIG_TABLE_SIZE];
    8282    float *pf2 = new float[TRIG_TABLE_SIZE];
     83    float *pf3 = new float[TRIG_TABLE_SIZE];
    8384
    8485    for (size_t run = 0; run < TRIG_RUNS; run++)
     
    144145        result[5] += timer.GetMs();
    145146
     147        /* Sin & cos */
     148        timer.GetMs();
     149        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
     150        {
     151            pf2[i] = __builtin_sinf(pf[i]);
     152            pf3[i] = __builtin_cosf(pf[i]);
     153        }
     154        result[6] += timer.GetMs();
     155
     156        /* Fast sin & cos */
     157        timer.GetMs();
     158        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
     159        {
     160#if defined HAVE_FASTMATH_H
     161            pf2[i] = f_sinf(pf[i]);
     162            pf3[i] = f_cosf(pf[i]);
     163#else
     164            pf2[i] = sinf(pf[i]);
     165            pf3[i] = cosf(pf[i]);
     166#endif
     167        }
     168        result[7] += timer.GetMs();
     169
     170        /* Lol sincos */
     171        timer.GetMs();
     172        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
     173            lol_sincos(pf[i], &pf2[i], &pf3[i]);
     174        result[8] += timer.GetMs();
     175
    146176        /* Tan */
    147177        timer.GetMs();
    148178        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
    149179            pf2[i] = __builtin_tanf(pf[i]);
    150         result[6] += timer.GetMs();
     180        result[9] += timer.GetMs();
     181
     182        /* Fast tan */
     183        timer.GetMs();
     184        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
     185#if defined HAVE_FASTMATH_H
     186            pf2[i] = f_tanf(pf[i]);
     187#else
     188            pf2[i] = tanf(pf[i]);
     189#endif
     190        result[10] += timer.GetMs();
     191
     192        /* Lol tan */
     193        timer.GetMs();
     194        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
     195            pf2[i] = lol_tan(pf[i]);
     196        result[11] += timer.GetMs();
    151197    }
    152198
    153199    delete[] pf;
    154200    delete[] pf2;
     201    delete[] pf3;
    155202
    156203    for (size_t i = 0; i < sizeof(result) / sizeof(*result); i++)
    157204        result[i] *= 1000000.0f / (TRIG_TABLE_SIZE * TRIG_RUNS);
    158205
    159     Log::Info("                          ns/elem\n");
    160     Log::Info("float = sinf(float)      %7.3f\n", result[0]);
    161     Log::Info("float = fastsinf(float)  %7.3f\n", result[1]);
    162     Log::Info("float = lol_sinf(float)  %7.3f\n", result[2]);
    163     Log::Info("float = cosf(float)      %7.3f\n", result[3]);
    164     Log::Info("float = fastcosf(float)  %7.3f\n", result[4]);
    165     Log::Info("float = lol_cosf(float)  %7.3f\n", result[5]);
    166     Log::Info("float = tanf(float)      %7.3f\n", result[6]);
     206    Log::Info("                              ns/elem\n");
     207    Log::Info("float = sinf(float)          %7.3f\n", result[0]);
     208    Log::Info("float = f_sinf(float)        %7.3f\n", result[1]);
     209    Log::Info("float = lol_sin(float)       %7.3f\n", result[2]);
     210    Log::Info("float = cosf(float)          %7.3f\n", result[3]);
     211    Log::Info("float = f_cosf(float)        %7.3f\n", result[4]);
     212    Log::Info("float = lol_cos(float)       %7.3f\n", result[5]);
     213    Log::Info("float = sinf,cosf(float)     %7.3f\n", result[6]);
     214    Log::Info("float = f_sinf,f_cosf(float) %7.3f\n", result[7]);
     215    Log::Info("float = lol_sincos(float)    %7.3f\n", result[8]);
     216    Log::Info("float = tanf(float)          %7.3f\n", result[9]);
     217    Log::Info("float = f_tanf(float)        %7.3f\n", result[10]);
     218    Log::Info("float = lol_tanf(float)      %7.3f\n", result[11]);
    167219}
    168220
  • trunk/test/trig.cpp

    r897 r899  
    7171            CPPUNIT_ASSERT(fabs(a - b) <= fabs(f) * 1e-11);
    7272        }
     73
     74        for (int i = -10000; i < 10000; i++)
     75        {
     76            double f = (double)i * (1.0 / 1000.0);
     77            double a1 = __builtin_sin(f);
     78            double a2 = __builtin_cos(f);
     79            double b1, b2;
     80            lol_sincos(f, &b1, &b2);
     81            CPPUNIT_ASSERT(fabs(a1 - b1) <= fabs(f) * 1e-11);
     82            CPPUNIT_ASSERT(fabs(a2 - b2) <= fabs(f) * 1e-11);
     83        }
     84
     85        for (int i = -10000; i < 10000; i++)
     86        {
     87            double f = (double)i * (1.0 / 100000.0);
     88            double a1 = __builtin_sin(f);
     89            double a2 = __builtin_cos(f);
     90            double b1, b2;
     91            lol_sincos(f, &b1, &b2);
     92            CPPUNIT_ASSERT(fabs(a1 - b1) <= fabs(f) * 1e-11);
     93            CPPUNIT_ASSERT(fabs(a2 - b2) <= fabs(f) * 1e-11);
     94        }
     95
     96        for (int i = -10000; i < 10000; i++)
     97        {
     98            double f = (double)i * (1.0 / 1000.0);
     99            double a = __builtin_tan(f);
     100            double b = lol_tan(f);
     101            CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     102        }
     103
     104        for (int i = -10000; i < 10000; i++)
     105        {
     106            double f = (double)i * (1.0 / 100000.0);
     107            double a = __builtin_tan(f);
     108            double b = lol_tan(f);
     109            CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     110        }
    73111    }
    74112};
Note: See TracChangeset for help on using the changeset viewer.