diff --git a/src/trig.cpp b/src/trig.cpp index d587f4d7..1ea92291 100644 --- a/src/trig.cpp +++ b/src/trig.cpp @@ -24,6 +24,7 @@ namespace lol { static const double PI = 3.14159265358979323846; +static const double NEG_PI = -3.14159265358979323846; static const double PI_2 = PI / 2.0; static const double PI_4 = PI / 4.0; static const double INV_PI = 1.0 / PI; @@ -35,20 +36,20 @@ static const double NEG_ONE = -1.0; static const double HALF = 0.5; static const double TWO = 2.0; static const double VERY_SMALL_NUMBER = 0x1.0p-128; -static const double VERY_LARGE_NUMBER = 4503599627370496.0; -static const double EVEN_LARGER_NUMBER = 9007199254740992.0; +static const double TWO_EXP_52 = 4503599627370496.0; +static const double TWO_EXP_54 = 18014398509481984.0; /** sin Taylor series coefficients. */ static const double SC[] = { - -1.6666666666666666666667e-01, // 1/3! - +8.3333333333333333333333e-03, // 1/5! - -1.9841269841269841269841e-04, // 1/7! - +2.7557319223985890652557e-06, // 1/9! - -2.5052108385441718775052e-08, // 1/11! - +1.6059043836821614599392e-10, // 1/13! - -7.6471637318198164759011e-13, // 1/15! - +2.8114572543455207631989e-15, // 1/17! + -1.6449340668482264364724e-0, // pi^2/3! + +8.1174242528335364363700e-1, // pi^4/5! + -1.9075182412208421369647e-1, // pi^6/7! + +2.6147847817654800504653e-2, // pi^8/9! + -2.3460810354558236375089e-3, // pi^10/11! + +1.4842879303107100368487e-4, // pi^12/13! + -6.9758736616563804745344e-6, // pi^14/15! + +2.5312174041370276513517e-7, // pi^16/17! }; /* Custom intrinsics */ @@ -165,32 +166,30 @@ static inline double lol_trunc(double x) double lol_sin(double x) { - double absx = lol_fabs(x); - double sign = lol_fsel(x, ONE, NEG_ONE); + double absx = lol_fabs(x * INV_PI); + double sign = lol_fsel(x, PI, NEG_PI); #if defined __CELLOS_LV2__ - double num_cycles = lol_round(absx * INV_PI); + double num_cycles = lol_round(absx); double is_even = lol_trunc(num_cycles * HALF) - (num_cycles * HALF); + sign = lol_fsel(is_even, sign, -sign); #else - double num_cycles = absx * INV_PI + VERY_LARGE_NUMBER; - __asm__("" : "+m" (num_cycles)); - num_cycles -= VERY_LARGE_NUMBER; + double num_cycles = absx + TWO_EXP_52; + __asm__("" : "+m" (num_cycles)); num_cycles -= TWO_EXP_52; - double is_even = num_cycles - HALF; - __asm__("" : "+m" (is_even)); - is_even += EVEN_LARGER_NUMBER; + 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 -= EVEN_LARGER_NUMBER; - __asm__("" : "+m" (is_even)); - is_even -= num_cycles; + is_even -= TWO * num_cycles - ONE; + sign *= is_even; #endif - double norm_x = absx - PI * num_cycles; + double norm_x = absx - num_cycles; double y = norm_x * norm_x; - double taylor = ((((((((SC[7] * y + SC[6]) * y + SC[5]) - * y + SC[4]) * y + SC[3]) - * y + SC[2]) * y + SC[1]) - * y + SC[0]) * y); - sign = lol_fsel(is_even, sign, -sign); - double result = norm_x + norm_x * taylor; + double taylor = (((((((SC[7] * y + SC[6]) * y + SC[5]) + * y + SC[4]) * y + SC[3]) + * y + SC[2]) * y + SC[1]) + * y + SC[0]) * y + ONE; + double result = norm_x * taylor; return result * sign; } diff --git a/test/lol-bench.cpp b/test/lol-bench.cpp index 392aa96b..eb8c813c 100644 --- a/test/lol-bench.cpp +++ b/test/lol-bench.cpp @@ -49,6 +49,11 @@ int main(int argc, char **argv) Log::Info("------------------------\n"); bench_trig(2); + Log::Info("----------------------------\n"); + Log::Info(" Trigonometry [-1e-2, 1e-2]\n"); + Log::Info("----------------------------\n"); + bench_trig(3); + Log::Info("----------------------------\n"); Log::Info(" Float matrices [-2.0, 2.0]\n"); Log::Info("----------------------------\n"); @@ -88,6 +93,10 @@ static void bench_trig(int mode) for (size_t i = 0; i < TRIG_TABLE_SIZE; i++) pf[i] = RandF(-M_PI, M_PI); break; + case 3: + for (size_t i = 0; i < TRIG_TABLE_SIZE; i++) + pf[i] = RandF(-1e-2f, 1e-2f); + break; } /* Sin */