AK: Introduce Math.h

This is to implement constexpr template based implementations for
mathematical functions

This also changes math.cpp to use these implementations.

Also adds a fastpath for floating point trucation for values smaller
than the signed 64 bit limit.

AK/Math.h

@@ -0,0 +1,467 @@
+/*
+ * Copyright (c) 2021, Leon Albrecht <leon2002.la@gmail.com>.
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#pragma once
+
+#include <AK/Concepts.h>
+#include <AK/StdLibExtraDetails.h>
+#include <AK/Types.h>
+
+namespace AK {
+
+template<FloatingPoint T>
+constexpr T NaN = __builtin_nan("");
+template<FloatingPoint T>
+constexpr T Pi = 3.141592653589793238462643383279502884L;
+template<FloatingPoint T>
+constexpr T E = 2.718281828459045235360287471352662498L;
+
+namespace Details {
+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>(); }
+}
+
+#define CONSTEXPR_STATE(function, args...)        \
+    if (is_constant_evaluated()) {                \
+        if (IsSame<T, long double>)               \
+            return __builtin_##function##l(args); \
+        if (IsSame<T, double>)                    \
+            return __builtin_##function(args);    \
+        if (IsSame<T, float>)                     \
+            return __builtin_##function##f(args); \
+    }
+
+#define INTEGER_BUILTIN(name)                         \
+    template<Integral T>                              \
+    constexpr T name(T x)                             \
+    {                                                 \
+        if constexpr (sizeof(T) == sizeof(long long)) \
+            return __builtin_##name##ll(x);           \
+        if constexpr (sizeof(T) == sizeof(long))      \
+            return __builtin_##name##l(x);            \
+        return __builtin_##name(x);                   \
+    }
+
+INTEGER_BUILTIN(clz);
+INTEGER_BUILTIN(ctz);
+INTEGER_BUILTIN(popcnt);
+
+namespace Division {
+template<FloatingPoint T>
+constexpr T fmod(T x, T y)
+{
+    CONSTEXPR_STATE(fmod, x, y);
+    T res;
+    asm(
+        "fprem"
+        : "=t"(res)
+        : "0"(x), "u"(y));
+    return res;
+}
+template<FloatingPoint T>
+constexpr T remainder(T x, T y)
+{
+    CONSTEXPR_STATE(remainder, x, y);
+    T res;
+    asm(
+        "fprem1"
+        : "=t"(res)
+        : "0"(x), "u"(y));
+    return res;
+}
+}
+
+using Division::fmod;
+using Division::remainder;
+
+template<FloatingPoint T>
+constexpr T sqrt(T x)
+{
+    CONSTEXPR_STATE(sqrt, x);
+    T res;
+    asm("fsqrt"
+        : "=t"(res)
+        : "0"(x));
+    return res;
+}
+
+template<FloatingPoint T>
+constexpr T cbrt(T x)
+{
+    CONSTEXPR_STATE(cbrt, x);
+    if (__builtin_isinf(x) || x == 0)
+        return x;
+    if (x < 0)
+        return -cbrt(-x);
+
+    T r = x;
+    T ex = 0;
+
+    while (r < 0.125l) {
+        r *= 8;
+        ex--;
+    }
+    while (r > 1.0l) {
+        r *= 0.125l;
+        ex++;
+    }
+
+    r = (-0.46946116l * r + 1.072302l) * r + 0.3812513l;
+
+    while (ex < 0) {
+        r *= 0.5l;
+        ex++;
+    }
+    while (ex > 0) {
+        r *= 2.0l;
+        ex--;
+    }
+
+    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
+    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
+    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
+    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
+
+    return r;
+}
+
+template<FloatingPoint T>
+constexpr T fabs(T x)
+{
+    if (is_constant_evaluated())
+        return x < 0 ? -x : x;
+    asm(
+        "fabs"
+        : "+t"(x));
+    return x;
+}
+
+namespace Trigonometry {
+
+template<FloatingPoint T>
+constexpr T hypot(T x, T y)
+{
+    return sqrt(x * x + y * y);
+}
+
+template<FloatingPoint T>
+constexpr T sin(T angle)
+{
+    CONSTEXPR_STATE(sin, angle);
+    T ret;
+    asm(
+        "fsin"
+        : "=t"(ret)
+        : "0"(angle));
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T cos(T angle)
+{
+    CONSTEXPR_STATE(cos, angle);
+    T ret;
+    asm(
+        "fcos"
+        : "=t"(ret)
+        : "0"(angle));
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T tan(T angle)
+{
+    CONSTEXPR_STATE(tan, angle);
+    double ret, one;
+    asm(
+        "fptan"
+        : "=t"(one), "=u"(ret)
+        : "0"(angle));
+
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T atan(T value)
+{
+    CONSTEXPR_STATE(atan, value);
+
+    T ret;
+    asm(
+        "fld1\n"
+        "fpatan\n"
+        : "=t"(ret)
+        : "0"(value));
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T asin(T x)
+{
+    CONSTEXPR_STATE(asin, x);
+    if (x > 1 || x < -1)
+        return NaN<T>;
+    if (x > (T)0.5 || x < (T)-0.5)
+        return 2 * atan<T>(x / (1 + sqrt<T>(1 - x * x)));
+    T squared = x * x;
+    T value = x;
+    T i = x * squared;
+    value += i * Details::product_odd<1>() / Details::product_even<2>() / 3;
+    i *= squared;
+    value += i * Details::product_odd<3>() / Details::product_even<4>() / 5;
+    i *= squared;
+    value += i * Details::product_odd<5>() / Details::product_even<6>() / 7;
+    i *= squared;
+    value += i * Details::product_odd<7>() / Details::product_even<8>() / 9;
+    i *= squared;
+    value += i * Details::product_odd<9>() / Details::product_even<10>() / 11;
+    i *= squared;
+    value += i * Details::product_odd<11>() / Details::product_even<12>() / 13;
+    i *= squared;
+    value += i * Details::product_odd<13>() / Details::product_even<14>() / 15;
+    i *= squared;
+    value += i * Details::product_odd<15>() / Details::product_even<16>() / 17;
+    return value;
+}
+
+template<FloatingPoint T>
+constexpr T acos(T value)
+{
+    CONSTEXPR_STATE(acos, value);
+
+    // FIXME: I am naive
+    return Pi<T> + asin(value);
+}
+
+template<FloatingPoint T>
+constexpr T atan2(T y, T x)
+{
+    CONSTEXPR_STATE(atan2, y, x);
+
+    T ret;
+    asm("fpatan"
+        : "=t"(ret)
+        : "0"(x), "u"(y)
+        : "st(1)");
+    return ret;
+}
+
+}
+
+using Trigonometry::acos;
+using Trigonometry::asin;
+using Trigonometry::atan;
+using Trigonometry::atan2;
+using Trigonometry::cos;
+using Trigonometry::hypot;
+using Trigonometry::sin;
+using Trigonometry::tan;
+
+namespace Exponentials {
+
+template<FloatingPoint T>
+constexpr T log(T x)
+{
+    CONSTEXPR_STATE(log, x);
+
+    T ret;
+    asm(
+        "fldln2\n"
+        "fxch %%st(1)\n"
+        "fyl2x\n"
+        : "=t"(ret)
+        : "0"(x));
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T log2(T x)
+{
+    CONSTEXPR_STATE(log2, x);
+
+    T ret;
+    asm(
+        "fld1\n"
+        "fxch %%st(1)\n"
+        "fyl2x\n"
+        : "=t"(ret)
+        : "0"(x));
+    return ret;
+}
+
+template<Integral T>
+constexpr T log2(T x)
+{
+    return x ? 8 * sizeof(T) - clz(x) : 0;
+}
+
+template<FloatingPoint T>
+constexpr T log10(T x)
+{
+    CONSTEXPR_STATE(log10, x);
+
+    T ret;
+    asm(
+        "fldlg2\n"
+        "fxch %%st(1)\n"
+        "fyl2x\n"
+        : "=t"(ret)
+        : "0"(x));
+    return ret;
+}
+
+template<FloatingPoint T>
+constexpr T exp(T exponent)
+{
+    CONSTEXPR_STATE(exp, exponent);
+
+    T res;
+    asm("fldl2e\n"
+        "fmulp\n"
+        "fld1\n"
+        "fld %%st(1)\n"
+        "fprem\n"
+        "f2xm1\n"
+        "faddp\n"
+        "fscale\n"
+        "fstp %%st(1)"
+        : "=t"(res)
+        : "0"(exponent));
+    return res;
+}
+
+template<FloatingPoint T>
+constexpr T exp2(T exponent)
+{
+    CONSTEXPR_STATE(exp2, exponent);
+
+    T res;
+    asm("fld1\n"
+        "fld %%st(1)\n"
+        "fprem\n"
+        "f2xm1\n"
+        "faddp\n"
+        "fscale\n"
+        "fstp %%st(1)"
+        : "=t"(res)
+        : "0"(exponent));
+    return res;
+}
+template<Integral T>
+constexpr T exp2(T exponent)
+{
+    return 1u << exponent;
+}
+}
+
+using Exponentials::exp;
+using Exponentials::exp2;
+using Exponentials::log;
+using Exponentials::log10;
+using Exponentials::log2;
+
+namespace Hyperbolic {
+
+template<FloatingPoint T>
+constexpr T sinh(T x)
+{
+    T exponentiated = exp<T>(x);
+    if (x > 0)
+        return (exponentiated * exponentiated - 1) / 2 / exponentiated;
+    return (exponentiated - 1 / exponentiated) / 2;
+}
+
+template<FloatingPoint T>
+constexpr T cosh(T x)
+{
+    CONSTEXPR_STATE(cosh, x);
+
+    T exponentiated = exp(-x);
+    if (x < 0)
+        return (1 + exponentiated * exponentiated) / 2 / exponentiated;
+    return (1 / exponentiated + exponentiated) / 2;
+}
+
+template<FloatingPoint T>
+constexpr T tanh(T x)
+{
+    if (x > 0) {
+        T exponentiated = exp<T>(2 * x);
+        return (exponentiated - 1) / (exponentiated + 1);
+    }
+    T plusX = exp<T>(x);
+    T minusX = 1 / plusX;
+    return (plusX - minusX) / (plusX + minusX);
+}
+
+template<FloatingPoint T>
+constexpr T asinh(T x)
+{
+    return log<T>(x + sqrt<T>(x * x + 1));
+}
+
+template<FloatingPoint T>
+constexpr T acosh(T x)
+{
+    return log<T>(x + sqrt<T>(x * x - 1));
+}
+
+template<FloatingPoint T>
+constexpr T atanh(T x)
+{
+    return log<T>((1 + x) / (1 - x)) / (T)2.0l;
+}
+
+}
+
+using Hyperbolic::acosh;
+using Hyperbolic::asinh;
+using Hyperbolic::atanh;
+using Hyperbolic::cosh;
+using Hyperbolic::sinh;
+using Hyperbolic::tanh;
+
+template<FloatingPoint T>
+constexpr T pow(T x, T y)
+{
+    CONSTEXPR_STATE(pow, x, y);
+    // fixme I am naive
+    if (__builtin_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;
+    if (y == (T)y_as_int) {
+        T result = x;
+        for (int i = 0; i < fabs<T>(y) - 1; ++i)
+            result *= x;
+        if (y < 0)
+            result = 1.0l / result;
+        return result;
+    }
+
+    return exp2<T>(y * log2<T>(x));
+}
+
+#undef CONSTEXPR_STATE
+#undef INTEGER_BUILTIN
+
+}

Userland/Libraries/LibM/math.cpp

@@ -6,6 +6,7 @@
  */
 
 #include <AK/ExtraMathConstants.h>
