data/qcsrc/common/mathlib.qc

   1 /*
   2 Copyright (c) 2009 Rudolf Polzer
   3
   4 Permission is hereby granted, free of charge, to any person obtaining a copy of
   5 this software and associated documentation files (the "Software"), to deal in
   6 the Software without restriction, including without limitation the rights to
   7 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
   8 of the Software, and to permit persons to whom the Software is furnished to do
   9 so, subject to the following conditions:
  10
  11 The above copyright notice and this permission notice shall be included in all
  12 copies or substantial portions of the Software.
  13
  14 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  15 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  16 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  17 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  18 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  19 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  20 SOFTWARE.
  21 */
  22
  23 int fpclassify(float x)
  24 {
  25         if(isnan(x))
  26                 return FP_NAN;
  27         if(isinf(x))
  28                 return FP_INFINITE;
  29         if(x == 0)
  30                 return FP_ZERO;
  31         return FP_NORMAL;
  32 }
  33 int isfinite(float x)
  34 {
  35         return !(isnan(x) || isinf(x));
  36 }
  37 int isinf(float x)
  38 {
  39         return (x != 0) && (x + x == x);
  40 }
  41 int isnan(float x)
  42 {
  43         return isunordered(x, x);
  44 }
  45 int isnormal(float x)
  46 {
  47         return isfinite(x);
  48 }
  49 int signbit(float x)
  50 {
  51         return (x < 0);
  52 }
  53
  54 float acosh(float x)
  55 {
  56         return log(x + sqrt(x*x - 1));
  57 }
  58 float asinh(float x)
  59 {
  60         return log(x + sqrt(x*x + 1));
  61 }
  62 float atanh(float x)
  63 {
  64         return 0.5 * log((1+x) / (1-x));
  65 }
  66 float cosh(float x)
  67 {
  68         return 0.5 * (exp(x) + exp(-x));
  69 }
  70 float sinh(float x)
  71 {
  72         return 0.5 * (exp(x) - exp(-x));
  73 }
  74 float tanh(float x)
  75 {
  76         return sinh(x) / cosh(x);
  77 }
  78
  79 float exp(float x)
  80 {
  81         return pow(M_E, x);
  82 }
  83 float exp2(float x)
  84 {
  85         return pow(2, x);
  86 }
  87 float expm1(float x)
  88 {
  89         return exp(x) - 1;
  90 }
  91
  92 vector frexp(float x)
  93 {
  94         vector v;
  95         v_z = 0;
  96         v_y = ilogb(x);
  97         v_x = x / exp2(v_y);
  98         return v;
  99 }
 100 int ilogb(float x)
 101 {
 102         return floor(log2(x));
 103 }
 104 float ldexp(float x, int e)
 105 {
 106         return x * pow(2, e);
 107 }
 108 float log(float x)
 109 {
 110         // TODO improve speed
 111         float i;
 112         float r, r0;
 113         if(x <= 0)
 114                 return nan("log");
 115         if(!isfinite(x))
 116                 return x;
 117         if(x >= 8)
 118                 return -log(1 / x); // faster
 119         if(x < 0.0009765625)
 120                 return 2 * log(sqrt(x)); // faster
 121         r = 1;
 122         r0 = 0;
 123         for(i = 1; fabs(r - r0) >= 0.0000001; ++i)
 124         {
 125                 // Newton iteration on exp(r) = x:
 126                 //   r <- r - (exp(r) - x) / (exp(r))
 127                 //   r <- r - 1 + x / exp(r)
 128                 r0 = r;
 129                 r = r0 - 1 + x / exp(r0);
 130         }
 131         print(ftos(i), "\n");
 132         return r;
 133 }
 134 float log10(float x)
 135 {
 136         return log(x) * M_LOG10E;
 137 }
 138 float log1p(float x)
 139 {
 140         return log(x + 1);
 141 }
 142 float log2(float x)
 143 {
 144         return log(x) * M_LOG2E;
 145 }
 146 float logb(float x)
 147 {
 148         return floor(log2(x));
 149 }
 150 vector modf(float f)
 151 {
 152         return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
 153 }
 154
 155 float scalbn(float x, int n)
 156 {
 157         return x * pow(2, n);
 158 }
 159
 160 float cbrt(float x)
 161 {
 162         return pow(x, 1.0/3.0);
 163 }
 164 float hypot(float x, float y)
 165 {
 166         return sqrt(x*x + y*y);
 167 }
 168
 169 float erf(float x)
 170 {
 171         // approximation taken from wikipedia
 172         float y;
 173         y = x*x;
 174         return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
 175 }
 176 float erfc(float x)
 177 {
 178         return 1.0 - erf(x);
 179 }
 180 vector lgamma(float x)
 181 {
 182         // TODO improve accuracy
 183         if(x < 1 && x == floor(x))
 184                 return nan("gamma") * '1 1 1';
 185         if(x < 0.1)
 186         {
 187                 vector v;
 188                 v = lgamma(1.0 - x);
 189                 // reflection formula:
 190                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
 191                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
 192                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
 193                 v_z = sin(M_PI * x);
 194                 v_x = log(M_PI) - log(fabs(v_z)) - v_x;
 195                 if(v_z < 0)
 196                         v_y = -v_y;
 197                 v_z = 0;
 198                 return v;
 199         }
 200         if(x < 1.1)
 201                 return lgamma(x + 1) - log(x) * '1 0 0';
 202         x -= 1;
 203         return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
 204 }
 205 float tgamma(float x)
 206 {
 207         vector v;
 208         v = lgamma(x);
 209         return exp(v_x) * v_y;
 210 }
 211
 212 float nearbyint(float x)
 213 {
 214         return rint(x);
 215 }
 216 float trunc(float x)
 217 {
 218         return (x>=0) ? floor(x) : ceil(x);
 219 }
 220
 221 float fmod(float x, float y)
 222 {
 223         return x - y * trunc(x / y);
 224 }
 225 float remainder(float x, float y)
 226 {
 227         return x - y * rint(x / y);
 228 }
 229 vector remquo(float x, float y)
 230 {
 231         vector v;
 232         v_z = 0;
 233         v_y = rint(x / y);
 234         v_x = x - y * v_y;
 235         return v;
 236 }
 237
 238 float copysign(float x, float y)
 239 {
 240         return fabs(x) * ((y>0) ? 1 : -1);
 241 }
 242 float nan(string tag)
 243 {
 244         return sqrt(-1);
 245 }
 246 float nextafter(float x, float y)
 247 {
 248         // TODO very crude
 249         if(x == y)
 250                 return nan("nextafter");
 251         if(x > y)
 252                 return -nextafter(-x, -y);
 253         // now we know that x < y
 254         // so we need the next number > x
 255         float d, a, b;
 256         d = max(fabs(x), 0.00000000000000000000001);
 257         a = x + d;
 258         do
 259         {
 260                 d *= 0.5;
 261                 b = a;
 262                 a = x + d;
 263         }
 264         while(a != x);
 265         return b;
 266 }
 267 float nexttoward(float x, float y)
 268 {
 269         return nextafter(x, y);
 270 }
 271
 272 float fdim(float x, float y)
 273 {
 274         return max(x-y, 0);
 275 }
 276 float fmax(float x, float y)
 277 {
 278         return max(x, y);
 279 }
 280 float fmin(float x, float y)
 281 {
 282         return min(x, y);
 283 }
 284 float fma(float x, float y, float z)
 285 {
 286         return x * y + z;
 287 }
 288
 289 int isgreater(float x, float y)
 290 {
 291         return x > y;
 292 }
 293 int isgreaterequal(float x, float y)
 294 {
 295         return x >= y;
 296 }
 297 int isless(float x, float y)
 298 {
 299         return x < y;
 300 }
 301 int islessequal(float x, float y)
 302 {
 303         return x <= y;
 304 }
 305 int islessgreater(float x, float y)
 306 {
 307         return x < y || x > y;
 308 }
 309 int isunordered(float x, float y)
 310 {
 311         return !(x < y || x == y || x > y);
 312 }