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

src/core.h

@@ -24,6 +24,10 @@
 #   define LOL_FEATURE_CHEAP_BRANCHES
 #endif
 
+// Optimisation helpers
+#define __likely(x)   __builtin_expect(!!(x), 1)
+#define __unlikely(x) __builtin_expect(!!(x), 0)
+
 // Base types
 #include "trig.h"
 #include "half.h"

src/trig.cpp

@@ -65,6 +65,25 @@ static const double CC[] =
     +4.3030695870329470072978e-6, // pi^16/16!
 };
 
+/* These coefficients use Sloane’s http://oeis.org/A002430 and
+ * http://oeis.org/A036279 sequences for the Taylor series of tan().
+ * Note: the last value should be 443861162*pi^18/1856156927625, ie.
+ * 2.12485922978838540352881e5, but we tweak it in order to get
+ * sub 1e-11 precision in a larger range. */
+static const double TC[] =
+{
+    3.28986813369645287294483e0, // pi^2/3
+    1.29878788045336582981920e1, // 2*pi^4/15
+    5.18844961612069061254404e1, // 17*pi^6/315
+    2.07509320280908496804928e2, // 62*pi^8/2835
+    8.30024701695986756361561e2, // 1382*pi^10/155925
+    3.32009324029001216460018e3, // 21844*pi^12/6081075
+    1.32803704909665483598490e4, // 929569*pi^14/638512875
+    5.31214808666037709352112e4, // 6404582*pi^16/10854718875
+    2.373e5, // XXX: last value tweaked to improve precision
+    //2.12485922978838540352881e5, // 443861162*pi^18/1856156927625
+};
+
 /* Custom intrinsics */
 #define INLINEATTR __attribute__((always_inline))
 
@@ -236,6 +255,7 @@ double lol_sin(double x)
     }
 #endif
 
+    /* Compute a Tailor series for sin() and combine sign information. */
     sign *= (x >= 0.0) ? PI : NEG_PI;
 
     double x2 = absx * absx;
@@ -313,5 +333,140 @@ double lol_cos(double x)
     return taylor * sign;
 }
 
