implement horners rule

src/math_func_checker.c

@@ -5,9 +5,15 @@
 
 typedef double f64;
 typedef float f32;
+typedef uint64_t u64;
 typedef f64 (*math_func)(f64);
 #define ARRAY_LEN(x) (sizeof(x)/sizeof(x[0]))
 
+struct tester {
+    f64 max_error;
+    f64 max_error_at;
+};
+
 struct math_spec {
     const char *name;
     math_func mine;
@@ -80,6 +86,11 @@ static f64 my_sin_quarter_fitted(f64 num)
     return result * sign;
 }
 
+static f64 my_cos_using_sin(f64 num)
+{
+    return my_sin_quarter_fitted(num + M_PI/2);
+}
+
 static f64 my_cos(f64 num)
 {
     if (num < 0) {
@@ -113,9 +124,107 @@ static f64 my_asin(f64 num)
     return num;
 }
 
+static f64 factorial(size_t n)
+{
+    f64 result = 1;
+    for (size_t i = 1; i <= n; i++) {
+        result *= i;
+    }
+    return result;
+}
+
+static void create_taylor_series_coeffs(f64 *coeffs, size_t coeffs_len)
+{
+    for (size_t i = 0; i < coeffs_len; i++) {
+        size_t n = 1 + (i*2);
+        coeffs[i] = 1 / factorial(n);
+        if (i % 2 == 1) {
+            coeffs[i] *= -1;
+        }
+    }
+}
+
+static f64 sin_using_taylor(f64 input, f64 *coeffs, size_t coeffs_len)
+{
+    f64 result = 0;
+    for (size_t i = 0; i < coeffs_len; i++) {
+        f64 power = 1 + (i*2);
+        result += coeffs[i] * powf(input, power);
+    }
+    return result;
+}
+
+static f64 sin_using_taylor_horner(f64 input, f64 *coeffs, size_t coeffs_len)
+{
+    f64 input_sqr = input * input;
+
+    f64 result = coeffs[coeffs_len - 1];
+    for (size_t j = 1; j < coeffs_len; j++) {
+        size_t i = coeffs_len - j - 1;
+
+        result *= input_sqr;
+        result += coeffs[i];
+    }
+
+    result *= input;
+
+    return result;
+}
+
+static void tester_record_result(struct tester *tester, f64 input, f64 my_output, f64 reference_output)
+{
+    f64 error = reference_output - my_output;
+    if (error < 0) {
+        error = -error;
+    }
+
+    if (error > tester->max_error) {
+        tester->max_error = error;
+        tester->max_error_at = input;
+    }
+}
+
+static void tester_show_result(struct tester *tester, const char *name)
+{
+    printf("%20s max error %.16f at %f\n", name, tester->max_error, tester->max_error_at);
+}
+
+static f64 lerp(f64 x, f64 from, f64 to)
+{
+    return (to - from) * x + from;
+}
+
+static void test_taylor_series(u64 samples, f64 *coeffs, size_t coeffs_len)
+{
+    create_taylor_series_coeffs(coeffs, coeffs_len);
+
+    struct tester tester = { };
+    struct tester tester_horner = { };
+
+    for (int j = 0; j < samples; j++) {
+        f64 input = lerp((f64)j / samples, 0, M_PI/2);
+        tester_record_result(&tester, input, sin_using_taylor(input, coeffs, coeffs_len), sin(input));
+        tester_record_result(&tester_horner, input, sin_using_taylor_horner(input, coeffs, coeffs_len), sin(input));
+    }
+
+    {
+        char name[64] = { 0 };
+        snprintf(name, sizeof(name), "taylor   %zu", coeffs_len);
+        tester_show_result(&tester, name);
+    }
+
+    {
+        char name[64] = { 0 };
+        snprintf(name, sizeof(name), "taylor h %zu", coeffs_len);
+        tester_show_result(&tester_horner, name);
+    }
+
+    printf("            diff               %.16f\n", tester.max_error - tester_horner.max_error);
+}
+
 int main()
 {
-    uint64_t samples = 10000;
+    uint64_t samples = 100000;
     struct math_spec specs[] = {
         {
             .name = "sin",
@@ -138,6 +247,13 @@ int main()
             .from = -M_PI,
             .to = M_PI
         },
+        {
+            .name = "my_cos_using_sin",
+            .mine = my_cos_using_sin,
+            .reference = cos,
+            .from = -M_PI,
+            .to = M_PI
+        },
         {
             .name = "cos",
             .mine = my_cos,
@@ -171,26 +287,19 @@ int main()
     for (int i = 0; i < ARRAY_LEN(specs); i++) {
         struct math_spec spec = specs[i];
 
-        f64 max_error = 0;
-        f64 max_error_at = 0;
+        struct tester tester = { };
 
         for (int j = 0; j < samples; j++) {
-            f64 input = ((f64)j / samples) * (spec.to - spec.from) + spec.from;
-            f64 reference_output = spec.reference(input);
-            f64 my_output = spec.mine(input);
-
-            f64 error = reference_output - my_output;
-            if (error < 0) {
-                error = -error;
-            }
-
-            if (error > max_error) {
-                max_error = error;
-                max_error_at = input;
-            }
+            f64 input = lerp((f64)j / samples, spec.from, spec.to);
+            tester_record_result(&tester, input, spec.mine(input), spec.reference(input));
         }
 
-        printf("%20s max error %f at %f\n", spec.name, max_error, max_error_at);
+        tester_show_result(&tester, spec.name);
+    }
+
+    f64 taylor_coeffs[10];
+    for (int i = 1; i < ARRAY_LEN(taylor_coeffs); i++) {
+        test_taylor_series(samples, taylor_coeffs, i);
     }
 
     return 0;