module std::math;

// TODO Define these using quad precision.
const E = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466;
const LOG2E = 1.44269504088896340735992468100189214; // log2(e)
const LOG10E = 0.434294481903251827651128918916605082; // log10(e)
const LN2 = 0.693147180559945309417232121458176568; // ln(2)
const LN10 = 2.30258509299404568401799145468436421; // ln(10)
const PI = 3.14159265358979323846264338327950288419716939937510; // pi
const PI_2 = 1.57079632679489661923132169163975144; // pi / 2
const PI_4 = 0.785398163397448309615660845819875721; // pi / 4
const DIV_PI = 0.318309886183790671537767526745028724; // 1 / pi
const DIV_2_PI = 0.636619772367581343075535053490057448; // 2 / pi
const DIV_2_SQRTPI = 1.12837916709551257389615890312154517;   // 2/sqrt(pi)
const SQRT2 = 1.41421356237309504880168872420969808; // sqrt(2)
const DIV_1_SQRT2 = 0.707106781186547524400844362104849039; // 1 / sqrt(2)

const HALF_MAX = 6.5504e+4;
const HALF_MIN = 6.103515625e-5;
const HALF_DENORM_MIN = 5.9604644775390625e-8;
const HALF_DIG = 3;
const HALF_DEC_DIGITS = 5;
const HALF_MANT_DIG = 11;
const HALF_MAX_10_EXP = 4;
const HALF_MIN_10_EXP = -4;
const HALF_MAX_EXP = 16;
const HALF_MIN_EXP = -13;
const HALF_EPSILON = 9.765625e-4;

const FLOAT_MAX = 0x1.fffffep+127;
const FLOAT_MIN = 1.17549435e-38;
const FLOAT_DENORM_MIN = 1.40129846432481707092e-45;
const FLOAT_DIG = 6;
const FLOAT_DEC_DIGITS = 9;
const FLOAT_MANT_DIG = 24;
const FLOAT_MAX_10_EXP = 38;
const FLOAT_MIN_10_EXP = -37;
const FLOAT_MAX_EXP = 128;
const FLOAT_MIN_EXP = -125;
const FLOAT_EPSILON = 1.1920928955078125e-07;

const DOUBLE_MAX = 1.79769313486231570815e+308;
const DOUBLE_MIN = 2.2250738585072014e-308;
const DOUBLE_DENORM_MIN = 4.94065645841246544177e-324;
const DOUBLE_DIG = 15;
const DOUBLE_DEC_DIGITS = 17;
const DOUBLE_MANT_DIG = 53;
const DOUBLE_MAX_10_EXP = 308;
const DOUBLE_MIN_10_EXP = -307;
const DOUBLE_MAX_EXP = 1024;
const DOUBLE_MIN_EXP = -1021;
const DOUBLE_EPSILON = 2.22044604925031308085e-16;

const QUAD_MAX = 1.18973149535723176508575932662800702e+4932;
const QUAD_MIN = 3.36210314311209350626267781732175260e-4932;
const QUAD_DENORM_MIN = 6.47517511943802511092443895822764655e-4966;
const QUAD_DIG = 33;
const QUAD_DEC_DIGITS = 36;
const QUAD_MANT_DIG = 113;
const QUAD_MAX_10_EXP = 4932;
const QUAD_MIN_10_EXP = -4931;
const QUAD_MAX_EXP = 16384;
const QUAD_MIN_EXP = -16481;
const QUAD_EPSILON = 1.92592994438723585305597794258492732e-34;

private union DoubleLong
{
    double f;
    ulong i;
}

