LibM: Implement various trig functions

Patch from Anonymous.

Libraries/LibM/TestMath.cpp

@@ -0,0 +1,54 @@
+#include <AK/TestSuite.h>
+ 
+#include <math.h>
+
+#define EXPECT_CLOSE(a, b) { EXPECT(fabs(a - b) < 0.000001); }
+
+TEST_CASE(trig)
+{
+    EXPECT_CLOSE(sin(1234), 0.653316);
+    EXPECT_CLOSE(cos(1234), -0.830914);
+    EXPECT_CLOSE(tan(1234), -0.786262);
+    EXPECT_CLOSE(sqrt(1234), 35.128336)
+    EXPECT_CLOSE(sin(-1), -0.867955);
+    EXPECT_CLOSE(cos(-1), 0.594715);
+    EXPECT_CLOSE(tan(-1), -1.459446);
+    EXPECT(isnan(sqrt(-1)));
+}
+
+TEST_CASE(other)
+{
+    EXPECT_EQ(trunc(9999999999999.5), 9999999999999.0);
+    EXPECT_EQ(trunc(-9999999999999.5), -9999999999999.0);
+}
+
+TEST_CASE(exponents)
+{
+    struct values {
+        double x;
+        double exp;
+        double sinh;
+        double cosh;
+        double tanh;
+    };
+
+    values values[8] {
+        { 1.500000, 4.481626, 2.129246, 2.352379, 0.905148},
+        { 20.990000, 1304956710.432035, 652478355.216017, 652478355.216017, 1.000000},
+        { 20.010000, 490041186.687082, 245020593.343541, 245020593.343541, 1.000000},
+        { 0.000000, 1.000000, 0.000000, 1.000000, 0.000000},
+        { 0.010000, 1.010050, 0.010000, 1.000050, 0.010000},
+        { -0.010000, 0.990050, -0.010000, 1.000050, -0.010000},
+        { -1.000000, 0.367879, -1.175201, 1.543081, -0.761594},
+        { -17.000000, 0.000000, -12077476.376788, 12077476.376788, -1.000000},
+    };
+    for (auto& v : values) {
+        EXPECT_CLOSE(exp(v.x), v.exp);
+        EXPECT_CLOSE(sinh(v.x), v.sinh);
+        EXPECT_CLOSE(cosh(v.x), v.cosh);
+        EXPECT_CLOSE(tanh(v.x), v.tanh);
+    }
+    EXPECT_EQ(exp(1000), std::numeric_limits<double>::infinity());
+}
+
+TEST_MAIN(Math)

Libraries/LibM/math.cpp

@@ -1,10 +1,21 @@
 #include <LibC/assert.h>
 #include <LibM/math.h>
+#include <limits>
+#include <stdint.h>
+#include <stdlib.h>
+
+template<size_t> constexpr double e_to_power();
+template<> constexpr double e_to_power<0>() { return 1; }
+template<size_t exponent> constexpr double e_to_power() { return M_E * e_to_power<exponent - 1>(); }
+
+template<size_t> constexpr size_t factorial();
+template<> constexpr size_t factorial<0>() { return 1; }
+template<size_t value> constexpr size_t factorial() { return value * factorial<value - 1>(); }
 
 extern "C" {
 double trunc(double x)
 {
-    return (int)x;
+    return (int64_t)x;
 }
 
 double cos(double angle)
@@ -40,17 +51,25 @@ double pow(double x, double y)
     (void)x;
     (void)y;
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double ldexp(double, int exp)
 {
     (void)exp;
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
-double tanh(double)
+double tanh(double x)
 {
-    ASSERT_NOT_REACHED();
+    if (x > 0) {
+        double exponentiated = exp(2 * x);
+        return (exponentiated - 1) / (exponentiated + 1);
+    }
+    double plusX = exp(x);
+    double minusX = exp(-x);
+    return (plusX - minusX) / (plusX + minusX);
 }
 
 double tan(double angle)
@@ -65,19 +84,25 @@ double sqrt(double x)
     return res;
 }
 
-double sinh(double)
+double sinh(double x)
 {
-    ASSERT_NOT_REACHED();
+    if (x > 0) {
+        double exponentiated = exp(x);
+        return (exponentiated * exponentiated - 1) / 2 / exponentiated;
+    }
+    return (exp(x) - exp(-x)) / 2;
 }
 
 double log10(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double log(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double fmod(double index, double period)
@@ -85,67 +110,107 @@ double fmod(double index, double period)
     return index - trunc(index / period) * period;
 }
 
-double exp(double)
-{
-    ASSERT_NOT_REACHED();
-}
-
-double cosh(double)
-{
-    ASSERT_NOT_REACHED();
+double exp(double exponent)
+{
+    double result = 1;
+    if (exponent >= 1) {
+        size_t integer_part = (size_t)exponent;
+        if (integer_part & 1) result *= e_to_power<1>();
+        if (integer_part & 2) result *= e_to_power<2>();
+        if (integer_part > 3) {
+            if (integer_part & 4) result *= e_to_power<4>();
+            if (integer_part & 8) result *= e_to_power<8>();
+            if (integer_part & 16) result *= e_to_power<16>();
+            if (integer_part & 32) result *= e_to_power<32>();
+            if (integer_part >= 64) return std::numeric_limits<double>::infinity();
+        }
+        exponent -= integer_part;
+    } else if (exponent < 0)
+        return 1 / exp(-exponent);
+    double taylor_series_result = 1 + exponent;
+    double taylor_series_numerator = exponent * exponent;
+    taylor_series_result += taylor_series_numerator / factorial<2>();
+    taylor_series_numerator *= exponent;
+    taylor_series_result += taylor_series_numerator / factorial<3>();
+    taylor_series_numerator *= exponent;
+    taylor_series_result += taylor_series_numerator / factorial<4>();
+    taylor_series_numerator *= exponent;
+    taylor_series_result += taylor_series_numerator / factorial<5>();
+    return result * taylor_series_result;
+}
+
+double cosh(double x)
+{
+    if (x < 0) {
+        double exponentiated = exp(-x);
+        return (1 + exponentiated * exponentiated) / 2 / exponentiated;
+    }
+    return (exp(x) + exp(-x)) / 2;
 }
 
 double atan2(double, double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double atan(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double asin(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double acos(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double fabs(double value)
 {
     return value < 0 ? -value : value;
 }
+
 double log2(double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 float log2f(float)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 long double log2l(long double)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 double frexp(double, int*)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 float frexpf(float, int*)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
 
 long double frexpl(long double, int*)
 {
     ASSERT_NOT_REACHED();
+    return 0;
 }
+
 }

Libraries/LibM/math.h

@@ -5,6 +5,7 @@
 __BEGIN_DECLS
 
 #define HUGE_VAL 1e10000
+#define M_E 2.718281828459045
 #define M_PI 3.141592653589793
 #define M_PI_2 (M_PI / 2)
 #define M_TAU (M_PI * 2)