diff --git a/Userspace/libc/include/math.h b/Userspace/libc/include/math.h
index 0fa7dcde..d2a7f0c2 100644
--- a/Userspace/libc/include/math.h
+++ b/Userspace/libc/include/math.h
@@ -255,7 +255,7 @@ typedef long double double_t;
 	double round(double x);
 	float roundf(float x);
 	long double roundl(long double x);
-	double scalb(double x, double y);
+	double scalb(double x, double n);
 	double scalbln(double x, long n);
 	float scalblnf(float x, long n);
 	long double scalblnl(long double x, long n);
diff --git a/Userspace/libc/src/std/math.c b/Userspace/libc/src/std/math.c
index c11ec62f..a18a99f5 100644
--- a/Userspace/libc/src/std/math.c
+++ b/Userspace/libc/src/std/math.c
@@ -19,184 +19,1683 @@
 #include <sys/types.h>
 #include <float.h>
 #include <errno.h>
+#include <inttypes.h>
+#include <fenv.h>
 
 export int signgam;
 
-export double acos(double x);
-export float acosf(float x);
-export double acosh(double x);
-export float acoshf(float x);
-export long double acoshl(long double x);
-export long double acosl(long double x);
-export double asin(double x);
-export float asinf(float x);
-export double asinh(double x);
-export float asinhf(float x);
-export long double asinhl(long double x);
-export long double asinl(long double x);
-export double atan(double x);
-export double atan2(double y, double x);
-export float atan2f(float y, float x);
-export long double atan2l(long double y, long double x);
-export float atanf(float x);
-export double atanh(double x);
-export float atanhf(float x);
-export long double atanhl(long double x);
-export long double atanl(long double x);
-export double cbrt(double x);
-export float cbrtf(float x);
-export long double cbrtl(long double x);
-export double ceil(double x);
-export float ceilf(float x);
-export long double ceill(long double x);
-export double copysign(double x, double y);
-export float copysignf(float x, float y);
-export long double copysignl(long double x, long double y);
-export double cos(double x);
-export float cosf(float x);
-export double cosh(double x);
-export float coshf(float x);
-export long double coshl(long double x);
-export long double cosl(long double x);
-export double erf(double x);
-export double erfc(double x);
-export float erfcf(float x);
-export long double erfcl(long double x);
-export float erff(float x);
-export long double erfl(long double x);
-export double exp(double x);
-export double exp2(double x);
-export float exp2f(float x);
-export long double exp2l(long double x);
-export float expf(float x);
-export long double expl(long double x);
-export double expm1(double x);
-export float expm1f(float x);
-export long double expm1l(long double x);
-export double fabs(double x);
-export float fabsf(float x);
-export long double fabsl(long double x);
-export double fdim(double x, double y);
-export float fdimf(float x, float y);
-export long double fdiml(long double x, long double y);
-export double floor(double x);
-export float floorf(float x);
-export long double floorl(long double x);
-export double fma(double x, double y, double z);
-export float fmaf(float x, float y, float z);
-export long double fmal(long double x, long double y, long double z);
-export double fmax(double x, double y);
-export float fmaxf(float x, float y);
-export long double fmaxl(long double x, long double y);
-export double fmin(double x, double y);
-export float fminf(float x, float y);
-export long double fminl(long double x, long double y);
-export double fmod(double x, double y);
-export float fmodf(float x, float y);
-export long double fmodl(long double x, long double y);
-export double frexp(double num, int *exp);
-export float frexpf(float num, int *exp);
-export long double frexpl(long double num, int *exp);
-export double hypot(double x, double y);
-export float hypotf(float x, float y);
-export long double hypotl(long double x, long double y);
-export int ilogb(double x);
-export int ilogbf(float x);
-export int ilogbl(long double x);
-export double j0(double x);
-export double j1(double x);
-export double jn(int n, double x);
-export double ldexp(double x, int exp);
-export float ldexpf(float x, int exp);
-export long double ldexpl(long double x, int exp);
-export double lgamma(double x);
-export float lgammaf(float x);
-export long double lgammal(long double x);
-export long long llrint(double x);
-export long long llrintf(float x);
-export long long llrintl(long double x);
-export long long llround(double x);
-export long long llroundf(float x);
-export long long llroundl(long double x);
-export double log(double x);
-export double log10(double x);
-export float log10f(float x);
-export long double log10l(long double x);
-export double log1p(double x);
-export float log1pf(float x);
-export long double log1pl(long double x);
-export double log2(double x);
-export float log2f(float x);
-export long double log2l(long double x);
-export double logb(double x);
-export float logbf(float x);
-export long double logbl(long double x);
-export float logf(float x);
-export long double logl(long double x);
-export long lrint(double x);
-export long lrintf(float x);
-export long lrintl(long double x);
-export long lround(double x);
-export long lroundf(float x);
-export long lroundl(long double x);
-export double modf(double x, double *iptr);
-export float modff(float value, float *iptr);
-export long double modfl(long double x, long double *iptr);
-export double nan(const char *tagp);
-export float nanf(const char *tagp);
-export long double nanl(const char *tagp);
-export double nearbyint(double x);
-export float nearbyintf(float x);
-export long double nearbyintl(long double x);
-export double nextafter(double x, double y);
-export float nextafterf(float x, float y);
-export long double nextafterl(long double x, long double y);
-export double nexttoward(double x, long double y);
-export float nexttowardf(float x, long double y);
-export long double nexttowardl(long double x, long double y);
-export double pow(double x, double y);
-export float powf(float x, float y);
-export long double powl(long double x, long double y);
-export double remainder(double x, double y);
-export float remainderf(float x, float y);
-export long double remainderl(long double x, long double y);
-export double remquo(double x, double y, int *quo);
-export float remquof(float x, float y, int *quo);
-export long double remquol(long double x, long double y, int *quo);
-export double rint(double x);
-export float rintf(float x);
-export long double rintl(long double x);
-export double round(double x);
-export float roundf(float x);
-export long double roundl(long double x);
-export double scalb(double x, double y);
-export double scalbln(double x, long n);
-export float scalblnf(float x, long n);
-export long double scalblnl(long double x, long n);
-export double scalbn(double x, int n);
-export float scalbnf(float x, int n);
-export long double scalbnl(long double x, int n);
-export double sin(double x);
-export float sinf(float x);
-export double sinh(double x);
-export float sinhf(float x);
-export long double sinhl(long double x);
-export long double sinl(long double x);
-export double sqrt(double x);
-export float sqrtf(float x);
-export long double sqrtl(long double x);
-export double tan(double x);
-export float tanf(float x);
-export double tanh(double x);
-export float tanhf(float x);
-export long double tanhl(long double x);
-export long double tanl(long double x);
-export double tgamma(double x);
-export float tgammaf(float x);
-export long double tgammal(long double x);
-export double trunc(double x);
-export float truncf(float x);
-export long double truncl(long double x);
-export double y0(double x);
-export double y1(double x);
-export double yn(int n, double x);
+export double acos(double x)
+{
+	if (x < -1.0 || x > 1.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2(sqrt(1.0 - x * x), x);
+}
+
+export float acosf(float x)
+{
+	if (x < -1.0f || x > 1.0f)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2f(sqrtf(1.0f - x * x), x);
+}
+
+export double acosh(double x)
+{
+	if (x < 1.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return log(x + sqrt(x * x - 1.0));
+}
+
+export float acoshf(float x)
+{
+	if (x < 1.0f)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return logf(x + sqrtf(x * x - 1.0f));
+}
+
+export long double acoshl(long double x)
+{
+	if (x < 1.0L)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return logl(x + sqrtl(x * x - 1.0L));
+}
+
+export long double acosl(long double x)
+{
+	if (x < -1.0L || x > 1.0L)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2l(sqrtl(1.0L - x * x), x);
+}
+
+export double asin(double x)
+{
+	if (x < -1.0 || x > 1.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2(x, sqrt(1.0 - x * x));
+}
+
+export float asinf(float x)
+{
+	if (x < -1.0f || x > 1.0f)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2f(x, sqrtf(1.0f - x * x));
+}
+
+export double asinh(double x)
+{
+	return log(x + sqrt(x * x + 1.0));
+}
+
+export float asinhf(float x)
+{
+	return logf(x + sqrtf(x * x + 1.0f));
+}
+
+export long double asinhl(long double x)
+{
+	return logl(x + sqrtl(x * x + 1.0L));
+}
+
+export long double asinl(long double x)
+{
+	if (x < -1.0L || x > 1.0L)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return atan2l(x, sqrtl(1.0L - x * x));
+}
+
+export double atan(double x)
+{
+	return atan2(x, 1.0);
+}
+
+export double atan2(double y, double x)
+{
+	if (x > 0.0)
+	{
+		return atan(y / x);
+	}
+	else if (x < 0.0 && y >= 0.0)
+	{
+		return atan(y / x) + M_PI;
+	}
+	else if (x < 0.0 && y < 0.0)
+	{
+		return atan(y / x) - M_PI;
+	}
+	else if (x == 0.0 && y > 0.0)
+	{
+		return M_PI_2;
+	}
+	else if (x == 0.0 && y < 0.0)
+	{
+		return -M_PI_2;
+	}
+	else
+	{
+		errno = EDOM;
+		return NAN;
+	}
+}
+
+export float atan2f(float y, float x)
+{
+	if (x > 0.0f)
+	{
+		return atanf(y / x);
+	}
+	else if (x < 0.0f && y >= 0.0f)
+	{
+		return atanf(y / x) + (float)M_PI;
+	}
+	else if (x < 0.0f && y < 0.0f)
+	{
+		return atanf(y / x) - (float)M_PI;
+	}
+	else if (x == 0.0f && y > 0.0f)
+	{
+		return (float)M_PI_2;
+	}
+	else if (x == 0.0f && y < 0.0f)
+	{
+		return -(float)M_PI_2;
+	}
+	else
+	{
+		errno = EDOM;
+		return NAN;
+	}
+}
+
+export long double atan2l(long double y, long double x)
+{
+	if (x > 0.0L)
+	{
+		return atanl(y / x);
+	}
+	else if (x < 0.0L && y >= 0.0L)
+	{
+		return atanl(y / x) + M_PIl;
+	}
+	else if (x < 0.0L && y < 0.0L)
+	{
+		return atanl(y / x) - M_PIl;
+	}
+	else if (x == 0.0L && y > 0.0L)
+	{
+		return M_PI_2l;
+	}
+	else if (x == 0.0L && y < 0.0L)
+	{
+		return -M_PI_2l;
+	}
+	else
+	{
+		errno = EDOM;
+		return NAN;
+	}
+}
+
+export float atanf(float x)
+{
+	return atan2f(x, 1.0f);
+}
+
+export double atanh(double x)
+{
+	if (x <= -1.0 || x >= 1.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return 0.5 * log((1.0 + x) / (1.0 - x));
+}
+
+export float atanhf(float x)
+{
+	if (x <= -1.0f || x >= 1.0f)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return 0.5f * logf((1.0f + x) / (1.0f - x));
+}
+
+export long double atanhl(long double x)
+{
+	if (x <= -1.0L || x >= 1.0L)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return 0.5L * logl((1.0L + x) / (1.0L - x));
+}
+
+export long double atanl(long double x)
+{
+	return atan2l(x, 1.0L);
+}
+
+export double cbrt(double x)
+{
+	return pow(x, 1.0 / 3.0);
+}
+
+export float cbrtf(float x)
+{
+	return powf(x, 1.0f / 3.0f);
+}
+
+export long double cbrtl(long double x)
+{
+	return powl(x, 1.0L / 3.0L);
+}
+
+export double ceil(double x)
+{
+	if (x == (double)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0) ? (double)((long)x + 1) : (double)((long)x);
+}
+
+export float ceilf(float x)
+{
+	if (x == (float)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0f) ? (float)((long)x + 1) : (float)((long)x);
+}
+
+export long double ceill(long double x)
+{
+	if (x == (long double)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0L) ? (long double)((long)x + 1) : (long double)((long)x);
+}
+
+export double copysign(double x, double y)
+{
+	return (y < 0.0) ? -fabs(x) : fabs(x);
+}
+
+export float copysignf(float x, float y)
+{
+	return (y < 0.0f) ? -fabsf(x) : fabsf(x);
+}
+
+export long double copysignl(long double x, long double y)
+{
+	return (y < 0.0L) ? -fabsl(x) : fabsl(x);
+}
+
+export double cos(double x)
+{
+	return sin(x + M_PI_2);
+}
+
+export float cosf(float x)
+{
+	return sinf(x + (float)M_PI_2);
+}
+
+export double cosh(double x)
+{
+	return (exp(x) + exp(-x)) / 2.0;
+}
+
+export float coshf(float x)
+{
+	return (expf(x) + expf(-x)) / 2.0f;
+}
+
+export long double coshl(long double x)
+{
+	return (expl(x) + expl(-x)) / 2.0L;
+}
+
+export long double cosl(long double x)
+{
+	return sinl(x + M_PI_2l);
+}
+
+export double erf(double x)
+{
+	double t = 1.0 / (1.0 + 0.5 * fabs(x));
+	double tau = t * exp(-x * x - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277)))))))));
+	return (x >= 0.0) ? 1.0 - tau : tau - 1.0;
+}
+
+export double erfc(double x)
+{
+	return 1.0 - erf(x);
+}
+
+export float erfcf(float x)
+{
+	return 1.0f - erff(x);
+}
+
+export long double erfcl(long double x)
+{
+	return 1.0L - erfl(x);
+}
+
+export float erff(float x)
+{
+	float t = 1.0f / (1.0f + 0.5f * fabsf(x));
+	float tau = t * expf(-x * x - 1.26551223f + t * (1.00002368f + t * (0.37409196f + t * (0.09678418f + t * (-0.18628806f + t * (0.27886807f + t * (-1.13520398f + t * (1.48851587f + t * (-0.82215223f + t * 0.17087277f)))))))));
+	return (x >= 0.0f) ? 1.0f - tau : tau - 1.0f;
+}
+
+export long double erfl(long double x)
+{
+	long double t = 1.0L / (1.0L + 0.5L * fabsl(x));
+	long double tau = t * expl(-x * x - 1.26551223L + t * (1.00002368L + t * (0.37409196L + t * (0.09678418L + t * (-0.18628806L + t * (0.27886807L + t * (-1.13520398L + t * (1.48851587L + t * (-0.82215223L + t * 0.17087277L)))))))));
+	return (x >= 0.0L) ? 1.0L - tau : tau - 1.0L;
+}
+
+export double exp(double x)
+{
+	if (x == 0.0)
+	{
+		return 1.0;
+	}
+	if (x < 0.0)
+	{
+		return 1.0 / exp(-x);
+	}
+	double result = 1.0;
+	double term = 1.0;
+	for (int i = 1; term > DBL_EPSILON || term < -DBL_EPSILON; i++)
+	{
+		term *= x / i;
+		result += term;
+	}
+	return result;
+}
+
+export double exp2(double x)
+{
+	if (x == 0.0)
+	{
+		return 1.0;
+	}
+	if (x < 0.0)
+	{
+		return 1.0 / exp2(-x);
+	}
+	double result = 1.0;
+	double term = 1.0;
+	for (int i = 1; term > DBL_EPSILON || term < -DBL_EPSILON; i++)
+	{
+		term *= x * log(2.0) / i;
+		result += term;
+	}
+	return result;
+}
+
+export float exp2f(float x)
+{
+	if (x == 0.0f)
+	{
+		return 1.0f;
+	}
+	if (x < 0.0f)
+	{
+		return 1.0f / exp2f(-x);
+	}
+	float result = 1.0f;
+	float term = 1.0f;
+	for (int i = 1; term > FLT_EPSILON || term < -FLT_EPSILON; i++)
+	{
+		term *= x * logf(2.0f) / i;
+		result += term;
+	}
+	return result;
+}
+
+export long double exp2l(long double x)
+{
+	return powl(2.0L, x);
+}
+
+export float expf(float x)
+{
+	return powf(M_E, x);
+}
+
+export long double expl(long double x)
+{
+	return powl(M_E, x);
+}
+
+export double expm1(double x)
+{
+	return exp(x) - 1.0;
+}
+
+export float expm1f(float x)
+{
+	return expf(x) - 1.0f;
+}
+
+export long double expm1l(long double x)
+{
+	return expl(x) - 1.0L;
+}
+
+export double fabs(double x)
+{
+	return (x < 0.0) ? -x : x;
+}
+
+export float fabsf(float x)
+{
+	return (x < 0.0f) ? -x : x;
+}
+
+export long double fabsl(long double x)
+{
+	return (x < 0.0L) ? -x : x;
+}
+
+export double fdim(double x, double y)
+{
+	return (x > y) ? x - y : 0.0;
+}
+
+export float fdimf(float x, float y)
+{
+	return (x > y) ? x - y : 0.0f;
+}
+
+export long double fdiml(long double x, long double y)
+{
+	return (x > y) ? x - y : 0.0L;
+}
+
+export double floor(double x)
+{
+	if (x == (double)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0) ? (double)((long)x) : (double)((long)x - 1);
+}
+
+export float floorf(float x)
+{
+	if (x == (float)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0f) ? (float)((long)x) : (float)((long)x - 1);
+}
+
+export long double floorl(long double x)
+{
+	if (x == (long double)((long)x))
+	{
+		return x;
+	}
+	return (x > 0.0L) ? (long double)((long)x) : (long double)((long)x - 1);
+}
+
+export double fma(double x, double y, double z)
+{
+	return x * y + z;
+}
+
+export float fmaf(float x, float y, float z)
+{
+	return x * y + z;
+}
+
+export long double fmal(long double x, long double y, long double z)
+{
+	return x * y + z;
+}
+
+export double fmax(double x, double y)
+{
+	return (x > y) ? x : y;
+}
+
+export float fmaxf(float x, float y)
+{
+	return (x > y) ? x : y;
+}
+
+export long double fmaxl(long double x, long double y)
+{
+	return (x > y) ? x : y;
+}
+
+export double fmin(double x, double y)
+{
+	return (x < y) ? x : y;
+}
+
+export float fminf(float x, float y)
+{
+	return (x < y) ? x : y;
+}
+
+export long double fminl(long double x, long double y)
+{
+	return (x < y) ? x : y;
+}
+
+export double fmod(double x, double y)
+{
+	if (y == 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return x - y * floor(x / y);
+}
+
+export float fmodf(float x, float y)
+{
+	if (y == 0.0f)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return x - y * floorf(x / y);
+}
+
+export long double fmodl(long double x, long double y)
+{
+	if (y == 0.0L)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return x - y * floorl(x / y);
+}
+
+export double frexp(double num, int *exp)
+{
+	if (num == 0.0)
+	{
+		*exp = 0;
+		return 0.0;
+	}
+	double abs_num = fabs(num);
+	*exp = (int)log2(abs_num) + 1;
+	double mantissa = abs_num / pow(2.0, *exp);
+	if (num < 0.0)
+	{
+		mantissa = -mantissa;
+	}
+	return mantissa;
+}
+
+export float frexpf(float num, int *exp)
+{
+	if (num == 0.0f)
+	{
+		*exp = 0;
+		return 0.0f;
+	}
+	float abs_num = fabsf(num);
+	*exp = (int)log2f(abs_num) + 1;
+	float mantissa = abs_num / powf(2.0f, *exp);
+	if (num < 0.0f)
+	{
+		mantissa = -mantissa;
+	}
+	return mantissa;
+}
+
+export long double frexpl(long double num, int *exp)
+{
+	if (num == 0.0L)
+	{
+		*exp = 0;
+		return 0.0L;
+	}
+	long double abs_num = fabsl(num);
+	*exp = (int)log2l(abs_num) + 1;
+	long double mantissa = abs_num / powl(2.0L, *exp);
+	if (num < 0.0L)
+	{
+		mantissa = -mantissa;
+	}
+	return mantissa;
+}
+
+export double hypot(double x, double y)
+{
+	return sqrt(x * x + y * y);
+}
+
+export float hypotf(float x, float y)
+{
+	return sqrtf(x * x + y * y);
+}
+
+export long double hypotl(long double x, long double y)
+{
+	return sqrtl(x * x + y * y);
+}
+
+export int ilogb(double x)
+{
+	if (x == 0.0)
+	{
+		errno = EDOM;
+		return FP_ILOGB0;
+	}
+	if (isnan(x))
+	{
+		errno = EDOM;
+		return FP_ILOGBNAN;
+	}
+	return (int)log2(fabs(x));
+}
+
+export int ilogbf(float x)
+{
+	if (x == 0.0f)
+	{
+		errno = EDOM;
+		return FP_ILOGB0;
+	}
+	if (isnan(x))
+	{
+		errno = EDOM;
+		return FP_ILOGBNAN;
+	}
+	return (int)log2f(fabsf(x));
+}
+
+export int ilogbl(long double x)
+{
+	if (x == 0.0L)
+	{
+		errno = EDOM;
+		return FP_ILOGB0;
+	}
+	if (isnan(x))
+	{
+		errno = EDOM;
+		return FP_ILOGBNAN;
+	}
+	return (int)log2l(fabsl(x));
+}
+export double j0(double x)
+{
+	double ax = fabs(x);
+	if (ax < 8.0)
+	{
+		double y = x * x;
+		double ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7 + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
+		double ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y * 1.0))));
+		return ans1 / ans2;
+	}
+	else
+	{
+		double z = 8.0 / ax;
+		double y = z * z;
+		double xx = ax - 0.785398164;
+		double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
+		double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));
+		return sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
+	}
+}
+
+export double j1(double x)
+{
+	double ax = fabs(x);
+	if (ax < 8.0)
+	{
+		double y = x * x;
+		double ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1 + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
+		double ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74 + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
+		return ans1 / ans2;
+	}
+	else
+	{
+		double z = 8.0 / ax;
+		double y = z * z;
+		double xx = ax - 2.356194491;
+		double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
+		double ans2 = 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
+		double ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
+		if (x < 0.0)
+			ans = -ans;
+		return ans;
+	}
+}
+
+#define BIGNO 1.0e10
+#define BIGNI 1.0e-10
+
+export double jn(int n, double x)
+{
+	if (n == 0)
+		return j0(x);
+	if (n == 1)
+		return j1(x);
+	if (x == 0.0)
+		return 0.0;
+	double ax = fabs(x);
+	double tox = 2.0 / ax;
+	double bj, bjm, bjp, sum;
+	int j, m;
+	if (ax > (double)n)
+	{
+		bjm = j0(ax);
+		bj = j1(ax);
+		for (j = 1; j < n; j++)
+		{
+			bjp = j * tox * bj - bjm;
+			bjm = bj;
+			bj = bjp;
+		}
+		sum = bj;
+	}
+	else
+	{
+#define ACC 40
+		m = 2 * ((n + (int)sqrt(ACC * n)) / 2);
+		bjp = sum = 0.0;
+		bj = 1.0;
+		for (j = m; j > 0; j--)
+		{
+			bjm = j * tox * bj - bjp;
+			bjp = bj;
+			bj = bjm;
+			if (fabs(bj) > BIGNO)
+			{
+				bj *= BIGNI;
+				bjp *= BIGNI;
+				sum *= BIGNI;
+			}
+			if (j == n)
+				sum = bjp;
+		}
+		sum *= j0(ax) / bj;
+	}
+	return (x < 0.0 && (n & 1)) ? -sum : sum;
+}
+
+export double ldexp(double x, int exp)
+{
+	return x * pow(2.0, exp);
+}
+
+export float ldexpf(float x, int exp)
+{
+	return x * powf(2.0f, exp);
+}
+
+export long double ldexpl(long double x, int exp)
+{
+	return x * powl(2.0L, exp);
+}
+export double lgamma(double x)
+{
+	if (x <= 0.0 && floor(x) == x)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = 0.0;
+	if (x < 0.5)
+	{
+		result = log(M_PI) - log(sin(M_PI * x)) - lgamma(1.0 - x);
+	}
+	else
+	{
+		x -= 1.0;
+		double a = 0.99999999999980993;
+		double coefficients[] = {
+			676.5203681218851,
+			-1259.1392167224028,
+			771.32342877765313,
+			-176.61502916214059,
+			12.507343278686905,
+			-0.13857109526572012,
+			9.9843695780195716e-6,
+			1.5056327351493116e-7};
+		for (int i = 0; i < 8; i++)
+		{
+			a += coefficients[i] / (x + i + 1);
+		}
+		double t = x + 7.5;
+		result = log(2.5066282746310007 * a) - t + log(t) * (x + 0.5);
+	}
+	return result;
+}
+
+export float lgammaf(float x)
+{
+	return (float)lgamma((double)x);
+}
+
+export long double lgammal(long double x)
+{
+	return (long double)lgamma((double)x);
+}
+
+export long long llrint(double x)
+{
+	return (long long)round(x);
+}
+
+export long long llrintf(float x)
+{
+	return (long long)roundf(x);
+}
+
+export long long llrintl(long double x)
+{
+	return (long long)roundl(x);
+}
+
+export long long llround(double x)
+{
+	return (long long)round(x);
+}
+
+export long long llroundf(float x)
+{
+	return (long long)roundf(x);
+}
+
+export long long llroundl(long double x)
+{
+	return (long long)roundl(x);
+}
+
+export double log(double x)
+{
+	if (x <= 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = 0.0;
+	double y = (x - 1) / (x + 1);
+	double y2 = y * y;
+	double term = y;
+	for (int i = 1; term > DBL_EPSILON || term < -DBL_EPSILON; i += 2)
+	{
+		result += term / i;
+		term *= y2;
+	}
+	return 2 * result;
+}
+
+export double log10(double x)
+{
+	return log(x) / log(10.0);
+}
+
+export float log10f(float x)
+{
+	return (float)log10((double)x);
+}
+
+export long double log10l(long double x)
+{
+	return (long double)log10((double)x);
+}
+
+export double log1p(double x)
+{
+	if (x <= -1.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = 0.0;
+	double term = x;
+	for (int i = 1; term > DBL_EPSILON || term < -DBL_EPSILON; i++)
+	{
+		result += term / i;
+		term *= -x;
+	}
+	return result;
+}
+
+export float log1pf(float x)
+{
+	return (float)log1p((double)x);
+}
+
+export long double log1pl(long double x)
+{
+	return (long double)log1p((double)x);
+}
+
+export double log2(double x)
+{
+	return log(x) / log(2.0);
+}
+
+export float log2f(float x)
+{
+	return (float)log2((double)x);
+}
+
+export long double log2l(long double x)
+{
+	return (long double)log2((double)x);
+}
+
+export double logb(double x)
+{
+	if (x == 0.0)
+	{
+		errno = ERANGE;
+		return -HUGE_VAL;
+	}
+	if (isnan(x))
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	if (isinf(x))
+	{
+		return HUGE_VAL;
+	}
+	int exp;
+	frexp(x, &exp);
+	return (double)(exp - 1);
+}
+
+export float logbf(float x)
+{
+	return (float)logb((double)x);
+}
+
+export long double logbl(long double x)
+{
+	return (long double)logb((double)x);
+}
+
+export float logf(float x)
+{
+	return (float)log((double)x);
+}
+
+export long double logl(long double x)
+{
+	return (long double)log((double)x);
+}
+
+export long lrint(double x)
+{
+	if (isnan(x) || isinf(x))
+	{
+		errno = EDOM;
+		return 0;
+	}
+	return (long)rint(x);
+}
+
+export long lrintf(float x)
+{
+	if (isnan(x) || isinf(x))
+	{
+		errno = EDOM;
+		return 0;
+	}
+	return (long)rintf(x);
+}
+
+export long lrintl(long double x)
+{
+	if (isnan(x) || isinf(x))
+	{
+		errno = EDOM;
+		return 0;
+	}
+	return (long)rintl(x);
+}
+
+export long lround(double x)
+{
+	return (long)round(x);
+}
+
+export long lroundf(float x)
+{
+	return (long)roundf(x);
+}
+
+export long lroundl(long double x)
+{
+	return (long)roundl(x);
+}
+
+export double modf(double x, double *iptr)
+{
+	*iptr = floor(x);
+	return x - *iptr;
+}
+
+export float modff(float value, float *iptr)
+{
+	*iptr = floorf(value);
+	return value - *iptr;
+}
+
+export long double modfl(long double x, long double *iptr)
+{
+	*iptr = floorl(x);
+	return x - *iptr;
+}
+
+export double nan(const char *tagp)
+{
+	return NAN;
+}
+
+export float nanf(const char *tagp)
+{
+	return NAN;
+}
+
+export long double nanl(const char *tagp)
+{
+	return NAN;
+}
+
+export double nearbyint(double x)
+{
+	return round(x);
+}
+
+export float nearbyintf(float x)
+{
+	return roundf(x);
+}
+
+export long double nearbyintl(long double x)
+{
+	return roundl(x);
+}
+
+export double nextafter(double x, double y)
+{
+	if (isnan(x) || isnan(y))
+	{
+		return NAN;
+	}
+	if (x == y)
+	{
+		return y;
+	}
+	if (x == 0.0)
+	{
+		return (y > 0.0) ? DBL_MIN : -DBL_MIN;
+	}
+	union
+	{
+		double d;
+		uint64_t u;
+	} ux = {x};
+	if ((x > 0.0) == (y > x))
+	{
+		ux.u++;
+	}
+	else
+	{
+		ux.u--;
+	}
+	return ux.d;
+}
+
+export float nextafterf(float x, float y)
+{
+	if (isnan(x) || isnan(y))
+	{
+		return NAN;
+	}
+	if (x == y)
+	{
+		return y;
+	}
+	if (x == 0.0f)
+	{
+		return (y > 0.0f) ? FLT_MIN : -FLT_MIN;
+	}
+	union
+	{
+		float f;
+		uint32_t u;
+	} ux = {x};
+	if ((x > 0.0f) == (y > x))
+	{
+		ux.u++;
+	}
+	else
+	{
+		ux.u--;
+	}
+	return ux.f;
+}
+
+export long double nextafterl(long double x, long double y)
+{
+	if (isnan(x) || isnan(y))
+	{
+		return NAN;
+	}
+	if (x == y)
+	{
+		return y;
+	}
+	if (x == 0.0L)
+	{
+		return (y > 0.0L) ? LDBL_MIN : -LDBL_MIN;
+	}
+	union
+	{
+		long double ld;
+		uint64_t u[2];
+	} ux = {x};
+	if ((x > 0.0L) == (y > x))
+	{
+		ux.u[0]++;
+		if (ux.u[0] == 0)
+		{
+			ux.u[1]++;
+		}
+	}
+	else
+	{
+		if (ux.u[0] == 0)
+		{
+			ux.u[1]--;
+		}
+		ux.u[0]--;
+	}
+	return ux.ld;
+}
+
+export double nexttoward(double x, long double y)
+{
+	return nextafterl(x, y);
+}
+
+export float nexttowardf(float x, long double y)
+{
+	return nextafterl(x, y);
+}
+
+export long double nexttowardl(long double x, long double y)
+{
+	return nextafterl(x, y);
+}
+
+export double pow(double x, double y)
+{
+	if (x == 0.0 && y == 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return exp(y * log(x));
+}
+
+export float powf(float x, float y)
+{
+	return (float)pow((double)x, (double)y);
+}
+
+export long double powl(long double x, long double y)
+{
+	return (long double)pow((double)x, (double)y);
+}
+
+export double remainder(double x, double y)
+{
+	if (y == 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	return x - y * round(x / y);
+}
+
+export float remainderf(float x, float y)
+{
+	return (float)remainder((double)x, (double)y);
+}
+
+export long double remainderl(long double x, long double y)
+{
+	return (long double)remainder((double)x, (double)y);
+}
+
+export double remquo(double x, double y, int *quo)
+{
+	if (y == 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = remainder(x, y);
+	*quo = (int)((x - result) / y);
+	return result;
+}
+
+export float remquof(float x, float y, int *quo)
+{
+	return (float)remquo((double)x, (double)y, quo);
+}
+
+export long double remquol(long double x, long double y, int *quo)
+{
+	return (long double)remquo((double)x, (double)y, quo);
+}
+
+export double rint(double x)
+{
+	int current_rounding_mode = fegetround();
+	double result;
+
+	switch (current_rounding_mode)
+	{
+	case FE_TONEAREST:
+		result = nearbyint(x);
+		break;
+	case FE_DOWNWARD:
+		result = floor(x);
+		break;
+	case FE_UPWARD:
+		result = ceil(x);
+		break;
+	case FE_TOWARDZERO:
+		result = trunc(x);
+		break;
+	default:
+		result = nearbyint(x);
+		break;
+	}
+
+	if (result != x)
+	{
+		feraiseexcept(FE_INEXACT);
+	}
+
+	return result;
+}
+
+export float rintf(float x)
+{
+	int current_rounding_mode = fegetround();
+	float result;
+
+	switch (current_rounding_mode)
+	{
+	case FE_TONEAREST:
+		result = nearbyintf(x);
+		break;
+	case FE_DOWNWARD:
+		result = floorf(x);
+		break;
+	case FE_UPWARD:
+		result = ceilf(x);
+		break;
+	case FE_TOWARDZERO:
+		result = truncf(x);
+		break;
+	default:
+		result = nearbyintf(x);
+		break;
+	}
+
+	if (result != x)
+	{
+		feraiseexcept(FE_INEXACT);
+	}
+
+	return result;
+}
+
+export long double rintl(long double x)
+{
+	int current_rounding_mode = fegetround();
+	long double result;
+
+	switch (current_rounding_mode)
+	{
+	case FE_TONEAREST:
+		result = nearbyintl(x);
+		break;
+	case FE_DOWNWARD:
+		result = floorl(x);
+		break;
+	case FE_UPWARD:
+		result = ceill(x);
+		break;
+	case FE_TOWARDZERO:
+		result = truncl(x);
+		break;
+	default:
+		result = nearbyintl(x);
+		break;
+	}
+
+	if (result != x)
+	{
+		feraiseexcept(FE_INEXACT);
+	}
+
+	return result;
+}
+
+export double round(double x)
+{
+	return (x < 0.0) ? ceil(x - 0.5) : floor(x + 0.5);
+}
+
+export float roundf(float x)
+{
+	return (x < 0.0f) ? ceilf(x - 0.5f) : floorf(x + 0.5f);
+}
+
+export long double roundl(long double x)
+{
+	return (x < 0.0L) ? ceill(x - 0.5L) : floorl(x + 0.5L);
+}
+
+export double scalb(double x, double n)
+{
+	if (isnan(x) || isnan(n))
+	{
+		return NAN;
+	}
+	if (n == 0.0)
+	{
+		return x;
+	}
+	if (isinf(x) && !isinf(n))
+	{
+		return x;
+	}
+	if (x == 0.0 && !isinf(n))
+	{
+		return x;
+	}
+	if (x == 0.0 && isinf(n))
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	if (isinf(x) && isinf(n))
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	if (isinf(n))
+	{
+		return (n > 0.0) ? HUGE_VAL : 0.0;
+	}
+	if (isinf(-n))
+	{
+		return (n < 0.0) ? -HUGE_VAL : 0.0;
+	}
+
+	int exp = (int)n;
+	double result = x;
+	if (exp > 0)
+	{
+		while (exp > 0)
+		{
+			result *= 2.0;
+			exp--;
+		}
+	}
+	else
+	{
+		while (exp < 0)
+		{
+			result /= 2.0;
+			exp++;
+		}
+	}
+	return result;
+}
+
+export double scalbln(double x, long n)
+{
+	return x * pow(2.0, (double)n);
+}
+
+export float scalblnf(float x, long n)
+{
+	return x * powf(2.0f, (float)n);
+}
+
+export long double scalblnl(long double x, long n)
+{
+	return x * powl(2.0L, (long double)n);
+}
+
+export double scalbn(double x, int n)
+{
+	return x * pow(2.0, (double)n);
+}
+
+export float scalbnf(float x, int n)
+{
+	return x * powf(2.0f, (float)n);
+}
+
+export long double scalbnl(long double x, int n)
+{
+	return x * powl(2.0L, (long double)n);
+}
+
+export double sin(double x)
+{
+	double result = 0.0;
+	double term = x;
+	double x2 = x * x;
+	for (int i = 1; term > DBL_EPSILON || term < -DBL_EPSILON; i += 2)
+	{
+		result += term;
+		term *= -x2 / ((i + 1) * (i + 2));
+	}
+	return result;
+}
+
+export float sinf(float x)
+{
+	return (float)sin((double)x);
+}
+
+export double sinh(double x)
+{
+	return (exp(x) - exp(-x)) / 2.0;
+}
+
+export float sinhf(float x)
+{
+	return (float)sinh((double)x);
+}
+
+export long double sinhl(long double x)
+{
+	return (long double)sinh((double)x);
+}
+
+export long double sinl(long double x)
+{
+	return (long double)sin((double)x);
+}
+
+export double sqrt(double x)
+{
+	if (x < 0.0)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = x;
+	double last_result;
+	do
+	{
+		last_result = result;
+		result = (result + x / result) / 2.0;
+	} while (fabs(result - last_result) > DBL_EPSILON);
+	return result;
+}
+
+export float sqrtf(float x)
+{
+	return (float)sqrt((double)x);
+}
+
+export long double sqrtl(long double x)
+{
+	return (long double)sqrt((double)x);
+}
+
+export double tan(double x)
+{
+	return sin(x) / cos(x);
+}
+
+export float tanf(float x)
+{
+	return (float)tan((double)x);
+}
+
+export double tanh(double x)
+{
+	return (exp(x) - exp(-x)) / (exp(x) + exp(-x));
+}
+
+export float tanhf(float x)
+{
+	return (float)tanh((double)x);
+}
+
+export long double tanhl(long double x)
+{
+	return (long double)tanh((double)x);
+}
+
+export long double tanl(long double x)
+{
+	return (long double)tan((double)x);
+}
+
+export double tgamma(double x)
+{
+	if (x <= 0.0 && floor(x) == x)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+	double result = 1.0;
+	if (x < 0.5)
+	{
+		return M_PI / (sin(M_PI * x) * tgamma(1.0 - x));
+	}
+	x -= 1.0;
+	double coefficients[] = {
+		1.000000000190015,
+		76.18009172947146,
+		-86.50532032941677,
+		24.01409824083091,
+		-1.231739572450155,
+		0.001208650973866179,
+		-0.000005395239384953};
+	double y = x + 5.5;
+	y -= (x + 0.5) * log(y);
+	for (int i = 0; i < 7; i++)
+	{
+		result += coefficients[i] / (x + i + 1);
+	}
+	return sqrt(2.0 * M_PI) * pow(x + 5.5, x + 0.5) * exp(-y) * result;
+}
+
+export float tgammaf(float x)
+{
+	return (float)tgamma((double)x);
+}
+
+export long double tgammal(long double x)
+{
+	return (long double)tgamma((double)x);
+}
+
+export double trunc(double x)
+{
+	return (x < 0.0) ? ceil(x) : floor(x);
+}
+
+export float truncf(float x)
+{
+	return (x < 0.0f) ? ceilf(x) : floorf(x);
+}
+
+export long double truncl(long double x)
+{
+	return (x < 0.0L) ? ceill(x) : floorl(x);
+}
+
+export double y0(double x)
+{
+	if (x < 8.0)
+	{
+		double y = x * x;
+		double ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6 + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
+		double ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438 + y * (47447.26470 + y * (226.1030244 + y * 1.0))));
+		return (ans1 / ans2) + 0.636619772 * j0(x) * log(x);
+	}
+	else
+	{
+		double z = 8.0 / x;
+		double y = z * z;
+		double xx = x - 0.785398164;
+		double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
+		double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));
+		return sqrt(0.636619772 / x) * (sin(xx) * ans1 + z * cos(xx) * ans2);
+	}
+}
+
+export double y1(double x)
+{
+	if (x < 8.0)
+	{
+		double y = x * x;
+		double ans1 = x * (-0.4900604943e13 + y * (0.1275274390e13 + y * (-0.5153438139e11 + y * (0.7349264551e9 + y * (-0.4237922726e7 + y * 0.8511937935e4)))));
+		double ans2 = 0.2499580570e14 + y * (0.4244419664e12 + y * (0.3733650367e10 + y * (0.2245904002e8 + y * (0.1020426050e6 + y * (0.3549632885e3 + y)))));
+		return (ans1 / ans2) + 0.636619772 * (j1(x) * log(x) - 1.0 / x);
+	}
+	else
+	{
+		double z = 8.0 / x;
+		double y = z * z;
+		double xx = x - 2.356194491;
+		double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
+		double ans2 = 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
+		return sqrt(0.636619772 / x) * (sin(xx) * ans1 + z * cos(xx) * ans2);
+	}
+}
+
+export double yn(int n, double x)
+{
+	if (n == 0)
+		return y0(x);
+	if (n == 1)
+		return y1(x);
+	if (x == 0.0)
+		return -HUGE_VAL;
+	if (x < 0.0)
+		return NAN;
+	double tox = 2.0 / x;
+	double by, bym, byp;
+	int j;
+	bym = y0(x);
+	by = y1(x);
+	for (j = 1; j < n; j++)
+	{
+		byp = j * tox * by - bym;
+		bym = by;
+		by = byp;
+	}
+	return by;
+}