func double log10(double x)
{
    const double IVLN10HI = 4.34294481878168880939e-01; /* 0x3fdbcb7b, 0x15200000 */
    const double IVLN10LO = 2.50829467116452752298e-11; /* 0x3dbb9438, 0xca9aadd5 */
    const double LOG10_2HI = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
    const double LOG10_2LO = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
    const double LG1 = 6.666666666666735130e-01;  /* 3FE55555 55555593 */
    const double LG2 = 3.999999999940941908e-01;  /* 3FD99999 9997FA04 */
    const double LG3 = 2.857142874366239149e-01;  /* 3FD24924 94229359 */
    const double LG4 = 2.222219843214978396e-01;  /* 3FCC71C5 1D8E78AF */
    const double LG5 = 1.818357216161805012e-01;  /* 3FC74664 96CB03DE */
    const double LG6 = 1.531383769920937332e-01;  /* 3FC39A09 D078C69F */
    const double LG7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
    DoubleLong u = { .f = x };
	ulong hx = (uint)(u.i >> 32);
	int 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 = (uint)(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 = (ulong)(hx << 32) | (u.i & 0xffffffff);
	x = u.f;

	hx += 0x3ff00000 - 0x3fe6a09e;
	k += (int)(hx >> 20) - 0x3ff;
	hx = (hx & 0x000fffff) + 0x3fe6a09e;
	u.i = (ulong)(hx << 32) | (u.i & 0xffffffff);
	x = u.f;


	double f = x - 1.0;
	double hfsq = 0.5 * f * f;
	double s = f / (2.0+f);
	double z = s * s;
	double w = z * z;
	double t1 = w * (LG2 + w * (LG4 + w * LG6));
	double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
	double r = t2 + t1;

	/* See log2.c for details. */
	/* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
	double hi = f - hfsq;
	u.f = hi;
	// u.i &= (ulong)(-1) << 32;
	u.i &= 0xFFFFFFFF00000000;
	hi = u.f;
	double lo = f - hi - hfsq + s * (hfsq + r);

	/* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
	double val_hi = hi * IVLN10HI;
	double dk = k;
	double y = dk * LOG10_2HI;
	double val_lo = dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo*IVLN10HI;

	/*
	 * Extra precision in for adding y is not strictly needed
	 * since there is no very large cancellation near x = sqrt(2) or
	 * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
	 * with some parallelism and it reduces the error for many args.
	 */
	w = y + val_hi;
	val_lo += (y - w) + val_hi;
	val_hi = w;

	return val_lo + val_hi;
}

func double cos_limited(double x, double y)
{
    const double C1 = 4.16666666666666019037e-02; /* 0x3FA55555, 0x5555554C */
    const double C2 = -1.38888888888741095749e-03; /* 0xBF56C16C, 0x16C15177 */
    const double C3 = 2.48015872894767294178e-05; /* 0x3EFA01A0, 0x19CB1590 */
    const double C4 = -2.75573143513906633035e-07; /* 0xBE927E4F, 0x809C52AD */
    const double C5 = 2.08757232129817482790e-09; /* 0x3E21EE9E, 0xBDB4B1C4 */
    const double C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */

	double z = x * x;
	double w = z * z;
	double r = z * (C1+ z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
	double hz = 0.5 * z;
	w  = 1.0 - hz;
	return w + (((1.0 - w) - hz) + (z*r - x*y));
}

private func double sin_limited(double x, double y, bool iy)
{
    const double S1 = -1.66666666666666324348e-01;  // 0xBFC55555, 0x55555549
    const double S2 = 8.33333333332248946124e-03; // 0x3F811111, 0x1110F8A6
    const double S3 = -1.98412698298579493134e-04; // 0xBF2A01A0, 0x19C161D5
    const double S4 = 2.75573137070700676789e-06; // 0x3EC71DE3, 0x57B1FE7D
    const double S5 = -2.50507602534068634195e-08; // 0xBE5AE5E6, 0x8A2B9CEB
    const double S6 = 1.58969099521155010221e-10; // 0x3DE5D93A, 0x5ACFD57C

	double z = x * x;
	double w = z * z;
    double r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
	double v = z * x;
	if (!iy)
	{
		return x + v * (S1 + z * r);
	}
	else
	{
	    return x - ((z * (0.5 * y - v * r) - y) - v * S1);
	}
}
/*
public func double cos(double x)
{
	double[2] y;
	uint32_t ix;
	unsigned n;

	GET_HIGH_WORD(ix, x);
	ix &= 0x7fffffff;

	/* |x| ~< pi/4 */
	if (ix <= 0x3fe921fb)
	{
		if (ix < 0x3e46a09e)
		{
		    // |x| < 2**-27 * sqrt(2)
			/* raise inexact if x!=0 */
			FORCE_EVAL(x + 0x1p120f);
			return 1.0;
		}
		return cos_limited(x, 0);
	}

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

	/* argument reduction */
	n = __rem_pio2(x, y);
	switch (n&3)
	{
	    case 0: return  cos_limited(y[0], y[1]);
	    case 1: return -sin_limited(y[0], y[1], true);
	    case 2: return -cos_limited(y[0], y[1]);
	    default:
		    return  sin_limited(y[0], y[1], true);
	}
}
*/