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 //float pow(float x, float y);
 169 //float sqrt(float x, float y);
 170
 171 float erf(float x)
 172 {
 173         // approximation taken from wikipedia
 174         float y;
 175         y = x*x;
 176         return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
 177 }
 178 float erfc(float x)
 179 {
 180         return 1.0 - erf(x);
 181 }
 182 vector lgamma(float x)
 183 {
 184         // TODO improve accuracy
 185         if(x < 1 && x == floor(x))
 186                 return nan("gamma") * '1 1 1';
 187         if(x < 0.1)
 188         {
 189                 vector v;
 190                 v = lgamma(1.0 - x);
 191                 // reflection formula:
 192                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
 193                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
 194                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
 195                 v_z = sin(M_PI * x);
 196                 v_x = log(M_PI) - log(fabs(v_z)) - v_x;
 197                 if(v_z < 0)
 198                         v_y = -v_y;
 199                 v_z = 0;
 200                 return v;
 201         }
 202         if(x < 1.1)
 203                 return lgamma(x + 1) - log(x) * '1 0 0';
 204         x -= 1;
 205         return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
 206 }
 207 float tgamma(float x)
 208 {
 209         vector v;
 210         v = lgamma(x);
 211         return exp(v_x) * v_y;
 212 }
 213
 214 //float ceil(float x);
 215 //float floor(float x);
 216 float nearbyint(float x)
 217 {
 218         return rint(x);
 219 }
 220 //float rint(float x);
 221 //float round(float x);
 222 float trunc(float x)
 223 {
 224         return (x>=0) ? floor(x) : ceil(x);
 225 }
 226
 227 float fmod(float x, float y)
 228 {
 229         return x - y * trunc(x / y);
 230 }
 231 float remainder(float x, float y)
 232 {
 233         return x - y * rint(x / y);
 234 }
 235 vector remquo(float x, float y)
 236 {
 237         vector v;
 238         v_z = 0;
 239         v_y = rint(x / y);
 240         v_x = x - y * v_y;
 241         return v;
 242 }
 243
 244 float copysign(float x, float y)
 245 {
 246         return fabs(x) * ((y>0) ? 1 : -1);
 247 }
 248 float nan(string tag)
 249 {
 250         return sqrt(-1);
 251 }
 252 float nextafter(float x, float y)
 253 {
 254         // TODO very crude
 255         if(x == y)
 256                 return nan("nextafter");
 257         if(x > y)
 258                 return -nextafter(-x, -y);
 259         // now we know that x < y
 260         // so we need the next number > x
 261         float d, a, b;
 262         d = max(fabs(x), 0.00000000000000000000001);
 263         a = x + d;
 264         do
 265         {
 266                 d *= 0.5;
 267                 b = a;
 268                 a = x + d;
 269         }
 270         while(a != x);
 271         return b;
 272 }
 273 float nexttoward(float x, float y)
 274 {
 275         return nextafter(x, y);
 276 }
 277
 278 float fdim(float x, float y)
 279 {
 280         return max(x-y, 0);
 281 }
 282 float fmax(float x, float y)
 283 {
 284         return max(x, y);
 285 }
 286 float fmin(float x, float y)
 287 {
 288         return min(x, y);
 289 }
 290 float fma(float x, float y, float z)
 291 {
 292         return x * y + z;
 293 }
 294
 295 int isgreater(float x, float y)
 296 {
 297         return x > y;
 298 }
 299 int isgreaterequal(float x, float y)
 300 {
 301         return x >= y;
 302 }
 303 int isless(float x, float y)
 304 {
 305         return x < y;
 306 }
 307 int islessequal(float x, float y)
 308 {
 309         return x <= y;
 310 }
 311 int islessgreater(float x, float y)
 312 {
 313         return x < y || x > y;
 314 }
 315 int isunordered(float x, float y)
 316 {
 317         return !(x < y || x == y || x > y);
 318 }