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