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