Add benchmark coverage for parabola based cosine approximation

CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 3.12)
+cmake_policy(VERSION 3.12)
+
+project(cos_benchmarks
+    LANGUAGES C
+    VERSION 0.1
+    DESCRIPTION "cos implementation benchmarks"
+)
+
+add_executable(benchmarks benchmarks.c)
+target_link_libraries(benchmarks PRIVATE m)

benchmarks.c

@@ -17,12 +17,11 @@ double runtime(double (*func)(double))
 }
 
 // Finds the worst case for accuracy compared to math.h. Smaller number is better.
-double accuracy(double (*func)(double))
+double accuracy(double (*func)(double), double start, double stop)
 {
     double w = -1;
-    double start = 0;
-    double stop = CONST_2PI;
-    double step = 0.0000001;
+    const double steps = (2. * CONST_PI) / 0.0000001;
+    const double step = (stop - start) / steps;
     for (double i = start; i < stop; i += step)
     {
         double c = absd(func(i) - cos(i));
@@ -66,6 +65,11 @@ static struct {
     RTEST(cos_table_0_0001_LERP),
     /* Built-in cosine */
     RTEST(cos_math_h),
+    /* parabola based approximation */
+    RTEST(cos_parabola),
+    RTEST(cos_parabola_extra),
+    RTEST(cos_parabola_opt),
+    RTEST(cos_parabola_extra_opt),
 };
 
 const int num_tests = sizeof(tests) / sizeof(*tests);
@@ -75,6 +79,8 @@ int main(int argc, char *argv[])
 {
     int run_accuracy = 1;
     int run_runtime = 1;
+    double start = 0;
+    double stop = CONST_2PI;
     int i;
     int j;
 
@@ -94,7 +100,7 @@ int main(int argc, char *argv[])
             for (j = 0;j < num_tests;j++) {
                 if (!strcmp(argv[i], tests[j].name)) {
                     printf("ACCURACY\n");
-                    printf("%-35s %.16lf\n", tests[j].name, accuracy(tests[j].func));
+                    printf("%-35s %.16lf\n", tests[j].name, accuracy(tests[j].func, start, stop));
                     printf("\nTIME\n");            
                     printf("%-35s %.16lf\n", tests[j].name, runtime(tests[j].func));
                     return 0;
@@ -110,8 +116,16 @@ int main(int argc, char *argv[])
             }
             return 0;
         }
+        else if (!strcmp(argv[i], "-r")) {
+            start = -CONST_PI;
+            stop = CONST_PI;
+        }
+        else if (!strcmp(argv[i], "-R")) {
+            start = -CONST_PI * 10;
+            stop = CONST_PI * 10;
+        }
         else {
-            printf("Usage: %s [-na] [-nt] [-t <testname>]\n   -na - Don't run accuracy tests\n   -nt - Don't run speed tests.\n   -t <testname> - Run a particular test instead of all tests.\n   -p - Print all test names.\n", argv[0]);
+            printf("Usage: %s [-na] [-nt] [-r|-R] [-t <testname>]\n   -na - Don't run accuracy tests\n   -nt - Don't run speed tests.\n   -t <testname> - Run a particular test instead of all tests.\n   -p - Print all test names.\n   -r - Run accuracy test in range [-pi,pi] instead of [0,2pi].\n   -R - Run accuracy test in range [-10pi,10pi] instead of [0,2pi].\n", argv[0]);
             return 0;
         }
     }
@@ -121,7 +135,7 @@ int main(int argc, char *argv[])
     if (run_accuracy) {
         printf("ACCURACY\n");
         for (i = 0;i < num_tests;i++) {
-            printf("%-35s %.16lf\n", tests[i].name, accuracy(tests[i].func));
+            printf("%-35s %.16lf\n", tests[i].name, accuracy(tests[i].func, start, stop));
         }
     }
 

cosine.c

@@ -310,3 +310,69 @@ double cos_table_0_0001_LERP(double x)
 double cos_math_h(double x) {
     return cos(x);
 }
+
+///
+/// Parabola based cos approximations
+///
+/// Original by Nick for sin can be found at [1]. The cos
+/// version here is deduced from `cos(x) = sin(x + pi/2)`.
+/// Additionally, the original only accounts for input values in [-pi..pi],
+/// making its precision absolutely horrible in the benchmark tests.
+///
+/// The optimized versions have been deduced by Milian Wolff, see [2].
+/// They are quite quick and don't have abysmal precision, giving them
+/// a nice position in the middle for both metrics.
+///
+/// [1]: https://web.archive.org/web/20171228230531/http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648
+/// [2]: https://stackoverflow.com/a/28050328/35250
+///
+
+double cos_parabola(double x)
+{
+    // cos(x) = sin(x + pi/2)
+    x += (CONST_PI / 2.);
+
+    const double B = 4. / CONST_PI;
+    const double C = -4. / (CONST_PI * CONST_PI);
+
+    double y = B * x + C * x * fabs(x);
+
+    return y;
+}
+
+double cos_parabola_extra(double x)
+{
+    // cos(x) = sin(x + pi/2)
+    x += (CONST_PI / 2.);
+
+    const double B = 4. / CONST_PI;
+    const double C = -4. / (CONST_PI * CONST_PI);
+
+    double y = B * x + C * x * fabs(x);
+
+    // const float Q = 0.775;
+    const double P = 0.225;
+
+    y = P * (y * fabs(y) - y) + y;   // Q * y + P * y * abs(y)
+
+    return y;
+}
+
+double cos_parabola_opt(double x)
+{
+    const double tp = 1. / CONST_2PI;
+    x *= tp;
+    x -= .25 + floor(x + .25);
+    x *= 16. * (fabs(x) - .5);
+    return x;
+}
+
+double cos_parabola_extra_opt(double x)
+{
+    const double tp = 1. / CONST_2PI;
+    x *= tp;
+    x -= .25 + floor(x + .25);
+    x *= 16. * (fabs(x) - .5);
+    x += .225 * x * (fabs(x) - 1.);
+    return x;
+}