#include <math.h>
#include <immintrin.h>

#include <bits/ensure.h>

#include <stdint.h>

#include <mlibc/debug.hpp>

// Taken from musl. See musl for the license/copyright!
#define FORCE_EVAL(x) do {                        \
	if (sizeof(x) == sizeof(float)) {         \
		volatile float __x;               \
		__x = (x);                        \
	} else if (sizeof(x) == sizeof(double)) { \
		volatile double __x;              \
		__x = (x);                        \
	} else {                                  \
		volatile long double __x;         \
		__x = (x);                        \
	}                                         \
} while(0)

namespace ieee754 {

struct SoftDouble {
	typedef uint64_t Bits;
	typedef uint64_t Mantissa;
	typedef int16_t Exp;

	static constexpr int kMantissaBits = 52;
	static constexpr int kExpBits = 11;
	static constexpr int kBias = 1023;
	
	// this exponent represents zeros (when mantissa = 0) and subnormals (when mantissa != 0)
	static constexpr Exp kSubExp = -kBias;
	// this exponent represents infinities (when mantissa = 0) and NaNs (when mantissa != 0)
	static constexpr Exp kInfExp = ((Exp(1) << kExpBits) - 1) - kBias;

	static constexpr Bits kMantissaMask = (Bits(1) << kMantissaBits) - 1;
	static constexpr Bits kExpMask = ((Bits(1) << kExpBits) - 1) << kMantissaBits;
	static constexpr Bits kSignMask = Bits(1) << (kMantissaBits + kExpBits);

	SoftDouble(bool negative, Mantissa mantissa, Exp exp)
	: negative(negative), mantissa(mantissa), exp(exp) {
//		mlibc::infoLogger.log() << "(" << (int)negative << ", " << (void *)mantissa
//				<< ", " << exp << ")" << frg::end_log;
		__ensure(mantissa < (Mantissa(1) << kMantissaBits));
		__ensure((exp + kBias) >= 0);
		__ensure((exp + kBias) < (Exp(1) << kExpBits));
	}

	const bool negative;
	const Mantissa mantissa;
	const Exp exp;
};

template<typename F>
using Bits = typename F::Bits;

template<typename F>
using Mantissa = typename F::Mantissa;

template<typename F>
using Exp = typename F::Exp;

template<typename F>
bool isZero(F x) {
	return x.exp == F::kSubExp && x.mantissa == 0;
}

template<typename F>
bool isFinite(F x) {
	return x.exp != F::kInfExp;
}

// --------------------------------------------------------
// Soft float operations
// --------------------------------------------------------

template<typename F>
F constZero(bool negative) {
	return F(negative, 0, F::kSubExp);
}

template<typename F>
F constOne(bool negative) {
	return F(negative, 0, 0);
}

template<typename F>
F floor(F x) {
	if(!isFinite(x) || isZero(x)) // TODO: need exception for the not-finite case?
		return x;
	
	if(x.exp > F::kMantissaBits)
		return x; // x is already integral
	
	if(x.exp < 0) {
		// TODO: raise inexact
		// return -1 or +0
		return x.negative ? constOne<F>(true) : constZero<F>(false);
	}
	
	Mantissa<F> mask = F::kMantissaMask >> x.exp;
	if(!(x.mantissa & mask))
		return x; // x is already integral
	
	// TODO: raise inexact
	Mantissa<F> integral_position = (Mantissa<F>(1) << F::kMantissaBits) >> x.exp;
	if(x.negative)
		return F(true, (x.mantissa + integral_position) & (~mask), x.exp);
	return F(false, x.mantissa & (~mask), x.exp);
}

template<typename F>
F ceil(F x) {
	if(!isFinite(x) || isZero(x)) // TODO: need exception for the not-finite case?
		return x;
	
	if(x.exp > F::kMantissaBits)
		return x; // x is already integral
	
	if(x.exp < 0) {
		// TODO: raise inexact
		// return -0 or +1
		return x.negative ? constZero<F>(true) : constOne<F>(false);
	}
	
	Mantissa<F> mask = F::kMantissaMask >> x.exp;
	if(!(x.mantissa & mask))
		return x; // x is already integral
	
	// TODO: raise inexact
	Mantissa<F> integral_position = (Mantissa<F>(1) << F::kMantissaBits) >> x.exp;
	if(x.negative)
		return F(true, x.mantissa & (~mask), x.exp);
	return F(false, (x.mantissa + integral_position) & (~mask), x.exp);
}

// --------------------------------------------------------
// Soft float <-> bit string conversion functions
// --------------------------------------------------------

template<typename F>
uint64_t compileBits(F soft) {
	auto bits = Bits<F>(soft.mantissa) | ((Bits<F>(soft.exp) + F::kBias) << soft.kMantissaBits);
	return soft.negative ? (F::kSignMask | bits) : bits;
}

SoftDouble extractBits(uint64_t bits) {
	return SoftDouble(bits & SoftDouble::kSignMask, bits & SoftDouble::kMantissaMask,
			((bits & SoftDouble::kExpMask) >> SoftDouble::kMantissaBits) - SoftDouble::kBias);
}

// --------------------------------------------------------
// Soft float -> native float conversion functions
// --------------------------------------------------------

union DoubleBits {
	double fp;
	uint64_t bits;
};

double compileNative(SoftDouble soft) {
	DoubleBits word;
	word.bits = compileBits(soft);
	return word.fp;
}

SoftDouble extractNative(double native) {
	DoubleBits word;
	word.fp = native;
	return extractBits(word.bits);
}

} // namespace ieee754

