LibM: Reimplement some functions using FPU

Libraries/LibM/math.cpp

@@ -102,8 +102,12 @@ float sinf(float angle) NOEXCEPT
 double pow(double x, double y) NOEXCEPT
 {
     // FIXME: Please fix me. I am naive.
+    if (isnan(y))
+        return y;
     if (y == 0)
         return 1;
+    if (x == 0)
+        return 0;
     if (y == 1)
         return x;
     int y_as_int = (int)y;
@@ -115,40 +119,22 @@ double pow(double x, double y) NOEXCEPT
             result = 1.0 / result;
         return result;
     }
-    return exp(y * log(x));
+    return exp2(y * log2(x));
 }
 
 float powf(float x, float y) NOEXCEPT
 {
-    // FIXME: Please fix me. I am naive.
+    return (float)pow(x, y);
-    if (y == 0)
-        return 1;
-    if (y == 1)
-        return x;
-    int y_as_int = (int)y;
-    if (y == (float)y_as_int) {
-        float result = x;
-        for (int i = 0; i < fabs(y) - 1; ++i)
-            result *= x;
-        if (y < 0)
-            result = 1.0 / result;
-        return result;
-    }
-    return (float)exp((double)y * log((double)x));
 }
 
 double ldexp(double x, int exp) NOEXCEPT
 {
-    // FIXME: Please fix me. I am naive.
+    return x * exp2(exp);
-    double val = pow(2, exp);
-    return x * val;
 }
 
 float ldexpf(float x, int exp) NOEXCEPT
 {
-    // FIXME: Please fix me. I am naive.
+    return x * exp2f(exp);
-    float val = powf(2, exp);
-    return x * val;
 }
 
 double tanh(double x) NOEXCEPT
@@ -212,22 +198,28 @@ double sinh(double x) NOEXCEPT
 
 double log10(double x) NOEXCEPT
 {
-    return log(x) / M_LN10;
+    double ret = 0.0;
+    __asm__(
+        "fldlg2\n"
+        "fld %%st(1)\n"
+        "fyl2x\n"
+        "fstp %%st(1)"
+        : "=t"(ret)
+        : "0"(x));
+    return ret;
 }
 
 double log(double x) NOEXCEPT
 {
-    if (x < 0)
+    double ret = 0.0;
-        return NAN;
+    __asm__(
-    if (x == 0)
+        "fldln2\n"
-        return -INFINITY;
+        "fld %%st(1)\n"
-    double y = 1 + 2 * (x - 1) / (x + 1);
+        "fyl2x\n"
-    double exponentiated = exp(y);
+        "fstp %%st(1)"
-    y = y + 2 * (x - exponentiated) / (x + exponentiated);
+        : "=t"(ret)
-    exponentiated = exp(y);
+        : "0"(x));
-    y = y + 2 * (x - exponentiated) / (x + exponentiated);
+    return ret;
-    exponentiated = exp(y);
-    return y + 2 * (x - exponentiated) / (x + exponentiated);
 }
 
 float logf(float x) NOEXCEPT
@@ -247,38 +239,19 @@ float fmodf(float index, float period) NOEXCEPT
 
 double exp(double exponent) NOEXCEPT
 {
-    double result = 1;
+    double res = 0;
-    if (exponent >= 1) {
+    __asm__("fldl2e\n"
-        size_t integer_part = (size_t)exponent;
+            "fmulp\n"
-        if (integer_part & 1)
+            "fld1\n"
-            result *= e_to_power<1>();
+            "fld %%st(1)\n"
-        if (integer_part & 2)
+            "fprem\n"
-            result *= e_to_power<2>();
+            "f2xm1\n"
-        if (integer_part > 3) {
+            "faddp\n"
-            if (integer_part & 4)
+            "fscale\n"
-                result *= e_to_power<4>();
+            "fstp %%st(1)"
-            if (integer_part & 8)
+            : "=t"(res)
-                result *= e_to_power<8>();
+            : "0"(exponent));
-            if (integer_part & 16)
+    return res;
-                result *= e_to_power<16>();
-            if (integer_part & 32)
-                result *= e_to_power<32>();
-            if (integer_part >= 64)
-                return 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;
 }
 
 float expf(float exponent) NOEXCEPT
@@ -288,12 +261,22 @@ float expf(float exponent) NOEXCEPT
 
 double exp2(double exponent) NOEXCEPT
 {
-    return pow(2.0, exponent);
+    double res = 0;
+    __asm__("fld1\n"
+            "fld %%st(1)\n"
+            "fprem\n"
+            "f2xm1\n"
+            "faddp\n"
+            "fscale\n"
+            "fstp %%st(1)"
+            : "=t"(res)
+            : "0"(exponent));
+    return res;
 }
 
 float exp2f(float exponent) NOEXCEPT
 {
-    return pow(2.0f, exponent);
+    return (float)exp2(exponent);
 }
 
 double cosh(double x) NOEXCEPT
@@ -380,7 +363,15 @@ double fabs(double value) NOEXCEPT
 
 double log2(double x) NOEXCEPT
 {
-    return log(x) / M_LN2;
+    double ret = 0.0;
+    __asm__(
+        "fld1\n"
+        "fld %%st(1)\n"
+        "fyl2x\n"
+        "fstp %%st(1)"
+        : "=t"(ret)
+        : "0"(x));
+    return ret;
 }
 
 float log2f(float x) NOEXCEPT
@@ -490,16 +481,45 @@ double gamma(double x) NOEXCEPT
 
 double expm1(double x) NOEXCEPT
 {
-    return pow(M_E, x) - 1;
+    return exp(x) - 1;
 }
 
 double cbrt(double x) NOEXCEPT
 {
-    if (x > 0) {
+    if (isinf(x) || x == 0)
-        return pow(x, 1.0 / 3.0);
+        return x;
+    if (x < 0)
+        return -cbrt(-x);
+
+    double r = x;
+    double ex = 0;
+
+    while (r < 0.125) {
+        r *= 8;
+        ex--;
+    }
+    while (r > 1.0) {
+        r *= 0.125;
+        ex++;
     }
 
-    return -pow(-x, 1.0 / 3.0);
+    r = (-0.46946116 * r + 1.072302) * r + 0.3812513;
+
+    while (ex < 0) {
+        r *= 0.5;
+        ex++;
+    }
+    while (ex > 0) {
+        r *= 2;
+        ex--;
+    }
+
+    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
+    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
+    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
+    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
+
+    return r;
 }
 
 double log1p(double x) NOEXCEPT