LibM: Implement various functions.

Path from Anonymous.

3 files changed, 146 insertions, 52 deletions

diff --git a/Libraries/LibM/TestMath.cpp b/Libraries/LibM/TestMath.cpp
index 761ccb3150..48dc211267 100644
--- a/Libraries/LibM/TestMath.cpp
+++ b/Libraries/LibM/TestMath.cpp
@@ -14,6 +14,24 @@ TEST_CASE(trig)
     EXPECT_CLOSE(cos(-1), 0.594715);
     EXPECT_CLOSE(tan(-1), -1.459446);
     EXPECT(isnan(sqrt(-1)));
+    EXPECT(isnan(asin(1.1)));
+    EXPECT(isnan(asin(-1.1)));
+    EXPECT_CLOSE(asin(0), 0.0);
+    EXPECT_CLOSE(asin(0.01), 0.01);
+    EXPECT_CLOSE(asin(0.1), 0.100167);
+    EXPECT_CLOSE(asin(0.3), 0.304693);
+    EXPECT_CLOSE(asin(0.499), 0.522444);
+    EXPECT_CLOSE(asin(0.5), 0.523599);
+    EXPECT_CLOSE(asin(0.501), 0.524754);
+    EXPECT_CLOSE(asin(0.9), 1.119770);
+    EXPECT_CLOSE(asin(0.99), 1.429246);
+    EXPECT_CLOSE(asin(1.0), 1.570750);
+    EXPECT_CLOSE(atan(0), 0.0)
+    EXPECT_CLOSE(atan(0.5), 0.463648)
+    EXPECT_CLOSE(atan(-0.5), -0.463648)
+    EXPECT_CLOSE(atan(5.5), 1.390943)
+    EXPECT_CLOSE(atan(-5.5), -1.390943)
+    EXPECT_CLOSE(atan(555.5), 1.568996)
 }
 
 TEST_CASE(other)
@@ -51,4 +69,17 @@ TEST_CASE(exponents)
     EXPECT_EQ(exp(1000), std::numeric_limits<double>::infinity());
 }
 
+TEST_CASE(logarithms)
+{
+    EXPECT(isnan(log(-1)));
+    EXPECT(log(0) < -1000000);
+    EXPECT_CLOSE(log(0.5), -0.693233)
+    EXPECT_CLOSE(log(1.1), 0.095310)
+    EXPECT_CLOSE(log(5), 1.609480)
+    EXPECT_CLOSE(log(5.5), 1.704842)
+    EXPECT_CLOSE(log(500), 6.214104)
+    EXPECT_CLOSE(log2(5), 2.321989)
+    EXPECT_CLOSE(log10(5), 0.698988)
+}
+
 TEST_MAIN(Math)
diff --git a/Libraries/LibM/math.cpp b/Libraries/LibM/math.cpp
index a45e4fcc7d..0443a01f75 100644
--- a/Libraries/LibM/math.cpp
+++ b/Libraries/LibM/math.cpp
@@ -3,13 +3,33 @@
 #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>(); }
+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>(); }
+
+template<size_t>
+constexpr size_t product_even();
+template<>
+constexpr size_t product_even<2>() { return 2; }
+template<size_t value>
+constexpr size_t product_even() { return value * product_even<value - 2>(); }
+
+template<size_t>
+constexpr size_t product_odd();
+template<>
+constexpr size_t product_odd<1>() { return 1; }
+template<size_t value>
+constexpr size_t product_odd() { return value * product_odd<value - 2>(); }
 
 extern "C" {
 double trunc(double x)
@@ -67,7 +87,7 @@ double tanh(double x)
         return (exponentiated - 1) / (exponentiated + 1);
     }
     double plusX = exp(x);
-    double minusX = exp(-x);
+    double minusX = 1 / plusX;
     return (plusX - minusX) / (plusX + minusX);
 }
 
@@ -79,29 +99,38 @@ double tan(double angle)
 double sqrt(double x)
 {
     double res;
-    __asm__("fsqrt" : "=t"(res) : "0"(x));
+    __asm__("fsqrt"
+            : "=t"(res)
+            : "0"(x));
     return res;
 }
 
 double sinh(double x)
 {
-    if (x > 0) {
-        double exponentiated = exp(x);
+    double exponentiated = exp(x);
+    if (x > 0)
         return (exponentiated * exponentiated - 1) / 2 / exponentiated;
-    }
-    return (exp(x) - exp(-x)) / 2;
+    return (exponentiated - 1 / exponentiated) / 2;
 }
 
