int fpclassify(float x) { if(isnan(x)) return FP_NAN; if(isinf(x)) return FP_INFINITE; if(x == 0) return FP_ZERO; return FP_NORMAL; } int isfinite(float x) { return !(isnan(x) || isinf(x)); } int isinf(float x) { return (x != 0) && (x + x == x); } int isnan(float x) { return !(x + x == x + x); } int isnormal(float x) { return isfinite(x); } int signbit(float x) { return (x < 0); } float acosh(float x) { return log(x + sqrt(x*x - 1)); } float asinh(float x) { return log(x + sqrt(x*x + 1)); } float atanh(float x) { return 0.5 * log((1+x) / (1-x)); } float cosh(float x) { return 0.5 * (exp(x) + exp(-x)); } float sinh(float x) { return 0.5 * (exp(x) - exp(-x)); } float tanh(float x) { return sinh(x) / cosh(x); } float exp(float x) { return pow(M_E, x); } float exp2(float x) { return pow(2, x); } float expm1(float x) { return exp(x) - 1; } vector frexp(float x) { vector v; v_z = 0; v_y = ilogb(x) + 1; v_x = x / exp2(v_y); return v; } int ilogb(float x) { return floor(log2(fabs(x))); } float ldexp(float x, int e) { return x * pow(2, e); } float log10(float x) { return log(x) * M_LOG10E; } float log1p(float x) { return log(x + 1); } float log2(float x) { return log(x) * M_LOG2E; } float logb(float x) { return floor(log2(fabs(x))); } vector modf(float f) { return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f); } float scalbn(float x, int n) { return x * pow(2, n); } float cbrt(float x) { return copysign(pow(fabs(x), 1.0/3.0), x); } float hypot(float x, float y) { return sqrt(x*x + y*y); } float erf(float x) { // approximation taken from wikipedia float y; y = x*x; return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x); } float erfc(float x) { return 1.0 - erf(x); } vector lgamma(float x) { // TODO improve accuracy if(!isfinite(x)) return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0'; if(x < 1 && x == floor(x)) return nan("gamma") * '1 1 1'; if(x < 0.1) { vector v; v = lgamma(1.0 - x); // reflection formula: // gamma(1-z) * gamma(z) = pi / sin(pi*z) // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z)) // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z) v_z = sin(M_PI * x); v_x = log(M_PI) - log(fabs(v_z)) - v_x; if(v_z < 0) v_y = -v_y; v_z = 0; return v; } if(x < 1.1) return lgamma(x + 1) - log(x) * '1 0 0'; x -= 1; return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0'; } float tgamma(float x) { vector v; v = lgamma(x); return exp(v_x) * v_y; } float nearbyint(float x) { return rint(x); } float trunc(float x) { return (x>=0) ? floor(x) : ceil(x); } float fmod(float x, float y) { return x - y * trunc(x / y); } float remainder(float x, float y) { return x - y * rint(x / y); } vector remquo(float x, float y) { vector v; v_z = 0; v_y = rint(x / y); v_x = x - y * v_y; return v; } float copysign(float x, float y) { return fabs(x) * ((y>0) ? 1 : -1); } float nan(string tag) { return sqrt(-1); } float nextafter(float x, float y) { // TODO very crude if(x == y) return nan("nextafter"); if(x > y) return -nextafter(-x, -y); // now we know that x < y // so we need the next number > x float d, a, b; d = max(fabs(x), 0.00000000000000000000001); a = x + d; do { d *= 0.5; b = a; a = x + d; } while(a != x); return b; } float nexttoward(float x, float y) { return nextafter(x, y); } float fdim(float x, float y) { return max(x-y, 0); } float fmax(float x, float y) { return max(x, y); } float fmin(float x, float y) { return min(x, y); } float fma(float x, float y, float z) { return x * y + z; } int isgreater(float x, float y) { return x > y; } int isgreaterequal(float x, float y) { return x >= y; } int isless(float x, float y) { return x < y; } int islessequal(float x, float y) { return x <= y; } int islessgreater(float x, float y) { return x < y || x > y; } int isunordered(float x, float y) { return !(x < y || x == y || x > y); }