Added support for sine and cosine

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;
+}

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);
+}

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
+}

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 */

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
+}