+void lol_sincos(double x, double *sinx, double *cosx)
+{
+    double absx = lol_fabs(x * INV_PI);
+
+#if defined LOL_FEATURE_CHEAP_BRANCHES
+    if (absx < QUARTER)
+    {
+        double x2 = absx * absx;
+        double x4 = x2 * x2;
+
+        double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
+        double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
+        double taylors = subs2 * x2 + subs1;
+        *sinx = x * taylors;
+
+        double subc1 = (CC[5] * x4 + CC[3]) * x4 + CC[1];
+        double subc2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
+        double taylorc = (subc1 * x2 + subc2) * x2 + ONE;
+        *cosx = taylorc;
+
+        return;
+    }
+#endif
+
+#if defined __CELLOS_LV2__
+    double num_cycles = lol_round(absx);
+    double is_even = lol_trunc(num_cycles * HALF) - (num_cycles * HALF);
+
+    double sin_sign = lol_fsel(x, PI, NEG_PI);
+    sin_sign = lol_fsel(is_even, sin_sign, -sin_sign);
+    double cos_sign = lol_fsel(is_even, ONE, NEG_ONE);
+#else
+    double num_cycles = absx + TWO_EXP_52;
+    __asm__("" : "+m" (num_cycles)); num_cycles -= TWO_EXP_52;
+
+    double is_even = TWO * num_cycles - ONE;
+    __asm__("" : "+m" (is_even)); is_even += TWO_EXP_54;
+    __asm__("" : "+m" (is_even)); is_even -= TWO_EXP_54;
+    __asm__("" : "+m" (is_even));
+    is_even -= TWO * num_cycles - ONE;
+    double sin_sign = is_even;
+    double cos_sign = is_even;
+#endif
+    absx -= num_cycles;
+
+#if defined LOL_FEATURE_VERY_CHEAP_BRANCHES
+    if (lol_fabs(absx) > QUARTER)
+    {
+        sign = (x * absx >= 0.0) ? sign : -sign;
+
+        double x1 = HALF - lol_fabs(absx);
+        double x2 = x1 * x1;
+        double x4 = x2 * x2;
+
+        double subs1 = ((CC[5] * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
+        double subs2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
+        double taylors = subs2 * x2 + subs1;
+        *sinx = taylors * sign;
+
+        double subc1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
+        double subc2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
+        double taylorc = subc2 * x2 + subc1;
+        *cosx = x1 * taylorc * sign * PI;
+
+        return;
+    }
+#endif
+
+    sin_sign *= (x >= 0.0) ? PI : NEG_PI;
+
+    double x2 = absx * absx;
+    double x4 = x2 * x2;
+#if defined LOL_FEATURE_VERY_CHEAP_BRANCHES
+    double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
+    double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
+    double subc1 = ((CC[5] * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
+    double subc2 = (CC[4] * x4 + CC[2]) * x4 + CC[0];
+#else
+    double subs1 = (((SC[7] * x4 + SC[5]) * x4 + SC[3]) * x4 + SC[1]) * x4 + ONE;
+    double subs2 = ((SC[6] * x4 + SC[4]) * x4 + SC[2]) * x4 + SC[0];
+    double subc1 = (((CC[7] * x4 + CC[5]) * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
+    double subc2 = ((CC[6] * x4 + CC[4]) * x4 + CC[2]) * x4 + CC[0];
+#endif
+    double taylors = subs2 * x2 + subs1;
+    *sinx = absx * taylors * sin_sign;
+
+    double taylorc = subc2 * x2 + subc1;
+    *cosx = taylorc * cos_sign;
+}
+
+void lol_sincos(float x, float *sinx, float *cosx)
+{
+    double x2 = static_cast<double>(x);
+    double s2, c2;
+    lol_sincos(x2, &s2, &c2);
+    *sinx = static_cast<float>(s2);
+    *cosx = static_cast<float>(c2);
+}
+
+double lol_tan(double x)
+{
+#if defined LOL_FEATURE_CHEAP_BRANCHES
+    double absx = lol_fabs(x * INV_PI);
+
+    /* This value was determined empirically to ensure an error of no
+     * more than x * 1e-11 in this range. */
+    if (absx < 0.163)
+    {
+        double x2 = absx * absx;
+        double x4 = x2 * x2;
+        double sub1 = (((TC[7] * x4 + TC[5]) * x4
+                           + TC[3]) * x4 + TC[1]) * x4 + ONE;
+        double sub2 = (((TC[8] * x4 + TC[6]) * x4
+                           + TC[4]) * x4 + TC[2]) * x4 + TC[0];
+        double taylor = sub2 * x2 + sub1;
+        return x * taylor;
+    }
+#endif
+
+    double sinx, cosx;
+    lol_sincos(x, &sinx, &cosx);
+
+    /* Ensure cosx isn't zero. FIXME: we lose the cosx sign here. */
+    double absc = lol_fabs(cosx);
+#if defined __CELLOS_LV2__
+    double is_cos_not_zero = absc - VERY_SMALL_NUMBER;
+    cosx = lol_fsel(is_cos_not_zero, cosx, VERY_SMALL_NUMBER);
+    return __fdiv(sinx, cosx);
+#else
+    if (__unlikely(absc < VERY_SMALL_NUMBER))
+        cosx = VERY_SMALL_NUMBER;
+    return sinx / cosx;
+#endif
+}
+
 } /* namespace lol */
 

test/lol-bench.cpp

@@ -74,12 +74,13 @@ int main(int argc, char **argv)
 
 static void bench_trig(int mode)
 {
-    float result[7] = { 0.0f };
+    float result[12] = { 0.0f };
     Timer timer;
 
     /* Set up tables */
     float *pf = new float[TRIG_TABLE_SIZE];
     float *pf2 = new float[TRIG_TABLE_SIZE];
+    float *pf3 = new float[TRIG_TABLE_SIZE];
 
     for (size_t run = 0; run < TRIG_RUNS; run++)
     {
@@ -143,27 +144,78 @@ static void bench_trig(int mode)
             pf2[i] = lol_cos(pf[i]);
         result[5] += timer.GetMs();
 
+        /* Sin & cos */
+        timer.GetMs();
+        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
+        {
+            pf2[i] = __builtin_sinf(pf[i]);
+            pf3[i] = __builtin_cosf(pf[i]);
+        }
+        result[6] += timer.GetMs();
+
+        /* Fast sin & cos */
+        timer.GetMs();
+        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
+        {
+#if defined HAVE_FASTMATH_H
+            pf2[i] = f_sinf(pf[i]);
+            pf3[i] = f_cosf(pf[i]);
+#else
+            pf2[i] = sinf(pf[i]);
+            pf3[i] = cosf(pf[i]);
+#endif
+        }
+        result[7] += timer.GetMs();
+
+        /* Lol sincos */
+        timer.GetMs();
+        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
+            lol_sincos(pf[i], &pf2[i], &pf3[i]);
+        result[8] += timer.GetMs();
+
         /* Tan */
         timer.GetMs();
         for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
             pf2[i] = __builtin_tanf(pf[i]);
-        result[6] += timer.GetMs();
+        result[9] += timer.GetMs();
+
+        /* Fast tan */
+        timer.GetMs();
+        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
+#if defined HAVE_FASTMATH_H
+            pf2[i] = f_tanf(pf[i]);
+#else
+            pf2[i] = tanf(pf[i]);
+#endif
+        result[10] += timer.GetMs();
+
+        /* Lol tan */
+        timer.GetMs();
+        for (size_t i = 0; i < TRIG_TABLE_SIZE; i++)
+            pf2[i] = lol_tan(pf[i]);
+        result[11] += timer.GetMs();
     }
 
     delete[] pf;
     delete[] pf2;
+    delete[] pf3;
 
     for (size_t i = 0; i < sizeof(result) / sizeof(*result); i++)
         result[i] *= 1000000.0f / (TRIG_TABLE_SIZE * TRIG_RUNS);
 
-    Log::Info("                          ns/elem\n");
-    Log::Info("float = sinf(float)      %7.3f\n", result[0]);
-    Log::Info("float = fastsinf(float)  %7.3f\n", result[1]);
-    Log::Info("float = lol_sinf(float)  %7.3f\n", result[2]);
-    Log::Info("float = cosf(float)      %7.3f\n", result[3]);
-    Log::Info("float = fastcosf(float)  %7.3f\n", result[4]);
-    Log::Info("float = lol_cosf(float)  %7.3f\n", result[5]);
-    Log::Info("float = tanf(float)      %7.3f\n", result[6]);
+    Log::Info("                              ns/elem\n");
+    Log::Info("float = sinf(float)          %7.3f\n", result[0]);
+    Log::Info("float = f_sinf(float)        %7.3f\n", result[1]);
+    Log::Info("float = lol_sin(float)       %7.3f\n", result[2]);
+    Log::Info("float = cosf(float)          %7.3f\n", result[3]);
+    Log::Info("float = f_cosf(float)        %7.3f\n", result[4]);
+    Log::Info("float = lol_cos(float)       %7.3f\n", result[5]);
+    Log::Info("float = sinf,cosf(float)     %7.3f\n", result[6]);
+    Log::Info("float = f_sinf,f_cosf(float) %7.3f\n", result[7]);
+    Log::Info("float = lol_sincos(float)    %7.3f\n", result[8]);
+    Log::Info("float = tanf(float)          %7.3f\n", result[9]);
+    Log::Info("float = f_tanf(float)        %7.3f\n", result[10]);
+    Log::Info("float = lol_tanf(float)      %7.3f\n", result[11]);
 }
 
 static void bench_matrix(int mode)

test/trig.cpp

@@ -70,6 +70,44 @@ public:
             double b = lol_cos(f);
             CPPUNIT_ASSERT(fabs(a - b) <= fabs(f) * 1e-11);
         }
+
+        for (int i = -10000; i < 10000; i++)
+        {
+            double f = (double)i * (1.0 / 1000.0);
+            double a1 = __builtin_sin(f);
+            double a2 = __builtin_cos(f);
+            double b1, b2;
+            lol_sincos(f, &b1, &b2);
+            CPPUNIT_ASSERT(fabs(a1 - b1) <= fabs(f) * 1e-11);
+            CPPUNIT_ASSERT(fabs(a2 - b2) <= fabs(f) * 1e-11);
+        }
+
+        for (int i = -10000; i < 10000; i++)
+        {
+            double f = (double)i * (1.0 / 100000.0);
+            double a1 = __builtin_sin(f);
+            double a2 = __builtin_cos(f);
+            double b1, b2;
+            lol_sincos(f, &b1, &b2);
+            CPPUNIT_ASSERT(fabs(a1 - b1) <= fabs(f) * 1e-11);
+            CPPUNIT_ASSERT(fabs(a2 - b2) <= fabs(f) * 1e-11);
+        }
+
+        for (int i = -10000; i < 10000; i++)
+        {
+            double f = (double)i * (1.0 / 1000.0);
+            double a = __builtin_tan(f);
+            double b = lol_tan(f);
+            CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
+        }
+
+        for (int i = -10000; i < 10000; i++)
+        {
+            double f = (double)i * (1.0 / 100000.0);
+            double a = __builtin_tan(f);
+            double b = lol_tan(f);
+            CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
+        }
     }
 };