qcsrc/warpzonelib/mathlib.qc

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