qcsrc/lib/warpzone/mathlib.qc

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