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