data/qcsrc/common/mathlib.qc

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