qcsrc/lib/math.qh

   1 #pragma once
   2
   3 void mean_accumulate(entity e, .float a, .float c, float mean, float value, float weight)
   4 {
   5         if (weight == 0) return;
   6         if (mean == 0) e.(a) *= pow(value, weight);
   7         else e.(a) += pow(value, mean) * weight;
   8         e.(c) += weight;
   9 }
  10
  11 float mean_evaluate(entity e, .float a, .float c, float mean)
  12 {
  13         if (e.(c) == 0) return 0;
  14         if (mean == 0) return pow(e.(a), 1.0 / e.(c));
  15         else return pow(e.(a) / e.(c), 1.0 / mean);
  16 }
  17
  18 #define MEAN_ACCUMULATE(s, prefix, v, w) mean_accumulate(s, prefix##_accumulator, prefix##_count, prefix##_mean, v, w)
  19 #define MEAN_EVALUATE(s, prefix) mean_evaluate(s, prefix##_accumulator, prefix##_count, prefix##_mean)
  20 #define MEAN_DECLARE(prefix, m) float prefix##_mean = m; .float prefix##_count, prefix##_accumulator
  21
  22 /** Returns a random number between -1.0 and 1.0 */
  23 #define crandom() (2 * (random() - 0.5))
  24
  25
  26 /*
  27 ==================
  28 Angc used for animations
  29 ==================
  30 */
  31
  32
  33 float angc(float a1, float a2)
  34 {
  35         while (a1 > 180)
  36                 a1 -= 360;
  37         while (a1 < -179)
  38                 a1 += 360;
  39         while (a2 > 180)
  40                 a2 -= 360;
  41         while (a2 < -179)
  42                 a2 += 360;
  43         float a = a1 - a2;
  44         while (a > 180)
  45                 a -= 360;
  46         while (a < -179)
  47                 a += 360;
  48         return a;
  49 }
  50
  51 float fsnap(float val, float fsize)
  52 {
  53         return rint(val / fsize) * fsize;
  54 }
  55
  56 vector vsnap(vector point, float fsize)
  57 {
  58         vector vret;
  59
  60         vret.x = rint(point.x / fsize) * fsize;
  61         vret.y = rint(point.y / fsize) * fsize;
  62         vret.z = ceil(point.z / fsize) * fsize;
  63
  64         return vret;
  65 }
  66
  67 vector lerpv(float t0, vector v0, float t1, vector v1, float t)
  68 {
  69         return v0 + (v1 - v0) * ((t - t0) / (t1 - t0));
  70 }
  71
  72 vector bezier_quadratic_getpoint(vector a, vector b, vector c, float t)
  73 {
  74         return (c - 2 * b + a) * (t * t)
  75                + (b - a) * (2 * t)
  76                + a;
  77 }
  78
  79 vector bezier_quadratic_getderivative(vector a, vector b, vector c, float t)
  80 {
  81         return (c - 2 * b + a) * (2 * t)
  82                + (b - a) * 2;
  83 }
  84
  85 float cubic_speedfunc(float startspeedfactor, float endspeedfactor, float x)
  86 {
  87         return (((startspeedfactor + endspeedfactor - 2
  88                  ) * x - 2 * startspeedfactor - endspeedfactor + 3
  89                 ) * x + startspeedfactor
  90                ) * x;
  91 }
  92
  93 bool cubic_speedfunc_is_sane(float startspeedfactor, float endspeedfactor)
  94 {
  95         if (startspeedfactor < 0 || endspeedfactor < 0) return false;
  96
  97         /*
  98         // if this is the case, the possible zeros of the first derivative are outside
  99         // 0..1
 100         We can calculate this condition as condition
 101         if(se <= 3)
 102             return true;
 103         */
 104
 105         // better, see below:
 106         if (startspeedfactor <= 3 && endspeedfactor <= 3) return true;
 107
 108         // if this is the case, the first derivative has no zeros at all
 109         float se = startspeedfactor + endspeedfactor;
 110         float s_e = startspeedfactor - endspeedfactor;
 111         if (3 * (se - 4) * (se - 4) + s_e * s_e <= 12)  // an ellipse
 112                 return true;
 113
 114         // Now let s <= 3, s <= 3, s+e >= 3 (triangle) then we get se <= 6 (top right corner).
 115         // we also get s_e <= 6 - se
 116         // 3 * (se - 4)^2 + (6 - se)^2
 117         // is quadratic, has value 12 at 3 and 6, and value < 12 in between.
 118         // Therefore, above "better" check works!
 119
 120         return false;
 121
 122         // known good cases:
 123         // (0, [0..3])
 124         // (0.5, [0..3.8])
 125         // (1, [0..4])
 126         // (1.5, [0..3.9])
 127         // (2, [0..3.7])
 128         // (2.5, [0..3.4])
 129         // (3, [0..3])
 130         // (3.5, [0.2..2.3])
 131         // (4, 1)
 132
 133         /*
 134            On another note:
 135            inflection point is always at (2s + e - 3) / (3s + 3e - 6).
 136
 137            s + e - 2 == 0: no inflection
 138
 139            s + e > 2:
 140            0 < inflection < 1 if:
 141            0 < 2s + e - 3 < 3s + 3e - 6
 142            2s + e > 3 and 2e + s > 3
 143
 144            s + e < 2:
 145            0 < inflection < 1 if:
 146            0 > 2s + e - 3 > 3s + 3e - 6
 147            2s + e < 3 and 2e + s < 3
 148
 149            Therefore: there is an inflection point iff:
 150            e outside (3 - s)/2 .. 3 - s*2
 151
 152            in other words, if (s,e) in triangle (1,1)(0,3)(0,1.5) or in triangle (1,1)(3,0)(1.5,0)
 153         */
 154 }
 155
 156 /** continuous function mapping all reals into -1..1 */
 157 float float2range11(float f)
 158 {
 159         return f / (fabs(f) + 1);
 160 }
 161
 162 /** continuous function mapping all reals into 0..1 */
 163 float float2range01(float f)
 164 {
 165         return 0.5 + 0.5 * float2range11(f);
 166 }
 167
 168 float median(float a, float b, float c)
 169 {
 170         return (a < c) ? bound(a, b, c) : bound(c, b, a);
 171 }
 172
 173 float almost_equals(float a, float b)
 174 {
 175         float eps = (max(a, -a) + max(b, -b)) * 0.001;
 176         return a - b < eps && b - a < eps;
 177 }
 178
 179 float almost_in_bounds(float a, float b, float c)
 180 {
 181         float eps = (max(a, -a) + max(c, -c)) * 0.001;
 182         if (a > c) eps = -eps;
 183         return b == median(a - eps, b, c + eps);
 184 }
 185
 186 float ExponentialFalloff(float mindist, float maxdist, float halflifedist, float d)
 187 {
 188         if (halflifedist > 0) return pow(0.5, (bound(mindist, d, maxdist) - mindist) / halflifedist);
 189         else if (halflifedist < 0) return pow(0.5, (bound(mindist, d, maxdist) - maxdist) / halflifedist);
 190         else return 1;
 191 }
 192
 193 float power2of(float e)
 194 {
 195         return pow(2, e);
 196 }
 197
 198 float log2of(float x)
 199 {
 200         // NOTE: generated code
 201         if (x > 2048)
 202                 if (x > 131072)
 203                         if (x > 1048576)
 204                                 if (x > 4194304) return 23;
 205                                 else
 206                                         if (x > 2097152) return 22;
 207                                         else return 21;
 208                         else
 209                                 if (x > 524288) return 20;
 210                                 else
 211                                         if (x > 262144) return 19;
 212                                         else return 18;
 213                 else
 214                         if (x > 16384)
 215                                 if (x > 65536) return 17;
 216                                 else
 217                                         if (x > 32768) return 16;
 218                                         else return 15;
 219                         else
 220                                 if (x > 8192) return 14;
 221                                 else
 222                                         if (x > 4096) return 13;
 223                                         else return 12;
 224         else
 225                 if (x > 32)
 226                         if (x > 256)
 227                                 if (x > 1024) return 11;
 228                                 else
 229                                         if (x > 512) return 10;
 230                                         else return 9;
 231                         else
 232                                 if (x > 128) return 8;
 233                                 else
 234                                         if (x > 64) return 7;
 235                                         else return 6;
 236                 else
 237                         if (x > 4)
 238                                 if (x > 16) return 5;
 239                                 else
 240                                         if (x > 8) return 4;
 241                                         else return 3;
 242                         else
 243                                 if (x > 2) return 2;
 244                                 else
 245                                         if (x > 1) return 1;
 246                                         else return 0;
 247 }
 248
 249 /** ax^2 + bx + c = 0 */
 250 vector solve_quadratic(float a, float b, float c)
 251 {
 252         vector v;
 253         float D;
 254         v = '0 0 0';
 255         if (a == 0)
 256         {
 257                 if (b != 0)
 258                 {
 259                         v.x = v.y = -c / b;
 260                         v.z = 1;
 261                 }
 262                 else
 263                 {
 264                         if (c == 0)
 265                         {
 266                                 // actually, every number solves the equation!
 267                                 v.z = 1;
 268                         }
 269                 }
 270         }
 271         else
 272         {
 273                 D = b * b - 4 * a * c;
 274                 if (D >= 0)
 275                 {
 276                         D = sqrt(D);
 277                         if (a > 0)  // put the smaller solution first
 278                         {
 279                                 v.x = ((-b) - D) / (2 * a);
 280                                 v.y = ((-b) + D) / (2 * a);
 281                         }
 282                         else
 283                         {
 284                                 v.x = (-b + D) / (2 * a);
 285                                 v.y = (-b - D) / (2 * a);
 286                         }
 287                         v.z = 1;
 288                 }
 289                 else
 290                 {
 291                         // complex solutions!
 292                         D = sqrt(-D);
 293                         v.x = -b / (2 * a);
 294                         if (a > 0) v.y =  D / (2 * a);
 295                         else v.y = -D / (2 * a);
 296                         v.z = 0;
 297                 }
 298         }
 299         return v;
 300 }