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