+#include <AK/Math.h>
 #include <AK/Platform.h>
 #include <AK/StdLibExtras.h>
 #include <LibC/assert.h>
@@ -345,124 +346,196 @@ long double nanl(const char* s) NOEXCEPT
     return __builtin_nanl(s);
 }
 
-double trunc(double x) NOEXCEPT
+#define MAKE_AK_BACKED1(name)                     \
+    long double name##l(long double arg) NOEXCEPT \
+    {                                             \
+        return AK::name<long double>(arg);        \
+    }                                             \
+    double name(double arg) NOEXCEPT              \
+    {                                             \
+        return AK::name<double>(arg);             \
+    }                                             \
+    float name##f(float arg) NOEXCEPT             \
+    {                                             \
+        return AK::name<float>(arg);              \
+    }
+#define MAKE_AK_BACKED2(name)                                        \
+    long double name##l(long double arg1, long double arg2) NOEXCEPT \
+    {                                                                \
+        return AK::name<long double>(arg1, arg2);                    \
+    }                                                                \
+    double name(double arg1, double arg2) NOEXCEPT                   \
+    {                                                                \
+        return AK::name<double>(arg1, arg2);                         \
+    }                                                                \
+    float name##f(float arg1, float arg2) NOEXCEPT                   \
+    {                                                                \
+        return AK::name<float>(arg1, arg2);                          \
+    }
+
+MAKE_AK_BACKED1(sin);
+MAKE_AK_BACKED1(cos);
+MAKE_AK_BACKED1(tan);
+MAKE_AK_BACKED1(asin);
+MAKE_AK_BACKED1(acos);
+MAKE_AK_BACKED1(atan);
+MAKE_AK_BACKED1(sinh);
+MAKE_AK_BACKED1(cosh);
+MAKE_AK_BACKED1(tanh);
+MAKE_AK_BACKED1(asinh);
+MAKE_AK_BACKED1(acosh);
+MAKE_AK_BACKED1(atanh);
+MAKE_AK_BACKED1(sqrt);
+MAKE_AK_BACKED1(cbrt);
+MAKE_AK_BACKED1(log);
+MAKE_AK_BACKED1(log2);
+MAKE_AK_BACKED1(log10);
+MAKE_AK_BACKED1(exp);
+MAKE_AK_BACKED1(exp2);
+MAKE_AK_BACKED1(fabs);
+
+MAKE_AK_BACKED2(atan2);
+MAKE_AK_BACKED2(hypot);
+MAKE_AK_BACKED2(fmod);
+MAKE_AK_BACKED2(pow);
+MAKE_AK_BACKED2(remainder);
+
+long double truncl(long double x) NOEXCEPT
 {
+    if (fabsl(x) < LONG_LONG_MAX) {
+        // This is 1.6 times faster than the implemenation using the "internal_to_integer"
+        // helper (on x86_64)
+        // https://quick-bench.com/q/xBmxuY8am9qibSYVna90Y6PIvqA
+        u64 temp;
+        asm(
+            "fisttpq %[temp]\n"
+            "fildq %[temp]"
+            : "+t"(x)
+            : [temp] "m"(temp));
+        return x;
+    }
+
     return internal_to_integer(x, RoundingMode::ToZero);
 }
 
