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