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