qcsrc/lib/math.qh

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