From 3cab32b171eea1980d2b2069fe7e2fa2e3db24f0 Mon Sep 17 00:00:00 2001
From: Rokas Puzonas <rokas.puz@gmail.com>
Date: Thu, 30 Jan 2025 23:16:11 +0200
Subject: [PATCH] fit a 6-th degree curve to sin for coefficients

---
 fit_curve.py            | 30 ++++++++++++++
 src/math_func_checker.c | 86 +++++++++++++++++++++++++++++++++++++++--
 2 files changed, 113 insertions(+), 3 deletions(-)
 create mode 100644 fit_curve.py

diff --git a/fit_curve.py b/fit_curve.py
new file mode 100644
index 0000000..04d762d
--- /dev/null
+++ b/fit_curve.py
@@ -0,0 +1,30 @@
+import numpy as np
+from scipy.optimize import curve_fit
+from matplotlib import pyplot as plt
+import math
+
+x = np.linspace(0, math.pi / 2, num=1000)
+y = np.sin(x)
+
+
+def test(x, a, b, c, d, e, f):
+    return a * x**6 + b * x**5 + c * x**4 + d * x**3 + e * x**2 + f * x
+
+
+param, param_cov = curve_fit(test, x, y)
+
+print("Sine function coefficients:")
+print(*param)
+
+
+# ans stores the new y-data according to
+# the coefficients given by curve-fit() function
+ans = test(x, *param)
+
+print("Max error:")
+print(abs(y-ans).max())
+
+plt.plot(x, y, 'o', color='red', label="data")
+plt.plot(x, ans, '--', color='blue', label="optimized data")
+plt.legend()
+plt.show()
diff --git a/src/math_func_checker.c b/src/math_func_checker.c
index 699ffc4..d46ace6 100644
--- a/src/math_func_checker.c
+++ b/src/math_func_checker.c
@@ -16,14 +16,78 @@ struct math_spec {
     f64 to;
 };
 
+static f64 square(f64 num)
+{
+    return num * num;
+}
+
 static f64 my_sin(f64 num)
 {
-    return num;
+    f64 sign = 1;
+    if (num < 0) {
+        sign = -1;
+        num = -num;
+    }
+
+    f64 result =  (4 / M_PI) * num - square((2 / M_PI) * num);
+    return result * sign;
+}
+
+static f64 my_sin_quarter(f64 num)
+{
+    f64 sign = 1;
+    if (num < 0) {
+        sign *= -1;
+        num = -num;
+    }
+
+    if (num > M_PI / 2) {
+        num = M_PI - num;
+    }
+
+    f64 result =  (4 / M_PI) * num - square((2 / M_PI) * num);
+    return result * sign;
+}
+
+static f64 my_sin_quarter_fitted(f64 num)
+{
+    f64 sign = 1;
+    if (num < 0) {
+        sign *= -1;
+        num = -num;
+    }
+
+    if (num > M_PI / 2) {
+        num = M_PI - num;
+    }
+
+    f64 coeffs[] = {
+        -0.0009462013574375033,
+        0.010206999409602957,
+        -0.0018834448707133052,
+        -0.16568577654109595,
+        -0.00024161815691126504,
+        1.0000207681130697
+    };
+
+    f64 result = 0;
+    result += coeffs[0] * num * num * num * num * num * num;
+    result += coeffs[1] * num * num * num * num * num;
+    result += coeffs[2] * num * num * num * num;
+    result += coeffs[3] * num * num * num;
+    result += coeffs[4] * num * num;
+    result += coeffs[5] * num;
+    return result * sign;
 }
 
 static f64 my_cos(f64 num)
 {
-    return num;
+    if (num < 0) {
+        num *= -1;
+    }
+
+    f64 result = 1 - square((2 / M_PI) * num);
+    return result;
 }
 
 static f64 my_sqrt1(f64 num)
@@ -60,6 +124,20 @@ int main()
             .from = -M_PI,
             .to = M_PI
         },
+        {
+            .name = "sin_quater",
+            .mine = my_sin_quarter,
+            .reference = sin,
+            .from = -M_PI,
+            .to = M_PI
+        },
+        {
+            .name = "sin_quater_fit",
+            .mine = my_sin_quarter_fitted,
+            .reference = sin,
+            .from = -M_PI,
+            .to = M_PI
+        },
         {
             .name = "cos",
             .mine = my_cos,
@@ -94,6 +172,7 @@ int main()
         struct math_spec spec = specs[i];
 
         f64 max_error = 0;
+        f64 max_error_at = 0;
 
         for (int j = 0; j < samples; j++) {
             f64 input = ((f64)j / samples) * (spec.to - spec.from) + spec.from;
@@ -107,10 +186,11 @@ int main()
 
             if (error > max_error) {
                 max_error = error;
+                max_error_at = input;
             }
         }
 
-        printf("%s max error %f\n", spec.name, max_error);
+        printf("%20s max error %f at %f\n", spec.name, max_error, max_error_at);
     }
 
     return 0;