qcsrc/lib/warpzone/mathlib.qc

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