/*
 * The content of this file is free knowledge.
 * No copyright applies. No warranty is provided.
 *
 * The following functions are available (aliases in parentheses):
 * General utility: max, min, sgn, abs, floor, ceiling
 * Arithmetic: factorial(fact), binomial(binom, bin), gcd, lcm, totient(phi)
 * Calculus: exp, log(ln), pow, sin, cos, tan, cosh, sinh, tanh
 *           atan, atan2, asin, acos, atanh, asinh, acosh
 * Constants (as functions): e(), pi()
 *
 * Approximation of inverse trigonometric functions and of pi is not good.
 *
 * For functions returning non-integer values, remember to set the
 * desired scaled value before calling the function.
 */

/* General utility functions */

define max(x, y) {
	if (x >= y) return x
	return y
}

define min(x, y) {
	if (x <= y) return x
	return y
}

define sgn(c) {
	if (x > 0) return 1
	if (x < 0) return -1
	return 0
}

define abs(x) {
	return x * sgn(x)
}

define floor(x) {
	auto s, y
	if (x < 0) return -ceiling(-x)
	s = scale
	scale = 0
	y = x / 1
	scale = s
	return y
}

define ceiling(x) {
	if (x < 0) return -floor(-x)
	if (floor(x) == x) return x
	return 1 + floor(x)
}

/* Arithmetic */

define factorial(n) {
	auto i, res, s
	s = scale
	scale = 0
	res = 1
	for (i = 1; i <= n; i++) res *= i
	scale = s
	return res
}

define fact(n) {
	return factorial(n)
}

define binomial(n, k) {
	auto i, j, told[], tnew[], s
	if (k < 0) return -1
	if (k > n) return 0
	s = scale
	scale = 0
	tnew[0] = 1
	for (i = 1; i <= n; i++) {
		for (j = 0; j <= i; j++) told[j] = tnew[j]
		for (j = 1; j <= i; j++) tnew[j] = told[j] + told[j-1]
	}
	scale = s
	return tnew[k]
}

define binom(n, k) {
	return binomial(n, k)
}

define bin(n, k) {
	return binomial(n, k)
}

define gcd(a, b) {
	auto s, aux
	s = scale
	scale = 0
	a /= 1
	b /= 1
	while (b != 0) {
		aux = a
		a = b
		b = aux % b
	}
	scale = s
	return a
}

define lcm(a, b) {
	if (a == 0 || b == 0) return 0
	return a * b / gcd(a, b)
}

define phi(n) {
	auto i, j, f[], r, s
	s = scale
	scale = 0
	j = 0
	r = n
	for (i = 2; i*i <= r; i++) {
		if (n % i == 0) {
			f[j] = i
			j += 1
			while (n % i == 0) n /= i
		}
	}
	for (i = 0; i < j; i++) r -= r / f[i]
	scale = s
	return r
}

define totient(n) {
	return phi(n)
}

/* Calculus */

define exp(x) {
	auto i, n, d, series
	scale += 10
	i = 0
	n = 1
	d = 1
	series = 0
	while (n/d != 0) {
		series += n/d
		i += 1
		n *= x
		d *= i
	}
	scale -= 10
	return series / 1
}

define e() {
	return exp(1)
}

/* Log: first take some square roots to make x smaller, then
        use Newton's method to solve e^y-x=0 */
define log(x) {
	auto m, s, t, n, d
	scale += 10
	while (x > 3) {
		m += 1
		x = sqrt(x)
	}
	t = x
	d = exp(t)
	n = d - x
	while (n/d != 0) {
		t -= n/d
		d = exp(t)
		n = d - x
	}
	scale -= 10
	return (t * 2^m) / 1
}

define ln(x) {
	return log(x)
}

define pow(a, b) {
	if (b != floor(b)) return (2^floor(b)) * exp(log(a)*(b-floor(b)))
	return a ^ b
}

define trig(x, ii, in, id) {
	auto i, n, d, series
	scale += 10
	i = ii
	n = in
	d = id
	while (n/d != 0) {
		series += n/d
		i += 2
		n *= -x*x
		d *= i*(i-1)
	}
	scale -= 10
	return series / 1
}

define sin(x) {
	return trig(x, 1, x, 1)
}

define cos(x) {
	return trig(x, 0, 1, 1)
}

define tan(x) {
	return sin(x) / cos(x)
}

define sinh(x) {
	return 0.5*(exp(x) - exp(-x))
}

define cosh(x) {
	return 0.5*(exp(x) + exp(-x))
}

define tanh(x) {
	return sinh(x) / cosh(x)
}

/* Atan: approximate integral of 1/(1+x^2) */
define atan(x) {
	auto i, n, f1, f2, a, res, s
	s = scale
	scale = 10
	n = 10^4
	res = 0
	for (i = 0; i < n; i++) {
		f1 = 1 + (x*i/n)^2
		f2 = 1 + (x*(i+1)/n)^2
		a = (f1 + f2) / (2 * f1 * f2)
		res += a * x / n
	}
	scale = s
	return res/1
}

define pi() {
	return 4 * atan(1)
}

define atan2(y, x) {
	if (x > 0) return atan(y/x)
	if (x < 0 && y >= 0) return atan(y/x) + pi()
	if (x < 0 && y < 0) return atan(y/x) - pi()
	if (x == 0 && y > 0) return pi() / 2
	return -pi() / 2
}

define asin(x) {
	return atan2(x, sqrt((1+x)*(1-x)))
}

define acos(x) {
	return atan2(sqrt((1+x)*(1-x)), x)
}

define atanh(x) {
	return 0.5 * log((1+x)/(1-x))
}

define asinh(x) {
	return log(x + sqrt(x^2 + 1))
}

define acosh(x) {
	return log(x + sqrt(x^2 - 1))
}