-double log10(double)
+double log10(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    return log(x) / M_LN10;
 }
 
-double log(double)
+double log(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    if (x < 0)
+        return __builtin_nan("");
+    if (x == 0)
+        return -__builtin_huge_val();
+    double y = 1 + 2 * (x - 1) / (x + 1);
+    double exponentiated = exp(y);
+    y = y + 2 * (x - exponentiated) / (x + exponentiated);
+    exponentiated = exp(y);
+    y = y + 2 * (x - exponentiated) / (x + exponentiated);
+    exponentiated = exp(y);
+    return y + 2 * (x - exponentiated) / (x + exponentiated);
 }
 
 double fmod(double index, double period)
@@ -114,14 +143,21 @@ 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 & 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 __builtin_huge_val();
+            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 __builtin_huge_val();
         }
         exponent -= integer_part;
     } else if (exponent < 0)
@@ -140,35 +176,64 @@ double exp(double exponent)
 
 double cosh(double x)
 {
-    if (x < 0) {
-        double exponentiated = exp(-x);
+    double exponentiated = exp(-x);
+    if (x < 0)
         return (1 + exponentiated * exponentiated) / 2 / exponentiated;
-    }
-    return (exp(x) + exp(-x)) / 2;
+    return (1 / exponentiated + exponentiated) / 2;
 }
 
-double atan2(double, double)
+double atan2(double y, double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    if (x > 0)
+        return atan(y / x);
+    if (x == 0) {
+        if (y > 0)
+            return M_PI_2;
+        if (y < 0)
+            return -M_PI_2;
+        return 0;
+    }
+    if (y >= 0)
+        return atan(y / x) + M_PI;
+    return atan(y / x) - M_PI;
 }
 
-double atan(double)
+double atan(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    if (x < 0)
+        return -atan(-x);
+    if (x > 1)
+        return M_PI_2 - atan(1 / x);
+    double squared = x * x;
+    return x / (1 + 1 * 1 * squared / (3 + 2 * 2 * squared / (5 + 3 * 3 * squared / (7 + 4 * 4 * squared / (9 + 5 * 5 * squared / (11 + 6 * 6 * squared / (13 + 7 * 7 * squared)))))));
 }
 
-double asin(double)
+double asin(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    if (x > 1 || x < -1)
+        return __builtin_nan("");
+    if (x > 0.5 || x < -0.5)
+        return 2 * atan(x / (1 + sqrt(1 - x * x)));
+    double squared = x * x;
+    double value = x;
+    double i = x * squared;
+    value += i * product_odd<1>() / product_even<2>() / 3;
+    i *= squared;
+    value += i * product_odd<3>() / product_even<4>() / 5;
+    i *= squared;
+    value += i * product_odd<5>() / product_even<6>() / 7;
+    i *= squared;
+    value += i * product_odd<7>() / product_even<8>() / 9;
+    i *= squared;
+    value += i * product_odd<9>() / product_even<10>() / 11;
+    i *= squared;
+    value += i * product_odd<11>() / product_even<12>() / 13;
+    return value;
 }
 
-double acos(double)
+double acos(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    return M_PI_2 - asin(x);
 }
 
 double fabs(double value)
@@ -176,22 +241,19 @@ double fabs(double value)
     return value < 0 ? -value : value;
 }
 
-double log2(double)
+double log2(double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    return log(x) / M_LN2;
 }
 
-float log2f(float)
+float log2f(float x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    return log2(x);
 }
 
-long double log2l(long double)
+long double log2l(long double x)
 {
-    ASSERT_NOT_REACHED();
-    return 0;
+    return log2(x);
 }
 
 double frexp(double, int*)
@@ -211,5 +273,4 @@ long double frexpl(long double, int*)
     ASSERT_NOT_REACHED();
     return 0;
 }
-
 }
diff --git a/Libraries/LibM/math.h b/Libraries/LibM/math.h
index 548d49cc29..cf0eb99cab 100644
--- a/Libraries/LibM/math.h
+++ b/Libraries/LibM/math.h
@@ -9,6 +9,8 @@ __BEGIN_DECLS
 #define M_PI 3.141592653589793
 #define M_PI_2 (M_PI / 2)
 #define M_TAU (M_PI * 2)
+#define M_LN2 0.69314718055995
+#define M_LN10 2.30258509299405
 
 double acos(double);
 float acosf(float);