-float truncf(float x) NOEXCEPT
+double trunc(double x) NOEXCEPT
 {
+    if (fabs(x) < LONG_LONG_MAX) {
+        u64 temp;
+        asm(
+            "fisttpq %[temp]\n"
+            "fildq %[temp]"
+            : "+t"(x)
+            : [temp] "m"(temp));
+        return x;
+    }
+
     return internal_to_integer(x, RoundingMode::ToZero);
 }
 
-long double truncl(long double x) NOEXCEPT
+float truncf(float x) NOEXCEPT
 {
+    if (fabsf(x) < LONG_LONG_MAX) {
+        u64 temp;
+        asm(
+            "fisttpq %[temp]\n"
+            "fildq %[temp]"
+            : "+t"(x)
+            : [temp] "m"(temp));
+        return x;
+    }
+
     return internal_to_integer(x, RoundingMode::ToZero);
 }
 
-long double cosl(long double angle) NOEXCEPT
+long double rintl(long double value)
 {
-    long double ret = 0.0;
+    double res;
     asm(
-        "fcos"
-        : "=t"(ret)
-        : "0"(angle));
-    return ret;
+        "frndint\n"
+        : "=t"(res)
+        : "0"(value));
+    return res;
 }
-
-double cos(double angle) NOEXCEPT
+double rint(double value)
 {
-    double ret = 0.0;
+    double res;
     asm(
-        "fcos"
-        : "=t"(ret)
-        : "0"(angle));
-
-    return ret;
+        "frndint\n"
+        : "=t"(res)
+        : "0"(value));
+    return res;
 }
