qcsrc/warpzonelib/mathlib.qc

   1 int fpclassify(float x)
   2 {
   3         if(isnan(x))
   4                 return FP_NAN;
   5         if(isinf(x))
   6                 return FP_INFINITE;
   7         if(x == 0)
   8                 return FP_ZERO;
   9         return FP_NORMAL;
  10 }
  11 int isfinite(float x)
  12 {
  13         return !(isnan(x) || isinf(x));
  14 }
  15 int isinf(float x)
  16 {
  17         return (x != 0) && (x + x == x);
  18 }
  19 int isnan(float x)
  20 {
  21         return !(x + x == x + x);
  22 }
  23 int isnormal(float x)
  24 {
  25         return isfinite(x);
  26 }
  27 int signbit(float x)
  28 {
  29         return (x < 0);
  30 }
  31
  32 float acosh(float x)
  33 {
  34         return log(x + sqrt(x*x - 1));
  35 }
  36 float asinh(float x)
  37 {
  38         return log(x + sqrt(x*x + 1));
  39 }
  40 float atanh(float x)
  41 {
  42         return 0.5 * log((1+x) / (1-x));
  43 }
  44 float cosh(float x)
  45 {
  46         return 0.5 * (exp(x) + exp(-x));
  47 }
  48 float sinh(float x)
  49 {
  50         return 0.5 * (exp(x) - exp(-x));
  51 }
  52 float tanh(float x)
  53 {
  54         return sinh(x) / cosh(x);
  55 }
  56
  57 float exp(float x)
  58 {
  59         return pow(M_E, x);
  60 }
  61 float exp2(float x)
  62 {
  63         return pow(2, x);
  64 }
  65 float expm1(float x)
  66 {
  67         return exp(x) - 1;
  68 }
  69
  70 vector frexp(float x)
  71 {
  72         vector v;
  73         v_z = 0;
  74         v_y = ilogb(x) + 1;
  75         v_x = x / exp2(v_y);
  76         return v;
  77 }
  78 int ilogb(float x)
  79 {
  80         return floor(log2(fabs(x)));
  81 }
  82 float ldexp(float x, int e)
  83 {
  84         return x * pow(2, e);
  85 }
  86 float log10(float x)
  87 {
  88         return log(x) * M_LOG10E;
  89 }
  90 float log1p(float x)
  91 {
  92         return log(x + 1);
  93 }
  94 float log2(float x)
  95 {
  96         return log(x) * M_LOG2E;
  97 }
  98 float logb(float x)
  99 {
 100         return floor(log2(fabs(x)));
 101 }
 102 vector modf(float f)
 103 {
 104         return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
 105 }
 106
 107 float scalbn(float x, int n)
 108 {
 109         return x * pow(2, n);
 110 }
 111
 112 float cbrt(float x)
 113 {
 114         return copysign(pow(fabs(x), 1.0/3.0), x);
 115 }
 116 float hypot(float x, float y)
 117 {
 118         return sqrt(x*x + y*y);
 119 }
 120
 121 float erf(float x)
 122 {
 123         // approximation taken from wikipedia
 124         float y;
 125         y = x*x;
 126         return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
 127 }
 128 float erfc(float x)
 129 {
 130         return 1.0 - erf(x);
 131 }
 132 vector lgamma(float x)
 133 {
 134         // TODO improve accuracy
 135         if(!isfinite(x))
 136                 return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0';
 137         if(x < 1 && x == floor(x))
 138                 return nan("gamma") * '1 1 1';
 139         if(x < 0.1)
 140         {
 141                 vector v;
 142                 v = lgamma(1.0 - x);
 143                 // reflection formula:
 144                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
 145                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
 146                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
 147                 v_z = sin(M_PI * x);
 148                 v_x = log(M_PI) - log(fabs(v_z)) - v_x;
 149                 if(v_z < 0)
 150                         v_y = -v_y;
 151                 v_z = 0;
 152                 return v;
 153         }
 154         if(x < 1.1)
 155                 return lgamma(x + 1) - log(x) * '1 0 0';
 156         x -= 1;
 157         return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
 158 }
 159 float tgamma(float x)
 160 {
 161         vector v;
 162         v = lgamma(x);
 163         return exp(v_x) * v_y;
 164 }
 165
 166 float nearbyint(float x)
 167 {
 168         return rint(x);
 169 }
 170 float trunc(float x)
 171 {
 172         return (x>=0) ? floor(x) : ceil(x);
 173 }
 174
 175 float fmod(float x, float y)
 176 {
 177         return x - y * trunc(x / y);
 178 }
 179 float remainder(float x, float y)
 180 {
 181         return x - y * rint(x / y);
 182 }
 183 vector remquo(float x, float y)
 184 {
 185         vector v;
 186         v_z = 0;
 187         v_y = rint(x / y);
 188         v_x = x - y * v_y;
 189         return v;
 190 }
 191
 192 float copysign(float x, float y)
 193 {
 194         return fabs(x) * ((y>0) ? 1 : -1);
 195 }
 196 float nan(string tag)
 197 {
 198         return sqrt(-1);
 199 }
 200 float nextafter(float x, float y)
 201 {
 202         // TODO very crude
 203         if(x == y)
 204                 return nan("nextafter");
 205         if(x > y)
 206                 return -nextafter(-x, -y);
 207         // now we know that x < y
 208         // so we need the next number > x
 209         float d, a, b;
 210         d = max(fabs(x), 0.00000000000000000000001);
 211         a = x + d;
 212         do
 213         {
 214                 d *= 0.5;
 215                 b = a;
 216                 a = x + d;
 217         }
 218         while(a != x);
 219         return b;
 220 }
 221 float nexttoward(float x, float y)
 222 {
 223         return nextafter(x, y);
 224 }
 225
 226 float fdim(float x, float y)
 227 {
 228         return max(x-y, 0);
 229 }
 230 float fmax(float x, float y)
 231 {
 232         return max(x, y);
 233 }
 234 float fmin(float x, float y)
 235 {
 236         return min(x, y);
 237 }
 238 float fma(float x, float y, float z)
 239 {
 240         return x * y + z;
 241 }
 242
 243 int isgreater(float x, float y)
 244 {
 245         return x > y;
 246 }
 247 int isgreaterequal(float x, float y)
 248 {
 249         return x >= y;
 250 }
 251 int isless(float x, float y)
 252 {
 253         return x < y;
 254 }
 255 int islessequal(float x, float y)
 256 {
 257         return x <= y;
 258 }
 259 int islessgreater(float x, float y)
 260 {
 261         return x < y || x > y;
 262 }
 263 int isunordered(float x, float y)
 264 {
 265         return !(x < y || x == y || x > y);
 266 }