+++ /dev/null
-#include "mathlib.qh"
-#if defined(CSQC)
-#elif defined(MENUQC)
-#elif defined(SVQC)
-#endif
-
-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;
-}
-bool isfinite(float x)
-{
- return !(isnan(x) || isinf(x));
-}
-bool isinf(float x)
-{
- return (x != 0) && (x + x == x);
-}
-bool isnan(float x)
-{
- float y;
- y = x;
- return (x != y);
-}
-bool isnormal(float x)
-{
- return isfinite(x);
-}
-bool 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 logn(float x, float base)
-{
- return log(x) / log(base);
-}
-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;
-}
-
-/**
- * Pythonic mod:
- * TODO: %% operator?
- *
- * 1 % 2 == 1
- * -1 % 2 == 1
- * 1 % -2 == -1
- * -1 % -2 == -1
- */
-float pymod(float x, float y)
-{
- return x - y * floor(x / 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);
-}
-
-vector cross(vector a, vector b)
-{
- return
- '1 0 0' * (a.y * b.z - a.z * b.y)
- + '0 1 0' * (a.z * b.x - a.x * b.z)
- + '0 0 1' * (a.x * b.y - a.y * b.x);
-}