qcsrc/lib/math.qh

   1 #pragma once
   2
   3 #include "lib/float.qh"
   4
   5 ERASEABLE
   6 void mean_accumulate(entity e, .float a, .float c, float mean, float value, float weight)
   7 {
   8         if (weight == 0) return;
   9         if (mean == 0) e.(a) *= (value ** weight);
  10         else e.(a) += (value ** mean) * weight;
  11         e.(c) += weight;
  12 }
  13
  14 ERASEABLE
  15 float mean_evaluate(entity e, .float a, .float c, float mean)
  16 {
  17         if (e.(c) == 0) return 0;
  18         if (mean == 0) return (e.(a) ** (1.0 / e.(c)));
  19         else return ((e.(a) / e.(c)) ** (1.0 / mean));
  20 }
  21
  22 #define MEAN_ACCUMULATE(s, prefix, v, w) mean_accumulate(s, prefix##_accumulator, prefix##_count, prefix##_mean, v, w)
  23 #define MEAN_EVALUATE(s, prefix) mean_evaluate(s, prefix##_accumulator, prefix##_count, prefix##_mean)
  24 #define MEAN_DECLARE(prefix, m) float prefix##_mean = m; .float prefix##_count, prefix##_accumulator
  25
  26 /** Returns a random number between -1.0 and 1.0 */
  27 #define crandom() (2 * (random() - 0.5))
  28
  29
  30 /*
  31 ==================
  32 Angc used for animations
  33 ==================
  34 */
  35
  36
  37 ERASEABLE
  38 float angc(float a1, float a2)
  39 {
  40         while (a1 > 180)
  41                 a1 -= 360;
  42         while (a1 < -179)
  43                 a1 += 360;
  44         while (a2 > 180)
  45                 a2 -= 360;
  46         while (a2 < -179)
  47                 a2 += 360;
  48         float a = a1 - a2;
  49         while (a > 180)
  50                 a -= 360;
  51         while (a < -179)
  52                 a += 360;
  53         return a;
  54 }
  55
  56 ERASEABLE
  57 float fsnap(float val, float fsize)
  58 {
  59         return rint(val / fsize) * fsize;
  60 }
  61
  62 ERASEABLE
  63 vector vsnap(vector point, float fsize)
  64 {
  65         vector vret;
  66
  67         vret.x = rint(point.x / fsize) * fsize;
  68         vret.y = rint(point.y / fsize) * fsize;
  69         vret.z = ceil(point.z / fsize) * fsize;
  70
  71         return vret;
  72 }
  73
  74 ERASEABLE
  75 float lerpratio(float f0, float f1, float ratio)
  76 {
  77         return f0 * (1 - ratio) + f1 * ratio;
  78 }
  79
  80 ERASEABLE
  81 float lerp(float t0, float f0, float t1, float f1, float t)
  82 {
  83         return lerpratio(f0, f1, (t - t0) / (t1 - t0));
  84 }
  85
  86 ERASEABLE
  87 float lerp3ratio(float f0, float f1, float f2, float ratio)
  88 {
  89         float mid = 0.5;
  90         return ratio < mid ? lerpratio(f0, f1, ratio / mid) : ratio > mid ? lerpratio(f1, f2, (ratio - mid) / mid) : f1;
  91 }
  92
  93
  94 ERASEABLE
  95 vector lerpvratio(vector f0, vector f1, float ratio)
  96 {
  97         return f0 * (1 - ratio) + f1 * ratio;
  98 }
  99
 100 ERASEABLE
 101 vector lerpv3ratio(vector f0, vector f1, vector f2, float ratio)
 102 {
 103         float mid = 0.5;
 104         return ratio < mid ? lerpvratio(f0, f1, ratio / mid) : ratio > mid ? lerpvratio(f1, f2, (ratio - mid) / mid) : f1;
 105 }
 106
 107 ERASEABLE
 108 vector lerpv(float t0, vector v0, float t1, vector v1, float t)
 109 {
 110         return v0 + (v1 - v0) * ((t - t0) / (t1 - t0));
 111 }
 112
 113 ERASEABLE
 114 vector bezier_quadratic_getpoint(vector a, vector b, vector c, float t)
 115 {
 116         return (c - 2 * b + a) * (t * t)
 117                + (b - a) * (2 * t)
 118                + a;
 119 }
 120
 121 ERASEABLE
 122 vector bezier_quadratic_getderivative(vector a, vector b, vector c, float t)
 123 {
 124         return (c - 2 * b + a) * (2 * t)
 125                + (b - a) * 2;
 126 }
 127
 128 ERASEABLE
 129 float cubic_speedfunc(float startspeedfactor, float endspeedfactor, float spd)
 130 {
 131         return (((startspeedfactor + endspeedfactor - 2
 132                  ) * spd - 2 * startspeedfactor - endspeedfactor + 3
 133                 ) * spd + startspeedfactor
 134                ) * spd;
 135 }
 136
 137 ERASEABLE
 138 bool cubic_speedfunc_is_sane(float startspeedfactor, float endspeedfactor)
 139 {
 140         if (startspeedfactor < 0 || endspeedfactor < 0) return false;
 141
 142         /*
 143         // if this is the case, the possible zeros of the first derivative are outside
 144         // 0..1
 145         We can calculate this condition as condition
 146         if(se <= 3)
 147             return true;
 148         */
 149
 150         // better, see below:
 151         if (startspeedfactor <= 3 && endspeedfactor <= 3) return true;
 152
 153         // if this is the case, the first derivative has no zeros at all
 154         float se = startspeedfactor + endspeedfactor;
 155         float s_e = startspeedfactor - endspeedfactor;
 156         if (3 * (se - 4) * (se - 4) + s_e * s_e <= 12)  // an ellipse
 157                 return true;
 158
 159         // Now let s <= 3, s <= 3, s+e >= 3 (triangle) then we get se <= 6 (top right corner).
 160         // we also get s_e <= 6 - se
 161         // 3 * (se - 4)^2 + (6 - se)^2
 162         // is quadratic, has value 12 at 3 and 6, and value < 12 in between.
 163         // Therefore, above "better" check works!
 164
 165         return false;
 166
 167         // known good cases:
 168         // (0, [0..3])
 169         // (0.5, [0..3.8])
 170         // (1, [0..4])
 171         // (1.5, [0..3.9])
 172         // (2, [0..3.7])
 173         // (2.5, [0..3.4])
 174         // (3, [0..3])
 175         // (3.5, [0.2..2.3])
 176         // (4, 1)
 177
 178         /*
 179            On another note:
 180            inflection point is always at (2s + e - 3) / (3s + 3e - 6).
 181
 182            s + e - 2 == 0: no inflection
 183
 184            s + e > 2:
 185            0 < inflection < 1 if:
 186            0 < 2s + e - 3 < 3s + 3e - 6
 187            2s + e > 3 and 2e + s > 3
 188
 189            s + e < 2:
 190            0 < inflection < 1 if:
 191            0 > 2s + e - 3 > 3s + 3e - 6
 192            2s + e < 3 and 2e + s < 3
 193
 194            Therefore: there is an inflection point iff:
 195            e outside (3 - s)/2 .. 3 - s*2
 196
 197            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)
 198         */
 199 }
 200
 201 /** continuous function mapping all reals into -1..1 */
 202 ERASEABLE
 203 float float2range11(float f)
 204 {
 205         return f / (fabs(f) + 1);
 206 }
 207
 208 /** continuous function mapping all reals into 0..1 */
 209 ERASEABLE
 210 float float2range01(float f)
 211 {
 212         return 0.5 + 0.5 * float2range11(f);
 213 }
 214
 215 ERASEABLE
 216 float median(float a, float b, float c)
 217 {
 218         return (a < c) ? bound(a, b, c) : bound(c, b, a);
 219 }
 220
 221 ERASEABLE
 222 float almost_equals(float a, float b)
 223 {
 224         float eps = (max(a, -a) + max(b, -b)) * 0.001;
 225         return a - b < eps && b - a < eps;
 226 }
 227
 228 ERASEABLE
 229 float almost_equals_eps(float a, float b, float times_eps)
 230 {
 231         float eps = max(fabs(a), fabs(b)) * FLOAT_EPSILON * times_eps;
 232         return a - b < eps && b - a < eps;
 233 }
 234
 235 ERASEABLE
 236 float almost_in_bounds(float a, float b, float c)
 237 {
 238         float eps = (max(a, -a) + max(c, -c)) * 0.001;
 239         if (a > c) eps = -eps;
 240         return b == median(a - eps, b, c + eps);
 241 }
 242
 243 ERASEABLE
 244 float ExponentialFalloff(float mindist, float maxdist, float halflifedist, float d)
 245 {
 246         if (halflifedist > 0) return (0.5 ** ((bound(mindist, d, maxdist) - mindist) / halflifedist));
 247         else if (halflifedist < 0) return (0.5 ** ((bound(mindist, d, maxdist) - maxdist) / halflifedist));
 248         else return 1;
 249 }
 250
 251 #define power2of(e) (2 ** e)
 252
 253 ERASEABLE
 254 float log2of(float e)
 255 {
 256         // NOTE: generated code
 257         if (e > 2048)
 258                 if (e > 131072)
 259                         if (e > 1048576)
 260                                 if (e > 4194304) return 23;
 261                                 else
 262                                         if (e > 2097152) return 22;
 263                                         else return 21;
 264                         else
 265                                 if (e > 524288) return 20;
 266                                 else
 267                                         if (e > 262144) return 19;
 268                                         else return 18;
 269                 else
 270                         if (e > 16384)
 271                                 if (e > 65536) return 17;
 272                                 else
 273                                         if (e > 32768) return 16;
 274                                         else return 15;
 275                         else
 276                                 if (e > 8192) return 14;
 277                                 else
 278                                         if (e > 4096) return 13;
 279                                         else return 12;
 280         else
 281                 if (e > 32)
 282                         if (e > 256)
 283                                 if (e > 1024) return 11;
 284                                 else
 285                                         if (e > 512) return 10;
 286                                         else return 9;
 287                         else
 288                                 if (e > 128) return 8;
 289                                 else
 290                                         if (e > 64) return 7;
 291                                         else return 6;
 292                 else
 293                         if (e > 4)
 294                                 if (e > 16) return 5;
 295                                 else
 296                                         if (e > 8) return 4;
 297                                         else return 3;
 298                         else
 299                                 if (e > 2) return 2;
 300                                 else
 301                                         if (e > 1) return 1;
 302                                         else return 0;
 303 }
 304
 305 /** ax^2 + bx + c = 0 */
 306 ERASEABLE
 307 vector solve_quadratic(float a, float b, float c)
 308 {
 309         vector v;
 310         float D;
 311         v = '0 0 0';
 312         if (a == 0)
 313         {
 314                 if (b != 0)
 315                 {
 316                         v.x = v.y = -c / b;
 317                         v.z = 1;
 318                 }
 319                 else
 320                 {
 321                         if (c == 0)
 322                         {
 323                                 // actually, every number solves the equation!
 324                                 v.z = 1;
 325                         }
 326                 }
 327         }
 328         else
 329         {
 330                 D = b * b - 4 * a * c;
 331                 if (D >= 0)
 332                 {
 333                         D = sqrt(D);
 334                         if (a > 0)  // put the smaller solution first
 335                         {
 336                                 v.x = ((-b) - D) / (2 * a);
 337                                 v.y = ((-b) + D) / (2 * a);
 338                         }
 339                         else
 340                         {
 341                                 v.x = (-b + D) / (2 * a);
 342                                 v.y = (-b - D) / (2 * a);
 343                         }
 344                         v.z = 1;
 345                 }
 346                 else
 347                 {
 348                         // complex solutions!
 349                         D = sqrt(-D);
 350                         v.x = -b / (2 * a);
 351                         if (a > 0) v.y =  D / (2 * a);
 352                         else v.y = -D / (2 * a);
 353                         v.z = 0;
 354                 }
 355         }
 356         return v;
 357 }
 358
 359 /// Maps values between the src and dest range: src_min to dest_min, src_max to dest_max, values between them
 360 /// to the corresponding values between and extrapolates for values outside the range.
 361 ///
 362 /// src_min and src_max must not be the same or division by zero occurs.
 363 ///
 364 /// dest_max can be smaller than dest_min if you want the resulting range to be inverted, all values can be negative.
 365 ERASEABLE
 366 float map_ranges(float value, float src_min, float src_max, float dest_min, float dest_max) {
 367         float src_diff = src_max - src_min;
 368         float dest_diff = dest_max - dest_min;
 369         float ratio = (value - src_min) / src_diff;
 370         return dest_min + dest_diff * ratio;
 371 }
 372
 373 /// Same as `map_ranges` except that values outside the source range are clamped to min or max.
 374 ERASEABLE
 375 float map_bound_ranges(float value, float src_min, float src_max, float dest_min, float dest_max) {
 376         if (value <= src_min) return dest_min;
 377         if (value >= src_max) return dest_max;
 378         return map_ranges(value, src_min, src_max, dest_min, dest_max);
 379 }