-
-float cosf(float angle) NOEXCEPT
+float rintf(float value)
 {
-    float ret = 0.0;
+    double res;
     asm(
-        "fcos"
-        : "=t"(ret)
-        : "0"(angle));
-
-    return ret;
+        "frndint\n"
+        : "=t"(res)
+        : "0"(value));
+    return res;
 }
 
-long double sinl(long double angle) NOEXCEPT
+long lrintl(long double value)
 {
-    long double ret = 0.0;
+    long res;
     asm(
-        "fsin"
-        : "=t"(ret)
-        : "0"(angle));
-
-    return ret;
+        "fistpl %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
-
-// This can also be done with a taylor expansion, but for
-// now this works pretty well (and doesn't mess anything up
-// in quake in particular, which is very Floating-Point precision
-// heavy)
-double sin(double angle) NOEXCEPT
+long lrint(double value)
 {
-    double ret = 0.0;
+    long res;
     asm(
-        "fsin"
-        : "=t"(ret)
-        : "0"(angle));
-
-    return ret;
+        "fistpl %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
-
-float sinf(float angle) NOEXCEPT
+long lrintf(float value)
 {
-    float ret = 0.0f;
+    long res;
     asm(
-        "fsin"
-        : "=t"(ret)
-        : "0"(angle));
-    return ret;
+        "fistpl %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
 
-long double powl(long double x, long double y) NOEXCEPT
+long long llrintl(long double value)
 {
-    // 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;
-    if (y == (long double)y_as_int) {
-        long double result = x;
-        for (int i = 0; i < fabsl(y) - 1; ++i)
-            result *= x;
-        if (y < 0)
-            result = 1.0l / result;
-        return result;
-    }
-    if (x < 0) {
-        return 1.l / exp2l(y * log2l(-x));
-    }
-
-    return exp2l(y * log2l(x));
+    long long res;
+    asm(
+        "fistpq %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
-
-double pow(double x, double y) NOEXCEPT
+long long llrint(double value)
 {
-    return (double)powl(x, y);
+    long long res;
+    asm(
+        "fistpq %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
-
-float powf(float x, float y) NOEXCEPT
+long long llrintf(float value)
 {
-    return (float)powl(x, y);
+    long long res;
+    asm(
+        "fistpq %0\n"
+        : "+m"(res)
+        : "t"(value)
+        : "st");
+    return res;
 }
 
 // On systems where FLT_RADIX == 2, ldexp is equivalent to scalbn
@@ -481,27 +554,6 @@ float ldexpf(float x, int exp) NOEXCEPT
     return internal_scalbn(x, exp);
 }
 
-long double tanhl(long double x) NOEXCEPT
-{
-    if (x > 0) {
-        long double exponentiated = expl(2 * x);
-        return (exponentiated - 1) / (exponentiated + 1);
-    }
-    long double plusX = expl(x);
-    long double minusX = 1 / plusX;
-    return (plusX - minusX) / (plusX + minusX);
-}
-
-double tanh(double x) NOEXCEPT
-{
-    return (double)tanhl(x);
-}
-
-float tanhf(float x) NOEXCEPT
-{
-    return (float)tanhl(x);
-}
-
 [[maybe_unused]] static long double ampsin(long double angle) NOEXCEPT
 {
     long double looped_angle = fmodl(M_PI + angle, M_TAU) - M_PI;
@@ -519,345 +571,6 @@ float tanhf(float x) NOEXCEPT
     return quadratic_term + linear_term;
 }
 
-long double tanl(long double angle) NOEXCEPT
-{
-    long double ret = 0.0, one;
-    __asm__(
-        "fptan"
-        : "=t"(one), "=u"(ret)
-        : "0"(angle));
-
-    return ret;
-}
-
-double tan(double angle) NOEXCEPT
-{
-    return (double)tanl(angle);
-}
-
-float tanf(float angle) NOEXCEPT
-{
-    return (float)tanl(angle);
-}
-
-long double sqrtl(long double x) NOEXCEPT
-{
-    long double res;
-    asm("fsqrt"
-        : "=t"(res)
-        : "0"(x));
-    return res;
-}
-
-double sqrt(double x) NOEXCEPT
-{
-    double res;
-    __asm__("fsqrt"
-            : "=t"(res)
-            : "0"(x));
-    return res;
-}
-
-float sqrtf(float x) NOEXCEPT
-{
-    float res;
-    __asm__("fsqrt"
-            : "=t"(res)
-            : "0"(x));
-    return res;
-}
-
-long double sinhl(long double x) NOEXCEPT
-{
-    long double exponentiated = expl(x);
-    if (x > 0)
-        return (exponentiated * exponentiated - 1) / 2 / exponentiated;
-    return (exponentiated - 1 / exponentiated) / 2;
-}
-
-double sinh(double x) NOEXCEPT
-{
-    return (double)sinhl(x);
-}
-
-float sinhf(float x) NOEXCEPT
-{
-    return (float)sinhl(x);
-}
-
-long double log10l(long double x) NOEXCEPT
-{
-    long double ret = 0.0l;
-    __asm__(
-        "fldlg2\n"
-        "fld %%st(1)\n"
-        "fyl2x\n"
-        "fstp %%st(1)"
-        : "=t"(ret)
-        : "0"(x));
-    return ret;
-}
-
-double log10(double x) NOEXCEPT
-{
-    return (double)log10l(x);
-}
-
-float log10f(float x) NOEXCEPT
-{
-    return (float)log10l(x);
-}
-
-long double logl(long double x) NOEXCEPT
-{
-    long double ret = 0.0l;
-    asm(
-        "fldln2\n"
-        "fld %%st(1)\n"
-        "fyl2x\n"
-        "fstp %%st(1)"
-        : "=t"(ret)
-        : "0"(x));
-    return ret;
-}
-
-double log(double x) NOEXCEPT
-{
-    return (double)logl(x);
-}
-
-float logf(float x) NOEXCEPT
-{
-    return (float)logl(x);
-}
-
-long double fmodl(long double index, long double period) NOEXCEPT
-{
-    return index - truncl(index / period) * period;
-}
-
-double fmod(double index, double period) NOEXCEPT
-{
-    return index - trunc(index / period) * period;
-}
-
-float fmodf(float index, float period) NOEXCEPT
-{
-    return index - truncf(index / period) * period;
-}
-
-// FIXME: These aren't exactly like fmod, but these definitions are probably good enough for now
-long double remainderl(long double x, long double y) NOEXCEPT
-{
-    return fmodl(x, y);
-}
-
-double remainder(double x, double y) NOEXCEPT
-{
-    return fmod(x, y);
-}
-
-float remainderf(float x, float y) NOEXCEPT
-{
-    return fmodf(x, y);
-}
-
-long double expl(long double exponent) NOEXCEPT
-{
-    long double res = 0;
-    asm("fldl2e\n"
-        "fmulp\n"
-        "fld1\n"
-        "fld %%st(1)\n"
-        "fprem\n"
-        "f2xm1\n"
-        "faddp\n"
-        "fscale\n"
-        "fstp %%st(1)"
-        : "=t"(res)
-        : "0"(exponent));
-    return res;
-}
-
-double exp(double exponent) NOEXCEPT
-{
-    return (double)expl(exponent);
-}
-
-float expf(float exponent) NOEXCEPT
-{
-    return (float)expl(exponent);
-}
-
-long double exp2l(long double exponent) NOEXCEPT
-{
-    long 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;
-}
-
-double exp2(double exponent) NOEXCEPT
-{
-    return (double)exp2l(exponent);
-}
-
-float exp2f(float exponent) NOEXCEPT
-{
-    return (float)exp2l(exponent);
-}
-
-long double coshl(long double x) NOEXCEPT
-{
-    long double exponentiated = expl(-x);
-    if (x < 0)
-        return (1 + exponentiated * exponentiated) / 2 / exponentiated;
-    return (1 / exponentiated + exponentiated) / 2;
-}
-
-double cosh(double x) NOEXCEPT
-{
-    return (double)coshl(x);
-}
-
-float coshf(float x) NOEXCEPT
-{
-    return (float)coshl(x);
-}
-
-long double atan2l(long double y, long double x) NOEXCEPT
-{
-    long double result = 0; //atanl(y / x);
-    asm("fpatan"
-        : "=t"(result)
-        : "0"(x), "u"(y)
-        : "st(1)");
-    return result;
-}
-
-double atan2(double y, double x) NOEXCEPT
-{
-    double result = 0; //atanl(y / x);
-    asm("fpatan"
-        : "=t"(result)
-        : "0"(x), "u"(y)
-        : "st(1)");
-    return result;
-}
-float atan2f(float y, float x) NOEXCEPT
-{
-    float result = 0; //atanl(y / x);
-    asm("fpatan"
-        : "=t"(result)
-        : "0"(x), "u"(y)
-        : "st(1)");
-    return result;
-}
-
-long double atanl(long double x) NOEXCEPT
-{
-    asm(
-        "fld1\n"
-        "fpatan\n"
-        : "=t"(x)
-        : "0"(x));
-    return x;
-}
-
-double atan(double x) NOEXCEPT
-{
-    asm(
-        "fld1\n"
-        "fpatan\n"
-        : "=t"(x)
-        : "0"(x));
-    return x;
-}
-
-float atanf(float x) NOEXCEPT
-{
-    asm(
-        "fld1\n"
-        "fpatan\n"
-        : "=t"(x)
-        : "0"(x));
-    return x;
-}
-
-long double asinl(long double x) NOEXCEPT
-{
-    if (x > 1 || x < -1)
-        return NAN;
-    if (x > 0.5 || x < -0.5)
-        return 2 * atanl(x / (1 + sqrtl(1 - x * x)));
-    long double squared = x * x;
-    long double value = x;
-    long 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;
-    i *= squared;
-    value += i * product_odd<13>() / product_even<14>() / 15;
-    i *= squared;
-    value += i * product_odd<15>() / product_even<16>() / 17;
-    return value;
-}
-
-double asin(double x) NOEXCEPT
-{
-    return (double)asinl(x);
-}
-
-float asinf(float x) NOEXCEPT
-{
-    return (float)asinl(x);
-}
-
-long double acosl(long double x) NOEXCEPT
-{
-    return M_PI_2 - asinl(x);
-}
-
-double acos(double x) NOEXCEPT
-{
-    return M_PI_2 - asin(x);
-}
-
-float acosf(float x) NOEXCEPT
-{
-    return static_cast<float>(M_PI_2) - asinf(x);
-}
-
-long double fabsl(long double value) NOEXCEPT
-{
-    return value < 0 ? -value : value;
-}
-
-double fabs(double value) NOEXCEPT
-{
-    return value < 0 ? -value : value;
-}
-
-float fabsf(float value) NOEXCEPT
-{
-    return value < 0 ? -value : value;
-}
-
 int ilogbl(long double x) NOEXCEPT
 {
     return internal_ilogb(x);
@@ -888,29 +601,6 @@ float logbf(float x) NOEXCEPT
     return ilogbf(x);
 }
 
-long double log2l(long double x) NOEXCEPT
-{
-    long double ret = 0.0;
-    asm(
-        "fld1\n"
-        "fld %%st(1)\n"
-        "fyl2x\n"
-        "fstp %%st(1)"
-        : "=t"(ret)
-        : "0"(x));
-    return ret;
-}
-
-double log2(double x) NOEXCEPT
-{
-    return (double)log2l(x);
-}
-
-float log2f(float x) NOEXCEPT
-{
-    return (float)log2l(x);
-}
-
 double frexp(double x, int* exp) NOEXCEPT
 {
     *exp = (x == 0) ? 0 : (1 + ilogb(x));
@@ -989,51 +679,6 @@ long double floorl(long double value) NOEXCEPT
     return internal_to_integer(value, RoundingMode::Down);
 }
 
-long double rintl(long double value) NOEXCEPT
-{
-    return internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-double rint(double value) NOEXCEPT
-{
-    return internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-float rintf(float value) NOEXCEPT
-{
-    return internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long lrintl(long double value) NOEXCEPT
-{
-    return (long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long lrint(double value) NOEXCEPT
-{
-    return (long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long lrintf(float value) NOEXCEPT
-{
-    return (long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long long llrintl(long double value) NOEXCEPT
-{
-    return (long long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long long llrint(double value) NOEXCEPT
-{
-    return (long long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
-long long llrintf(float value) NOEXCEPT
-{
-    return (long long)internal_to_integer(value, RoundingMode { fegetround() });
-}
-
 float ceilf(float value) NOEXCEPT
 {
     return internal_to_integer(value, RoundingMode::Up);
@@ -1150,54 +795,6 @@ float expm1f(float x) NOEXCEPT
     return expf(x) - 1;
 }
 
-long double cbrtl(long double x) NOEXCEPT
-{
-    if (isinf(x) || x == 0)
-        return x;
-    if (x < 0)
-        return -cbrtl(-x);
-
-    long double r = x;
-    long double ex = 0;
-
-    while (r < 0.125l) {
-        r *= 8;
-        ex--;
-    }
-    while (r > 1.0l) {
-        r *= 0.125l;
-        ex++;
-    }
-
-    r = (-0.46946116l * r + 1.072302l) * r + 0.3812513l;
-
-    while (ex < 0) {
-        r *= 0.5l;
-        ex++;
-    }
-    while (ex > 0) {
-        r *= 2.0l;
-        ex--;
-    }
-
-    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
-    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
-    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
-    r = (2.0l / 3.0l) * r + (1.0l / 3.0l) * x / (r * r);
-
-    return r;
-}
-
-double cbrt(double x) NOEXCEPT
-{
-    return (double)cbrtl(x);
-}
-
-float cbrtf(float x) NOEXCEPT
-{
-    return (float)cbrtl(x);
-}
-
 long double log1pl(long double x) NOEXCEPT
 {
     return logl(1 + x);
@@ -1213,66 +810,6 @@ float log1pf(float x) NOEXCEPT
     return logf(1 + x);
 }
 
-long double acoshl(long double x) NOEXCEPT
-{
-    return logl(x + sqrtl(x * x - 1));
-}
-
-double acosh(double x) NOEXCEPT
-{
-    return log(x + sqrt(x * x - 1));
-}
-
-float acoshf(float x) NOEXCEPT
-{
-    return logf(x + sqrtf(x * x - 1));
-}
-
-long double asinhl(long double x) NOEXCEPT
-{
-    return logl(x + sqrtl(x * x + 1));
-}
-
-double asinh(double x) NOEXCEPT
-{
-    return log(x + sqrt(x * x + 1));
-}
-
-float asinhf(float x) NOEXCEPT
-{
-    return logf(x + sqrtf(x * x + 1));
-}
-
-long double atanhl(long double x) NOEXCEPT
-{
-    return logl((1 + x) / (1 - x)) / 2.0l;
-}
-
-double atanh(double x) NOEXCEPT
-{
-    return log((1 + x) / (1 - x)) / 2.0;
-}
-
-float atanhf(float x) NOEXCEPT
-{
-    return logf((1 + x) / (1 - x)) / 2.0f;
-}
-
-long double hypotl(long double x, long double y) NOEXCEPT
-{
-    return sqrtl(x * x + y * y);
-}
-
-double hypot(double x, double y) NOEXCEPT
-{
-    return sqrt(x * x + y * y);
-}
-
-float hypotf(float x, float y) NOEXCEPT
-{
-    return sqrtf(x * x + y * y);
-}
-
 long double erfl(long double x) NOEXCEPT
 {
     // algorithm taken from Abramowitz and Stegun (no. 26.2.17)