qcsrc/warpzonelib/mathlib.qc

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