From 72ef0eebc70202ae2b4cf1b8b8ec4100d7a8ba7e Mon Sep 17 00:00:00 2001
From: Cody <co.dy.0.i.s.0.he.re@gmail.com>
Date: Sun, 14 Dec 2025 14:14:07 -0500
Subject: [PATCH] Added support for sine and cosine

---
 lib/libm/src/__cosdf.c | 16 ++++++++++++++++
 lib/libm/src/__sindf.c | 16 ++++++++++++++++
 lib/libm/src/cosf.c    | 31 +++++++++++++++++++++++++++++++
 lib/libm/src/libmath.h |  7 +++++++
 lib/libm/src/sinf.c    | 31 +++++++++++++++++++++++++++++++
 5 files changed, 101 insertions(+)
 create mode 100644 lib/libm/src/__cosdf.c
 create mode 100644 lib/libm/src/__sindf.c
 create mode 100644 lib/libm/src/cosf.c
 create mode 100644 lib/libm/src/sinf.c

diff --git a/lib/libm/src/__cosdf.c b/lib/libm/src/__cosdf.c
new file mode 100644
index 0000000..1aacc51
--- /dev/null
+++ b/lib/libm/src/__cosdf.c
@@ -0,0 +1,16 @@
+#include "libmath.h"
+
+static const double
+C2 = 0x1.ffff5cp-2f,
+C4 = 0x1.5568aap-5f,
+C6 = 0x1.637e8cp-10f;
+
+float __cosdf(float y) {
+    float y2 = y*y;
+    float y4 = y2*y2;
+
+    float s1 = 1 - C2*y2;
+    float s2 = C4 - C6*y2;
+
+    return s1 + s2*y4;
+}
diff --git a/lib/libm/src/__sindf.c b/lib/libm/src/__sindf.c
new file mode 100644
index 0000000..2bd6d8c
--- /dev/null
+++ b/lib/libm/src/__sindf.c
@@ -0,0 +1,16 @@
+#include "libmath.h"
+
+static const double
+C3 = 0x1.5554fep-3f,
+C5 = 0x1.10f2dcp-7f,
+C7 = 0x1.95b32ep-13f;
+
+float __sindf(float y) {
+    float y2 = y*y;
+    float y4 = y2*y2;
+
+    float s1 = 1 - C3*y2;
+    float s2 = C5 - C7*y2;
+
+    return y * (s1 + s2*y4);
+}
diff --git a/lib/libm/src/cosf.c b/lib/libm/src/cosf.c
new file mode 100644
index 0000000..2843ab3
--- /dev/null
+++ b/lib/libm/src/cosf.c
@@ -0,0 +1,31 @@
+#include "libmath.h"
+
+float cosf(float x) {
+    float y = x;
+    union { float f; int i } u = { y };
+    u.i &= 0x7fffffff; //abs
+    if(u.f <= 0.0001) { return x; } // Tiny can get away with it
+    if(u.f >= 0x7f800000) { return x-x; } // NaN / Inf
+
+    if(u.f <= 8388607.5) { // Largest representable non-integer in f32
+        // Wrap value
+        const float MAGIC = 12582912.0; // 2^23 (big enough to exploit rounding errors)
+        float xf = y * 0x1.45f306p-1f; // 2/pi
+        union { float f; int i } u = { xf };
+        int sign = 0x80000000 & u.i;
+        u.i &= 0x7fffffff; //  abs
+        u.f = u.f + MAGIC;
+        int quad = u.i & 3; // Quadrant of value
+        u.f = u.f - MAGIC; // I love IEEE-754 (rounding)
+        u.i |= sign; // copy sign
+        y -= u.f * 0x1.921fb6p0f; // pi/2
+        
+        // Evaluate
+        switch(quad) {
+            case 0: return __cosdf(y);
+            case 1: return -__sindf(y);
+            case 2: return -__cosdf(y);
+            default: return __sindf(y);
+        }
+    } else { return 0.0; } //Too much
+}
diff --git a/lib/libm/src/libmath.h b/lib/libm/src/libmath.h
index 24d894d..56d0775 100644
--- a/lib/libm/src/libmath.h
+++ b/lib/libm/src/libmath.h
@@ -9,4 +9,11 @@ double exp(double x);
 double pow(double x, double y);
 double fabs(double x);
 
+float sinf(float x);
+float cosf(float x);
+
+__attribute__((visibility("hidden")))
+float __sindf(float x);
+float __cosdf(float x);
+
 #endif /* MATH_H */
diff --git a/lib/libm/src/sinf.c b/lib/libm/src/sinf.c
new file mode 100644
index 0000000..518f329
--- /dev/null
+++ b/lib/libm/src/sinf.c
@@ -0,0 +1,31 @@
+#include "libmath.h"
+
+float sinf(float x) {
+    float y = x;
+    union { float f; int i } u = { y };
+    u.i &= 0x7fffffff; //abs
+    if(u.f <= 0.0001) { return x; } // Tiny can get away with it
+    if(u.f >= 0x7f800000) { return x-x; } // NaN / Inf
+
+    if(u.f <= 8388607.5) { // Largest representable non-integer in f32
+        // Wrap value
+        const float MAGIC = 12582912.0; // 2^23 (big enough to exploit rounding errors)
+        float xf = y * 0x1.45f306p-1f; // 2/pi
+        union { float f; int i } u = { xf };
+        int sign = 0x80000000 & u.i;
+        u.i &= 0x7fffffff; //  abs
+        u.f = u.f + MAGIC;
+        int quad = u.i & 3; // Quadrant of value
+        u.f = u.f - MAGIC; // I love IEEE-754 (rounding)
+        u.i |= sign; // copy sign
+        y -= u.f * 0x1.921fb6p0f; // pi/2
+        
+        // Evaluate
+        switch(quad) {
+            case 0: return __sindf(y);
+            case 1: return __cosdf(y);
+            case 2: return -__sindf(y);
+            default: return -__cosdf(y);
+        }
+    } else { return 0.0; } //Too much
+}