int __mlibc_fpclassify(double x) {
	return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x);
}
int __mlibc_fpclassifyf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
int __mlibc_fpclassifyl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double acos(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float acosf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double acosl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double asin(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float asinf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double asinl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double atan(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float atanf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double atanl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double atan2(double x, double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float atan2f(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double atan2l(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

// Taken from musl. See musl for the license/copyright!
float __sindf(double x) {
	/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
	static const double S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */
			S2 = 0x111110896efbb2.0p-59, /*  0.0083333293858894631756 */
			S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */
			S4 = 0x16cd878c3b46a7.0p-71; /*  0.0000027183114939898219064 */

	double r, s, w, z;

	/* Try to optimize for parallel evaluation as in __tandf.c. */
	z = x*x;
	w = z*z;
	r = S3 + z*S4;
	s = z*x;
	return (x + s*(S1 + z*S2)) + s*w*r;
}

// Taken from musl. See musl for the license/copyright!
float __cosdf(double x) {
	/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
	static const double C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
			C1 = 0x155553e1053a42.0p-57, /*  0.0416666233237390631894 */
			C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
			C3 = 0x199342e0ee5069.0p-68; /*  0.0000243904487962774090654 */

	double r, w, z;

	/* Try to optimize for parallel evaluation as in __tandf.c. */
	z = x*x;
	w = z*z;
	r = C2+z*C3;
	return ((1.0+z*C0) + w*C1) + (w*z)*r;
}

float __tandf(double x, int odd) {
	/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
	static const double T[] = {
		0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */
		0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */
		0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */
		0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */
		0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */
		0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */
	};

	double z,r,w,s,t,u;

	z = x*x;
	/*
	 * Split up the polynomial into small independent terms to give
	 * opportunities for parallel evaluation.  The chosen splitting is
	 * micro-optimized for Athlons (XP, X64).  It costs 2 multiplications
	 * relative to Horner's method on sequential machines.
	 *
	 * We add the small terms from lowest degree up for efficiency on
	 * non-sequential machines (the lowest degree terms tend to be ready
	 * earlier).  Apart from this, we don't care about order of
	 * operations, and don't need to to care since we have precision to
	 * spare.  However, the chosen splitting is good for accuracy too,
	 * and would give results as accurate as Horner's method if the
	 * small terms were added from highest degree down.
	 */
	r = T[4] + z*T[5];
	t = T[2] + z*T[3];
	w = z*z;
	s = z*x;
	u = T[0] + z*T[1];
	r = (x + s*u) + (s*w)*(t + w*r);
	return odd ? -1.0/r : r;
}

#define DBL_EPSILON 2.22044604925031308085e-16
#define EPS DBL_EPSILON

/* Get a 32 bit int from a float.  */
#define GET_FLOAT_WORD(w,d)                       \
do {                                              \
  union {float f; uint32_t i;} __u;               \
  __u.f = (d);                                    \
  (w) = __u.i;                                    \
} while (0)

/* Get the more significant 32 bit int from a double.  */
#define GET_HIGH_WORD(hi,d)                       \
do {                                              \
  union {double f; uint64_t i;} __u;              \
  __u.f = (d);                                    \
  (hi) = __u.i >> 32;                             \
} while (0)

/* Get the less significant 32 bit int from a double.  */
#define GET_LOW_WORD(lo,d)                        \
do {                                              \
  union {double f; uint64_t i;} __u;              \
  __u.f = (d);                                    \
  (lo) = (uint32_t)__u.i;                         \
} while (0)

// Taken from musl. See musl for the license/copyright!
int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
{
	static const int init_jk[] = {3,4,4,6}; /* initial value for jk */

	/*
	 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
	 *
	 *              integer array, contains the (24*i)-th to (24*i+23)-th
	 *              bit of 2/pi after binary point. The corresponding
	 *              floating value is
	 *
	 *                      ipio2[i] * 2^(-24(i+1)).
	 *
	 * NB: This table must have at least (e0-3)/24 + jk terms.
	 *     For quad precision (e0 <= 16360, jk = 6), this is 686.
	 */
	static const int32_t ipio2[] = {
	0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
	0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
	0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
	0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
	0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
	0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
	0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
	0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
	0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
	0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
	0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,

	#if LDBL_MAX_EXP > 1024
	0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
	0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
	0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
	0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
	0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
	0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
	0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
	0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
	0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
	0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
	0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
	0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
	0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
	0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
	0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
	0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
	0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
	0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
	0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
	0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
	0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
	0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
	0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
	0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
	0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
	0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
	0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
	0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
	0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
	0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
	0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
	0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
	0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
	0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
	0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
	0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
	0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
	0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
	0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
	0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
	0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
	0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
	0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
	0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
	0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
	0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
	0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
	0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
	0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
	0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
	0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
	0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
	0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
	0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
	0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
	0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
	0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
	0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
	0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
	0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
	0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
	0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
	0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
	0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
	0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
	0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
	0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
	0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
	0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
	0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
	0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
	0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
	0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
	0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
	0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
	0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
	0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
	0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
	0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
	0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
	0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
	0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
	0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
	0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
	0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
	0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
	0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
	0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
	0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
	0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
	0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
	0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
	0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
	0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
	0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
	0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
	0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
	0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
	0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
	0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
	0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
	0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
	0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
	0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
	#endif
	};

	static const double PIo2[] = {
	  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
	  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
	  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
	  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
	  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
	  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
	  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
	  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
	};

	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
	double z,fw,f[20],fq[20],q[20];

	/* initialize jk*/
	jk = init_jk[prec];
	jp = jk;

	/* determine jx,jv,q0, note that 3>q0 */
	jx = nx-1;
	jv = (e0-3)/24;  if(jv<0) jv=0;
	q0 = e0-24*(jv+1);

	/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
	j = jv-jx; m = jx+jk;
	for (i=0; i<=m; i++,j++)
		f[i] = j<0 ? 0.0 : (double)ipio2[j];

	/* compute q[0],q[1],...q[jk] */
	for (i=0; i<=jk; i++) {
		for (j=0,fw=0.0; j<=jx; j++)
			fw += x[j]*f[jx+i-j];
		q[i] = fw;
	}

	jz = jk;
recompute:
	/* distill q[] into iq[] reversingly */
	for (i=0,j=jz,z=q[jz]; j>0; i++,j--) {
		fw    = (double)(int32_t)(0x1p-24*z);
		iq[i] = (int32_t)(z - 0x1p24*fw);
		z     = q[j-1]+fw;
	}

	/* compute n */
	z  = scalbn(z,q0);       /* actual value of z */
	z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
	n  = (int32_t)z;
	z -= (double)n;
	ih = 0;
	if (q0 > 0) {  /* need iq[jz-1] to determine n */
		i  = iq[jz-1]>>(24-q0); n += i;
		iq[jz-1] -= i<<(24-q0);
		ih = iq[jz-1]>>(23-q0);
	}
	else if (q0 == 0) ih = iq[jz-1]>>23;
	else if (z >= 0.5) ih = 2;

	if (ih > 0) {  /* q > 0.5 */
		n += 1; carry = 0;
		for (i=0; i<jz; i++) {  /* compute 1-q */
			j = iq[i];
			if (carry == 0) {
				if (j != 0) {
					carry = 1;
					iq[i] = 0x1000000 - j;
				}
			} else
				iq[i] = 0xffffff - j;
		}
		if (q0 > 0) {  /* rare case: chance is 1 in 12 */
			switch(q0) {
			case 1:
				iq[jz-1] &= 0x7fffff; break;
			case 2:
				iq[jz-1] &= 0x3fffff; break;
			}
		}
		if (ih == 2) {
			z = 1.0 - z;
			if (carry != 0)
				z -= scalbn(1.0,q0);
		}
	}

	/* check if recomputation is needed */
	if (z == 0.0) {
		j = 0;
		for (i=jz-1; i>=jk; i--) j |= iq[i];
		if (j == 0) {  /* need recomputation */
			for (k=1; iq[jk-k]==0; k++);  /* k = no. of terms needed */

			for (i=jz+1; i<=jz+k; i++) {  /* add q[jz+1] to q[jz+k] */
				f[jx+i] = (double)ipio2[jv+i];
				for (j=0,fw=0.0; j<=jx; j++)
					fw += x[j]*f[jx+i-j];
				q[i] = fw;
			}
			jz += k;
			goto recompute;
		}
	}

	/* chop off zero terms */
	if (z == 0.0) {
		jz -= 1;
		q0 -= 24;
		while (iq[jz] == 0) {
			jz--;
			q0 -= 24;
		}
	} else { /* break z into 24-bit if necessary */
		z = scalbn(z,-q0);
		if (z >= 0x1p24) {
			fw = (double)(int32_t)(0x1p-24*z);
			iq[jz] = (int32_t)(z - 0x1p24*fw);
			jz += 1;
			q0 += 24;
			iq[jz] = (int32_t)fw;
		} else
			iq[jz] = (int32_t)z;
	}

	/* convert integer "bit" chunk to floating-point value */
	fw = scalbn(1.0,q0);
	for (i=jz; i>=0; i--) {
		q[i] = fw*(double)iq[i];
		fw *= 0x1p-24;
	}

	/* compute PIo2[0,...,jp]*q[jz,...,0] */
	for(i=jz; i>=0; i--) {
		for (fw=0.0,k=0; k<=jp && k<=jz-i; k++)
			fw += PIo2[k]*q[i+k];
		fq[jz-i] = fw;
	}

	/* compress fq[] into y[] */
	switch(prec) {
	case 0:
		fw = 0.0;
		for (i=jz; i>=0; i--)
			fw += fq[i];
		y[0] = ih==0 ? fw : -fw;
		break;
	case 1:
	case 2:
		fw = 0.0;
		for (i=jz; i>=0; i--)
			fw += fq[i];
		// TODO: drop excess precision here once double_t is used
		fw = (double)fw;
		y[0] = ih==0 ? fw : -fw;
		fw = fq[0]-fw;
		for (i=1; i<=jz; i++)
			fw += fq[i];
		y[1] = ih==0 ? fw : -fw;
		break;
	case 3:  /* painful */
		for (i=jz; i>0; i--) {
			fw      = fq[i-1]+fq[i];
			fq[i]  += fq[i-1]-fw;
			fq[i-1] = fw;
		}
		for (i=jz; i>1; i--) {
			fw      = fq[i-1]+fq[i];
			fq[i]  += fq[i-1]-fw;
			fq[i-1] = fw;
		}
		for (fw=0.0,i=jz; i>=2; i--)
			fw += fq[i];
		if (ih==0) {
			y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
		} else {
			y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
		}
	}
	return n&7;
}

int __rem_pio2f(float x, double *y) {
	/*
	 * invpio2:  53 bits of 2/pi
	 * pio2_1:   first 25 bits of pi/2
	 * pio2_1t:  pi/2 - pio2_1
	 */
	static const double toint = 1.5/EPS,
			invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
			pio2_1  = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
			pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */

	union {float f; uint32_t i;} u = {x};
	double tx[1],ty[1];
	double fn;
	uint32_t ix;
	int n, sign, e0;

	ix = u.i & 0x7fffffff;
	/* 25+53 bit pi is good enough for medium size */
	if (ix < 0x4dc90fdb) {  /* |x| ~< 2^28*(pi/2), medium size */
		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
		fn = (double)x*invpio2 + toint - toint;
		n  = (int32_t)fn;
		*y = x - fn*pio2_1 - fn*pio2_1t;
		return n;
	}
	if(ix>=0x7f800000) {  /* x is inf or NaN */
		*y = x-x;
		return 0;
	}
	/* scale x into [2^23, 2^24-1] */
	sign = u.i>>31;
	e0 = (ix>>23) - (0x7f+23);  /* e0 = ilogb(|x|)-23, positive */
	u.i = ix - (e0<<23);
	tx[0] = u.f;
	n  =  __rem_pio2_large(tx,ty,e0,1,0);
	if (sign) {
		*y = -ty[0];
		return -n;
	}
	*y = ty[0];
	return n;
}

double cos(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
// Taken from musl. See musl for the license/copyright!
float cosf(float x) {
	/* Small multiples of pi/2 rounded to double precision. */
	static const double c1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
			c2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
			c3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
			c4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */

	double y;
	uint32_t ix;
	unsigned n, sign;

	GET_FLOAT_WORD(ix, x);
	sign = ix >> 31;
	ix &= 0x7fffffff;

	if (ix <= 0x3f490fda) {  /* |x| ~<= pi/4 */
		if (ix < 0x39800000) {  /* |x| < 2**-12 */
			/* raise inexact if x != 0 */
			FORCE_EVAL(x + 0x1p120f);
			return 1.0f;
		}
		return __cosdf(x);
	}
	if (ix <= 0x407b53d1) {  /* |x| ~<= 5*pi/4 */
		if (ix > 0x4016cbe3)  /* |x|  ~> 3*pi/4 */
			return -__cosdf(sign ? x+c2pio2 : x-c2pio2);
		else {
			if (sign)
				return __sindf(x + c1pio2);
			else
				return __sindf(c1pio2 - x);
		}
	}
	if (ix <= 0x40e231d5) {  /* |x| ~<= 9*pi/4 */
		if (ix > 0x40afeddf)  /* |x| ~> 7*pi/4 */
			return __cosdf(sign ? x+c4pio2 : x-c4pio2);
		else {
			if (sign)
				return __sindf(-x - c3pio2);
			else
				return __sindf(x - c3pio2);
		}
	}

	/* cos(Inf or NaN) is NaN */
	if (ix >= 0x7f800000)
		return x-x;

	/* general argument reduction needed */
	n = __rem_pio2f(x,&y);
	switch (n&3) {
	case 0: return  __cosdf(y);
	case 1: return  __sindf(-y);
	case 2: return -__cosdf(y);
	default:
		return  __sindf(y);
	}
}
long double cosl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double sin(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
// Taken from musl. See musl for the license/copyright!
float sinf(float x) {
	/* Small multiples of pi/2 rounded to double precision. */
	static const double s1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
		s2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
		s3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
		s4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */

	double y;
	uint32_t ix;
	int n, sign;

	GET_FLOAT_WORD(ix, x);
	sign = ix >> 31;
	ix &= 0x7fffffff;

	if (ix <= 0x3f490fda) {  /* |x| ~<= pi/4 */
		if (ix < 0x39800000) {  /* |x| < 2**-12 */
			/* raise inexact if x!=0 and underflow if subnormal */
			FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f);
			return x;
		}
		return __sindf(x);
	}
	if (ix <= 0x407b53d1) {  /* |x| ~<= 5*pi/4 */
		if (ix <= 0x4016cbe3) {  /* |x| ~<= 3pi/4 */
			if (sign)
				return -__cosdf(x + s1pio2);
			else
				return __cosdf(x - s1pio2);
		}
		return __sindf(sign ? -(x + s2pio2) : -(x - s2pio2));
	}
	if (ix <= 0x40e231d5) {  /* |x| ~<= 9*pi/4 */
		if (ix <= 0x40afeddf) {  /* |x| ~<= 7*pi/4 */
			if (sign)
				return __cosdf(x + s3pio2);
			else
				return -__cosdf(x - s3pio2);
		}
		return __sindf(sign ? x + s4pio2 : x - s4pio2);
	}

	/* sin(Inf or NaN) is NaN */
	if (ix >= 0x7f800000)
		return x - x;

	/* general argument reduction needed */
	n = __rem_pio2f(x, &y);
	switch (n&3) {
	case 0: return  __sindf(y);
	case 1: return  __cosdf(y);
	case 2: return  __sindf(-y);
	default:
		return -__cosdf(y);
	}
}
long double sinl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double tan(double x) {
	mlibc::infoLogger() << "mlibc: tan() is not precise" << frg::endlog;
	return tanf(x);
}
// Taken from musl. See musl for the license/copyright!
float tanf(float x) {
	/* Small multiples of pi/2 rounded to double precision. */
	static const double t1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
			t2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
			t3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
			t4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */

	double y;
	uint32_t ix;
	unsigned n, sign;

	GET_FLOAT_WORD(ix, x);
	sign = ix >> 31;
	ix &= 0x7fffffff;

	if (ix <= 0x3f490fda) {  /* |x| ~<= pi/4 */
		if (ix < 0x39800000) {  /* |x| < 2**-12 */
			/* raise inexact if x!=0 and underflow if subnormal */
			FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f);
			return x;
		}
		return __tandf(x, 0);
	}
	if (ix <= 0x407b53d1) {  /* |x| ~<= 5*pi/4 */
		if (ix <= 0x4016cbe3)  /* |x| ~<= 3pi/4 */
			return __tandf((sign ? x+t1pio2 : x-t1pio2), 1);
		else
			return __tandf((sign ? x+t2pio2 : x-t2pio2), 0);
	}
	if (ix <= 0x40e231d5) {  /* |x| ~<= 9*pi/4 */
		if (ix <= 0x40afeddf)  /* |x| ~<= 7*pi/4 */
			return __tandf((sign ? x+t3pio2 : x-t3pio2), 1);
		else
			return __tandf((sign ? x+t4pio2 : x-t4pio2), 0);
	}

	/* tan(Inf or NaN) is NaN */
	if (ix >= 0x7f800000)
		return x - x;

	/* argument reduction */
	n = __rem_pio2f(x, &y);
	return __tandf(y, n&1);
}
long double tanl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double acosh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float acoshf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double acoshl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double asinh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float asinhf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double asinhl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double atanh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float atanhf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double atanhl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double cosh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float coshf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double coshl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double sinh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float sinhf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double sinhl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double tanh(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float tanhf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double tanhl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double exp(double x) {
	static const double half[2] = {0.5,-0.5},
			ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
			ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
			invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
			P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
			P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
			P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
			P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
			P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */

	double hi, lo, c, xx, y;
	int k, sign;
	uint32_t hx;

	GET_HIGH_WORD(hx, x);
	sign = hx>>31;
	hx &= 0x7fffffff;  /* high word of |x| */

	/* special cases */
	if (hx >= 0x4086232b) {  /* if |x| >= 708.39... */
		if (isnan(x))
			return x;
		if (x > 709.782712893383973096) {
			/* overflow if x!=inf */
			x *= 0x1p1023;
			return x;
		}
		if (x < -708.39641853226410622) {
			/* underflow if x!=-inf */
			FORCE_EVAL((float)(-0x1p-149/x));
			if (x < -745.13321910194110842)
				return 0;
		}
	}

	/* argument reduction */
	if (hx > 0x3fd62e42) {  /* if |x| > 0.5 ln2 */
		if (hx >= 0x3ff0a2b2)  /* if |x| >= 1.5 ln2 */
			k = (int)(invln2*x + half[sign]);
		else
			k = 1 - sign - sign;
		hi = x - k*ln2hi;  /* k*ln2hi is exact here */
		lo = k*ln2lo;
		x = hi - lo;
	} else if (hx > 0x3e300000)  {  /* if |x| > 2**-28 */
		k = 0;
		hi = x;
		lo = 0;
	} else {
		/* inexact if x!=0 */
		FORCE_EVAL(0x1p1023 + x);
		return 1 + x;
	}

	/* x is now in primary range */
	xx = x*x;
	c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
	y = 1 + (x*c/(2-c) - lo + hi);
	if (k == 0)
		return y;
	return scalbn(y, k);
}
float expf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double expl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double exp2(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
// Taken from musl. See musl for the license/copyright!
float exp2f(float x) {
	constexpr int TBLSIZE = 16;

	constexpr float redux = 0x1.8p23f / TBLSIZE;
	constexpr float P1 = 0x1.62e430p-1f;
	constexpr float P2 = 0x1.ebfbe0p-3f;
	constexpr float P3 = 0x1.c6b348p-5f;
	constexpr float P4 = 0x1.3b2c9cp-7f;

	constexpr double exp2ft[TBLSIZE] = {
		0x1.6a09e667f3bcdp-1,
		0x1.7a11473eb0187p-1,
		0x1.8ace5422aa0dbp-1,
		0x1.9c49182a3f090p-1,
		0x1.ae89f995ad3adp-1,
		0x1.c199bdd85529cp-1,
		0x1.d5818dcfba487p-1,
		0x1.ea4afa2a490dap-1,
		0x1.0000000000000p+0,
		0x1.0b5586cf9890fp+0,
		0x1.172b83c7d517bp+0,
		0x1.2387a6e756238p+0,
		0x1.306fe0a31b715p+0,
		0x1.3dea64c123422p+0,
		0x1.4bfdad5362a27p+0,
		0x1.5ab07dd485429p+0,
	};

	double t, r, z;
	union {float f; uint32_t i;} u = {x};
	union {double f; uint64_t i;} uk;
	uint32_t ix, i0, k;

	/* Filter out exceptional cases. */
	ix = u.i & 0x7fffffff;
	if (ix > 0x42fc0000) {  /* |x| > 126 */
		if (ix > 0x7f800000) /* NaN */
			return x;
		if (u.i >= 0x43000000 && u.i < 0x80000000) {  /* x >= 128 */
			x *= 0x1p127f;
			return x;
		}
		if (u.i >= 0x80000000) {  /* x < -126 */
			if (u.i >= 0xc3160000 || (u.i & 0x0000ffff))
				FORCE_EVAL(-0x1p-149f/x);
			if (u.i >= 0xc3160000)  /* x <= -150 */
				return 0;
		}
	} else if (ix <= 0x33000000) {  /* |x| <= 0x1p-25 */
		return 1.0f + x;
	}

	/* Reduce x, computing z, i0, and k. */
	u.f = x + redux;
	i0 = u.i;
	i0 += TBLSIZE / 2;
	k = i0 / TBLSIZE;
	uk.i = (uint64_t)(0x3ff + k)<<52;
	i0 &= TBLSIZE - 1;
	u.f -= redux;
	z = x - u.f;
	/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
	r = exp2ft[i0];
	t = r * z;
	r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4);

	/* Scale by 2**k */
	return r * uk.f;
}
long double exp2l(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double expm1(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float expm1f(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double expm1l(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double frexp(double x, int *power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float frexpf(float x, int *power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double frexpl(long double x, int *power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double ilogb(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float ilogbf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double ilogbl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double ldexp(double x, int power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float ldexpf(float x, int power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double ldexpl(long double x, int power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double log(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float logf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double logl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double log10(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float log10f(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double log10l(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double log1p(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float log1pf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double log1pl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

// Taken from musl. See musl for the license/copyright!
double log2(double x) {
	static const double
		ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
		ivln2lo = 1.67517131648865118353e-10, /* 0x3de705fc, 0x2eefa200 */
		Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
		Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
		Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
		Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
		Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
		Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
		Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */

	union {double f; uint64_t i;} u = {x};
	double hfsq,f,s,z,R,w,t1,t2,y,hi,lo,val_hi,val_lo;
	uint32_t hx;
	int k;

	hx = u.i>>32;
	k = 0;
	if (hx < 0x00100000 || hx>>31) {
		if (u.i<<1 == 0)
			return -1/(x*x);  /* log(+-0)=-inf */
		if (hx>>31)
			return (x-x)/0.0; /* log(-#) = NaN */
		/* subnormal number, scale x up */
		k -= 54;
		x *= 0x1p54;
		u.f = x;
		hx = u.i>>32;
	} else if (hx >= 0x7ff00000) {
		return x;
	} else if (hx == 0x3ff00000 && u.i<<32 == 0)
		return 0;

	/* reduce x into [sqrt(2)/2, sqrt(2)] */
	hx += 0x3ff00000 - 0x3fe6a09e;
	k += (int)(hx>>20) - 0x3ff;
	hx = (hx&0x000fffff) + 0x3fe6a09e;
	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff);
	x = u.f;

	f = x - 1.0;
	hfsq = 0.5*f*f;
	s = f/(2.0+f);
	z = s*s;
	w = z*z;
	t1 = w*(Lg2+w*(Lg4+w*Lg6));
	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
	R = t2 + t1;

	/*
	 * f-hfsq must (for args near 1) be evaluated in extra precision
	 * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
	 * This is fairly efficient since f-hfsq only depends on f, so can
	 * be evaluated in parallel with R.  Not combining hfsq with R also
	 * keeps R small (though not as small as a true `lo' term would be),
	 * so that extra precision is not needed for terms involving R.
	 *
	 * Compiler bugs involving extra precision used to break Dekker's
	 * theorem for spitting f-hfsq as hi+lo, unless double_t was used
	 * or the multi-precision calculations were avoided when double_t
	 * has extra precision.  These problems are now automatically
	 * avoided as a side effect of the optimization of combining the
	 * Dekker splitting step with the clear-low-bits step.
	 *
	 * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
	 * precision to avoid a very large cancellation when x is very near
	 * these values.  Unlike the above cancellations, this problem is
	 * specific to base 2.  It is strange that adding +-1 is so much
	 * harder than adding +-ln2 or +-log10_2.
	 *
	 * This uses Dekker's theorem to normalize y+val_hi, so the
	 * compiler bugs are back in some configurations, sigh.  And I
	 * don't want to used double_t to avoid them, since that gives a
	 * pessimization and the support for avoiding the pessimization
	 * is not yet available.
	 *
	 * The multi-precision calculations for the multiplications are
	 * routine.
	 */

	/* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
	hi = f - hfsq;
	u.f = hi;
	u.i &= (uint64_t)-1<<32;
	hi = u.f;
	lo = f - hi - hfsq + s*(hfsq+R);

	val_hi = hi*ivln2hi;
	val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;

	/* spadd(val_hi, val_lo, y), except for not using double_t: */
	y = k;
	w = y + val_hi;
	val_lo += (y - w) + val_hi;
	val_hi = w;

	return val_lo + val_hi;
}
float log2f(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double log2l(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double logb(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float logbf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double logbl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double modf(double x, double *integral) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float modff(float x, float *integral) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double modfl(long double x, long double *integral) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double scalbn(double x, int n) {
	union {double f; uint64_t i;} u;
	double y = x;

	if (n > 1023) {
		y *= 0x1p1023;
		n -= 1023;
		if (n > 1023) {
			y *= 0x1p1023;
			n -= 1023;
			if (n > 1023)
				n = 1023;
		}
	} else if (n < -1022) {
		/* make sure final n < -53 to avoid double
		   rounding in the subnormal range */
		y *= 0x1p-1022 * 0x1p53;
		n += 1022 - 53;
		if (n < -1022) {
			y *= 0x1p-1022 * 0x1p53;
			n += 1022 - 53;
			if (n < -1022)
				n = -1022;
		}
	}
	u.i = (uint64_t)(0x3ff+n)<<52;
	x = y * u.f;
	return x;
}
float scalbnf(float x, int power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double scalbnl(long double x, int power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double scalbln(double x, long power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float scalblnf(float x, long power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double scalblnl(long double x, long power) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double cbrt(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float cbrtf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double cbrtl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double fabs(double x) {
	return signbit(x) ? -x : x;
}
float fabsf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double fabsl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double hypot(double x, double y) {
	__ensure(isfinite(x));
	__ensure(isfinite(y));
	// TODO: fix exception handling
	double u = fabs(x);
	double v = fabs(y);
	if(u > v)
		return u * sqrt(1 + (v / u) * (v / u));
	return v * sqrt(1 + (u / v) * (u / v));
}
float hypotf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double hypotl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double pow(double x, double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float powf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double powl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double sqrt(double x) {
	auto sse_x = _mm_set_sd(x);
	return _mm_cvtsd_f64(_mm_sqrt_sd(sse_x, sse_x));
}
float sqrtf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double sqrtl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double erf(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float erff(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double erfl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double erfc(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float erfcf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double erfcl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double lgamma(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float lgammaf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double lgammal(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double tgamma(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float tgammaf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double tgammal(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double ceil(double x) {
	auto soft_x = ieee754::extractNative(x);
	auto result = ieee754::ceil(soft_x);
	return ieee754::compileNative(result);
}
float ceilf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double ceill(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double floor(double x) {
	auto soft_x = ieee754::extractNative(x);
	auto result = ieee754::floor(soft_x);
	return ieee754::compileNative(result);
}
float floorf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double floorl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double nearbyint(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float nearbyintf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double nearbyintl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double rint(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float rintf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double rintl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

long lrint(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long lrintf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long lrintl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

long long llrint(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long long llrintf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long long llrintl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double round(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float roundf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double roundl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

long lround(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long lroundf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long lroundl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

long long llround(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long long llroundf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long long llroundl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double trunc(double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float truncf(float x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double truncl(long double x) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double fmod(double x, double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float fmodf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double fmodl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double remainder(double x, double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float remainderf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double remainderl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double remquo(double x, double y, int *quotient) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float remquof(float x, float y, int *quotient) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double remquol(long double x, long double y, int *quotient) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double copysign(double x, double sign) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float copysignf(float x, float sign) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double copysignl(long double x, long double sign) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double nan(const char *tag) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float nanf(const char *tag) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double nanl(const char *tag) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double nextafter(double x, double dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float nextafterf(float x, float dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double nextafterl(long double x, long double dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double nexttoward(double x, long double dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float nexttowardf(float x, long double dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double nexttowardl(long double x, long double dir) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double fdim(double x, double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
float fdimf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double fdiml(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double fmax(double x, double y) {
	__ensure(isfinite(x) && isfinite(y));
	return x < y ? y : x;
}
float fmaxf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double fmaxl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

double fmin(double x, double y) {
	__ensure(isfinite(x) && isfinite(y));
	return x < y ? x : y;
}
float fminf(float x, float y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}
long double fminl(long double x, long double y) {
	__ensure(!"Not implemented");
	__builtin_unreachable();
}

//gnu extension

void sincos(double x, double *sx, double *cx) {
	mlibc::infoLogger() << "mlibc: sincos() is not precise" << frg::endlog;
	float sxf;
	float cxf;
	sincosf(x, &sxf, &cxf);
	*sx = sxf;
	*cx = cxf;
}

void sincosf(float x, float *sx, float *cx) {
	// This is a lazy implementation.
	__ensure(sx);
	__ensure(cx);
	*sx = sinf(x);
	*cx = cosf(x);
}
void sincosl(long double, long double *, long double *) {
	__ensure(!"sincosl() not implemented");
	__builtin_unreachable();
}

double exp10(double) {
	__ensure(!"exp10() not implemented");
	__builtin_unreachable();
}
float exp10f(float) {
	__ensure(!"exp10f() not implemented");
	__builtin_unreachable();
}
long double exp10l(long double) {
	__ensure(!"exp10l() not implemented");
	__builtin_unreachable();
}

double pow10(double) {
	__ensure(!"pow10() not implemented");
	__builtin_unreachable();
}
float pow10f(float) {
	__ensure(!"pow10f() not implemented");
	__builtin_unreachable();
}
long double pow10l(long double) {
	__ensure(!"pow10l() not implemented");
	__builtin_unreachable();
}