math.c (20102B)
1 #include "external/cephes.c" 2 3 function void 4 fill_kronecker_sub_matrix(f32 *out, i32 out_stride, f32 scale, f32 *b, iv2 b_dim) 5 { 6 f32x4 vscale = dup_f32x4((f32)scale); 7 for (i32 i = 0; i < b_dim.y; i++) { 8 for (i32 j = 0; j < b_dim.x; j += 4, b += 4) { 9 f32x4 vb = load_f32x4(b); 10 store_f32x4(out + j, mul_f32x4(vscale, vb)); 11 } 12 out += out_stride; 13 } 14 } 15 16 /* NOTE: this won't check for valid space/etc and assumes row major order */ 17 function void 18 kronecker_product(f32 *out, f32 *a, iv2 a_dim, f32 *b, iv2 b_dim) 19 { 20 iv2 out_dim = {{a_dim.x * b_dim.x, a_dim.y * b_dim.y}}; 21 assert(out_dim.y % 4 == 0); 22 for (i32 i = 0; i < a_dim.y; i++) { 23 f32 *vout = out; 24 for (i32 j = 0; j < a_dim.x; j++, a++) { 25 fill_kronecker_sub_matrix(vout, out_dim.y, *a, b, b_dim); 26 vout += b_dim.y; 27 } 28 out += out_dim.y * b_dim.x; 29 } 30 } 31 32 /* NOTE/TODO: to support even more hadamard sizes use the Paley construction */ 33 function f32 * 34 make_hadamard_transpose(Arena *a, i32 dim) 35 { 36 read_only local_persist f32 hadamard_12_12_transpose[] = { 37 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 38 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 39 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 40 1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, 41 1, 1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 42 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 43 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, -1, 44 1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, 45 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 46 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 47 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 48 1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 49 }; 50 51 read_only local_persist f32 hadamard_20_20_transpose[] = { 52 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 53 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 54 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, 55 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 56 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 57 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 58 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 59 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 60 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 61 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 62 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 63 1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 64 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 65 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 66 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 67 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 68 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 69 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, 70 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 71 1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 72 }; 73 74 f32 *result = 0; 75 76 b32 power_of_2 = ISPOWEROF2(dim); 77 b32 multiple_of_12 = dim % 12 == 0; 78 b32 multiple_of_20 = dim % 20 == 0; 79 iz elements = dim * dim; 80 81 i32 base_dim = 0; 82 if (power_of_2) { 83 base_dim = dim; 84 } else if (multiple_of_20 && ISPOWEROF2(dim / 20)) { 85 base_dim = 20; 86 dim /= 20; 87 } else if (multiple_of_12 && ISPOWEROF2(dim / 12)) { 88 base_dim = 12; 89 dim /= 12; 90 } 91 92 if (ISPOWEROF2(dim) && base_dim && arena_capacity(a, f32) >= elements * (1 + (dim != base_dim))) { 93 result = push_array(a, f32, elements); 94 95 Arena tmp = *a; 96 f32 *m = dim == base_dim ? result : push_array(&tmp, f32, elements); 97 98 #define IND(i, j) ((i) * dim + (j)) 99 m[0] = 1; 100 for (i32 k = 1; k < dim; k *= 2) { 101 for (i32 i = 0; i < k; i++) { 102 for (i32 j = 0; j < k; j++) { 103 f32 val = m[IND(i, j)]; 104 m[IND(i + k, j)] = val; 105 m[IND(i, j + k)] = val; 106 m[IND(i + k, j + k)] = -val; 107 } 108 } 109 } 110 #undef IND 111 112 f32 *m2 = 0; 113 iv2 m2_dim; 114 switch (base_dim) { 115 case 12:{ m2 = hadamard_12_12_transpose; m2_dim = (iv2){{12, 12}}; }break; 116 case 20:{ m2 = hadamard_20_20_transpose; m2_dim = (iv2){{20, 20}}; }break; 117 } 118 if (m2) kronecker_product(result, m, (iv2){{dim, dim}}, m2, m2_dim); 119 } 120 121 return result; 122 } 123 124 function b32 125 u128_equal(u128 a, u128 b) 126 { 127 b32 result = a.U64[0] == b.U64[0] && a.U64[1] == b.U64[1]; 128 return result; 129 } 130 131 function b32 132 iv2_equal(iv2 a, iv2 b) 133 { 134 b32 result = a.x == b.x && a.y == b.y; 135 return result; 136 } 137 138 function b32 139 iv3_equal(iv3 a, iv3 b) 140 { 141 b32 result = a.x == b.x && a.y == b.y && a.z == b.z; 142 return result; 143 } 144 145 function v2 146 clamp_v2_rect(v2 v, Rect r) 147 { 148 v2 result = v; 149 result.x = CLAMP(v.x, r.pos.x, r.pos.x + r.size.x); 150 result.y = CLAMP(v.y, r.pos.y, r.pos.y + r.size.y); 151 return result; 152 } 153 154 function v2 155 v2_scale(v2 a, f32 scale) 156 { 157 v2 result; 158 result.x = a.x * scale; 159 result.y = a.y * scale; 160 return result; 161 } 162 163 function v2 164 v2_add(v2 a, v2 b) 165 { 166 v2 result; 167 result.x = a.x + b.x; 168 result.y = a.y + b.y; 169 return result; 170 } 171 172 function v2 173 v2_sub(v2 a, v2 b) 174 { 175 v2 result = v2_add(a, v2_scale(b, -1.0f)); 176 return result; 177 } 178 179 function v2 180 v2_mul(v2 a, v2 b) 181 { 182 v2 result; 183 result.x = a.x * b.x; 184 result.y = a.y * b.y; 185 return result; 186 } 187 188 function v2 189 v2_div(v2 a, v2 b) 190 { 191 v2 result; 192 result.x = a.x / b.x; 193 result.y = a.y / b.y; 194 return result; 195 } 196 197 function v2 198 v2_floor(v2 a) 199 { 200 v2 result; 201 result.x = (f32)((i32)a.x); 202 result.y = (f32)((i32)a.y); 203 return result; 204 } 205 206 function f32 207 v2_magnitude_squared(v2 a) 208 { 209 f32 result = a.x * a.x + a.y * a.y; 210 return result; 211 } 212 213 function f32 214 v2_magnitude(v2 a) 215 { 216 f32 result = sqrt_f32(a.x * a.x + a.y * a.y); 217 return result; 218 } 219 220 function v3 221 cross(v3 a, v3 b) 222 { 223 v3 result; 224 result.x = a.y * b.z - a.z * b.y; 225 result.y = a.z * b.x - a.x * b.z; 226 result.z = a.x * b.y - a.y * b.x; 227 return result; 228 } 229 230 function v3 231 v3_from_iv3(iv3 v) 232 { 233 v3 result; 234 result.E[0] = (f32)v.E[0]; 235 result.E[1] = (f32)v.E[1]; 236 result.E[2] = (f32)v.E[2]; 237 return result; 238 } 239 240 function v3 241 v3_from_f32_array(f32 v[3]) 242 { 243 v3 result; 244 result.E[0] = v[0]; 245 result.E[1] = v[1]; 246 result.E[2] = v[2]; 247 return result; 248 } 249 250 function v3 251 v3_abs(v3 a) 252 { 253 v3 result; 254 result.x = ABS(a.x); 255 result.y = ABS(a.y); 256 result.z = ABS(a.z); 257 return result; 258 } 259 260 function v3 261 v3_scale(v3 a, f32 scale) 262 { 263 v3 result; 264 result.x = scale * a.x; 265 result.y = scale * a.y; 266 result.z = scale * a.z; 267 return result; 268 } 269 270 function v3 271 v3_add(v3 a, v3 b) 272 { 273 v3 result; 274 result.x = a.x + b.x; 275 result.y = a.y + b.y; 276 result.z = a.z + b.z; 277 return result; 278 } 279 280 function v3 281 v3_sub(v3 a, v3 b) 282 { 283 v3 result = v3_add(a, v3_scale(b, -1.0f)); 284 return result; 285 } 286 287 function v3 288 v3_div(v3 a, v3 b) 289 { 290 v3 result; 291 result.x = a.x / b.x; 292 result.y = a.y / b.y; 293 result.z = a.z / b.z; 294 return result; 295 } 296 297 function f32 298 v3_dot(v3 a, v3 b) 299 { 300 f32 result = a.x * b.x + a.y * b.y + a.z * b.z; 301 return result; 302 } 303 304 function f32 305 v3_magnitude_squared(v3 a) 306 { 307 f32 result = v3_dot(a, a); 308 return result; 309 } 310 311 function f32 312 v3_magnitude(v3 a) 313 { 314 f32 result = sqrt_f32(v3_dot(a, a)); 315 return result; 316 } 317 318 function v3 319 v3_normalize(v3 a) 320 { 321 v3 result = v3_scale(a, 1.0f / v3_magnitude(a)); 322 return result; 323 } 324 325 function v4 326 v4_scale(v4 a, f32 scale) 327 { 328 v4 result; 329 result.x = scale * a.x; 330 result.y = scale * a.y; 331 result.z = scale * a.z; 332 result.w = scale * a.w; 333 return result; 334 } 335 336 function v4 337 v4_add(v4 a, v4 b) 338 { 339 v4 result; 340 result.x = a.x + b.x; 341 result.y = a.y + b.y; 342 result.z = a.z + b.z; 343 result.w = a.w + b.w; 344 return result; 345 } 346 347 function v4 348 v4_sub(v4 a, v4 b) 349 { 350 v4 result = v4_add(a, v4_scale(b, -1)); 351 return result; 352 } 353 354 function f32 355 v4_dot(v4 a, v4 b) 356 { 357 f32 result = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; 358 return result; 359 } 360 361 function v4 362 v4_lerp(v4 a, v4 b, f32 t) 363 { 364 v4 result = v4_add(a, v4_scale(v4_sub(b, a), t)); 365 return result; 366 } 367 368 function m4 369 m4_identity(void) 370 { 371 m4 result; 372 result.c[0] = (v4){{1, 0, 0, 0}}; 373 result.c[1] = (v4){{0, 1, 0, 0}}; 374 result.c[2] = (v4){{0, 0, 1, 0}}; 375 result.c[3] = (v4){{0, 0, 0, 1}}; 376 return result; 377 } 378 379 function v4 380 m4_row(m4 a, u32 row) 381 { 382 v4 result; 383 result.E[0] = a.c[0].E[row]; 384 result.E[1] = a.c[1].E[row]; 385 result.E[2] = a.c[2].E[row]; 386 result.E[3] = a.c[3].E[row]; 387 return result; 388 } 389 390 function m4 391 m4_mul(m4 a, m4 b) 392 { 393 m4 result; 394 for (u32 i = 0; i < 4; i++) { 395 for (u32 j = 0; j < 4; j++) { 396 result.c[i].E[j] = v4_dot(m4_row(a, j), b.c[i]); 397 } 398 } 399 return result; 400 } 401 402 /* NOTE(rnp): based on: 403 * https://web.archive.org/web/20131215123403/ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf 404 * TODO(rnp): redo with SIMD as given in the link (but need to rewrite for column-major) 405 */ 406 function m4 407 m4_inverse(m4 m) 408 { 409 m4 result; 410 result.E[ 0] = m.E[5] * m.E[10] * m.E[15] - m.E[5] * m.E[11] * m.E[14] - m.E[9] * m.E[6] * m.E[15] + m.E[9] * m.E[7] * m.E[14] + m.E[13] * m.E[6] * m.E[11] - m.E[13] * m.E[7] * m.E[10]; 411 result.E[ 4] = -m.E[4] * m.E[10] * m.E[15] + m.E[4] * m.E[11] * m.E[14] + m.E[8] * m.E[6] * m.E[15] - m.E[8] * m.E[7] * m.E[14] - m.E[12] * m.E[6] * m.E[11] + m.E[12] * m.E[7] * m.E[10]; 412 result.E[ 8] = m.E[4] * m.E[ 9] * m.E[15] - m.E[4] * m.E[11] * m.E[13] - m.E[8] * m.E[5] * m.E[15] + m.E[8] * m.E[7] * m.E[13] + m.E[12] * m.E[5] * m.E[11] - m.E[12] * m.E[7] * m.E[ 9]; 413 result.E[12] = -m.E[4] * m.E[ 9] * m.E[14] + m.E[4] * m.E[10] * m.E[13] + m.E[8] * m.E[5] * m.E[14] - m.E[8] * m.E[6] * m.E[13] - m.E[12] * m.E[5] * m.E[10] + m.E[12] * m.E[6] * m.E[ 9]; 414 result.E[ 1] = -m.E[1] * m.E[10] * m.E[15] + m.E[1] * m.E[11] * m.E[14] + m.E[9] * m.E[2] * m.E[15] - m.E[9] * m.E[3] * m.E[14] - m.E[13] * m.E[2] * m.E[11] + m.E[13] * m.E[3] * m.E[10]; 415 result.E[ 5] = m.E[0] * m.E[10] * m.E[15] - m.E[0] * m.E[11] * m.E[14] - m.E[8] * m.E[2] * m.E[15] + m.E[8] * m.E[3] * m.E[14] + m.E[12] * m.E[2] * m.E[11] - m.E[12] * m.E[3] * m.E[10]; 416 result.E[ 9] = -m.E[0] * m.E[ 9] * m.E[15] + m.E[0] * m.E[11] * m.E[13] + m.E[8] * m.E[1] * m.E[15] - m.E[8] * m.E[3] * m.E[13] - m.E[12] * m.E[1] * m.E[11] + m.E[12] * m.E[3] * m.E[ 9]; 417 result.E[13] = m.E[0] * m.E[ 9] * m.E[14] - m.E[0] * m.E[10] * m.E[13] - m.E[8] * m.E[1] * m.E[14] + m.E[8] * m.E[2] * m.E[13] + m.E[12] * m.E[1] * m.E[10] - m.E[12] * m.E[2] * m.E[ 9]; 418 result.E[ 2] = m.E[1] * m.E[ 6] * m.E[15] - m.E[1] * m.E[ 7] * m.E[14] - m.E[5] * m.E[2] * m.E[15] + m.E[5] * m.E[3] * m.E[14] + m.E[13] * m.E[2] * m.E[ 7] - m.E[13] * m.E[3] * m.E[ 6]; 419 result.E[ 6] = -m.E[0] * m.E[ 6] * m.E[15] + m.E[0] * m.E[ 7] * m.E[14] + m.E[4] * m.E[2] * m.E[15] - m.E[4] * m.E[3] * m.E[14] - m.E[12] * m.E[2] * m.E[ 7] + m.E[12] * m.E[3] * m.E[ 6]; 420 result.E[10] = m.E[0] * m.E[ 5] * m.E[15] - m.E[0] * m.E[ 7] * m.E[13] - m.E[4] * m.E[1] * m.E[15] + m.E[4] * m.E[3] * m.E[13] + m.E[12] * m.E[1] * m.E[ 7] - m.E[12] * m.E[3] * m.E[ 5]; 421 result.E[14] = -m.E[0] * m.E[ 5] * m.E[14] + m.E[0] * m.E[ 6] * m.E[13] + m.E[4] * m.E[1] * m.E[14] - m.E[4] * m.E[2] * m.E[13] - m.E[12] * m.E[1] * m.E[ 6] + m.E[12] * m.E[2] * m.E[ 5]; 422 result.E[ 3] = -m.E[1] * m.E[ 6] * m.E[11] + m.E[1] * m.E[ 7] * m.E[10] + m.E[5] * m.E[2] * m.E[11] - m.E[5] * m.E[3] * m.E[10] - m.E[ 9] * m.E[2] * m.E[ 7] + m.E[ 9] * m.E[3] * m.E[ 6]; 423 result.E[ 7] = m.E[0] * m.E[ 6] * m.E[11] - m.E[0] * m.E[ 7] * m.E[10] - m.E[4] * m.E[2] * m.E[11] + m.E[4] * m.E[3] * m.E[10] + m.E[ 8] * m.E[2] * m.E[ 7] - m.E[ 8] * m.E[3] * m.E[ 6]; 424 result.E[11] = -m.E[0] * m.E[ 5] * m.E[11] + m.E[0] * m.E[ 7] * m.E[ 9] + m.E[4] * m.E[1] * m.E[11] - m.E[4] * m.E[3] * m.E[ 9] - m.E[ 8] * m.E[1] * m.E[ 7] + m.E[ 8] * m.E[3] * m.E[ 5]; 425 result.E[15] = m.E[0] * m.E[ 5] * m.E[10] - m.E[0] * m.E[ 6] * m.E[ 9] - m.E[4] * m.E[1] * m.E[10] + m.E[4] * m.E[2] * m.E[ 9] + m.E[ 8] * m.E[1] * m.E[ 6] - m.E[ 8] * m.E[2] * m.E[ 5]; 426 427 f32 determinant = m.E[0] * result.E[0] + m.E[1] * result.E[4] + m.E[2] * result.E[8] + m.E[3] * result.E[12]; 428 determinant = 1.0f / determinant; 429 for(i32 i = 0; i < 16; i++) 430 result.E[i] *= determinant; 431 return result; 432 } 433 434 function m4 435 m4_translation(v3 delta) 436 { 437 m4 result; 438 result.c[0] = (v4){{1, 0, 0, 0}}; 439 result.c[1] = (v4){{0, 1, 0, 0}}; 440 result.c[2] = (v4){{0, 0, 1, 0}}; 441 result.c[3] = (v4){{delta.x, delta.y, delta.z, 1}}; 442 return result; 443 } 444 445 function m4 446 m4_scale(v3 scale) 447 { 448 m4 result; 449 result.c[0] = (v4){{scale.x, 0, 0, 0}}; 450 result.c[1] = (v4){{0, scale.y, 0, 0}}; 451 result.c[2] = (v4){{0, 0, scale.z, 0}}; 452 result.c[3] = (v4){{0, 0, 0, 1}}; 453 return result; 454 } 455 456 function m4 457 m4_rotation_about_z(f32 turns) 458 { 459 f32 sa = sin_f32(turns * 2 * PI); 460 f32 ca = cos_f32(turns * 2 * PI); 461 m4 result; 462 result.c[0] = (v4){{ca, -sa, 0, 0}}; 463 result.c[1] = (v4){{sa, ca, 0, 0}}; 464 result.c[2] = (v4){{0, 0, 1, 0}}; 465 result.c[3] = (v4){{0, 0, 0, 1}}; 466 return result; 467 } 468 469 function m4 470 m4_rotation_about_y(f32 turns) 471 { 472 f32 sa = sin_f32(turns * 2 * PI); 473 f32 ca = cos_f32(turns * 2 * PI); 474 m4 result; 475 result.c[0] = (v4){{ca, 0, -sa, 0}}; 476 result.c[1] = (v4){{0, 1, 0, 0}}; 477 result.c[2] = (v4){{sa, 0, ca, 0}}; 478 result.c[3] = (v4){{0, 0, 0, 1}}; 479 return result; 480 } 481 482 function m4 483 y_aligned_volume_transform(v3 extent, v3 translation, f32 rotation_turns) 484 { 485 m4 T = m4_translation(translation); 486 m4 R = m4_rotation_about_y(rotation_turns); 487 m4 S = m4_scale(extent); 488 m4 result = m4_mul(T, m4_mul(R, S)); 489 return result; 490 } 491 492 function v4 493 m4_mul_v4(m4 a, v4 v) 494 { 495 v4 result; 496 result.x = v4_dot(m4_row(a, 0), v); 497 result.y = v4_dot(m4_row(a, 1), v); 498 result.z = v4_dot(m4_row(a, 2), v); 499 result.w = v4_dot(m4_row(a, 3), v); 500 return result; 501 } 502 503 function m4 504 orthographic_projection(f32 n, f32 f, f32 t, f32 r) 505 { 506 m4 result; 507 f32 a = -2 / (f - n); 508 f32 b = - (f + n) / (f - n); 509 result.c[0] = (v4){{1 / r, 0, 0, 0}}; 510 result.c[1] = (v4){{0, 1 / t, 0, 0}}; 511 result.c[2] = (v4){{0, 0, a, 0}}; 512 result.c[3] = (v4){{0, 0, b, 1}}; 513 return result; 514 } 515 516 function m4 517 perspective_projection(f32 n, f32 f, f32 fov, f32 aspect) 518 { 519 m4 result; 520 f32 t = tan_f32(fov / 2.0f); 521 f32 r = t * aspect; 522 f32 a = -(f + n) / (f - n); 523 f32 b = -2 * f * n / (f - n); 524 result.c[0] = (v4){{1 / r, 0, 0, 0}}; 525 result.c[1] = (v4){{0, 1 / t, 0, 0}}; 526 result.c[2] = (v4){{0, 0, a, -1}}; 527 result.c[3] = (v4){{0, 0, b, 0}}; 528 return result; 529 } 530 531 function m4 532 camera_look_at(v3 camera, v3 point) 533 { 534 v3 orthogonal = {{0, 1.0f, 0}}; 535 v3 normal = v3_normalize(v3_sub(camera, point)); 536 v3 right = cross(orthogonal, normal); 537 v3 up = cross(normal, right); 538 539 v3 translate; 540 camera = v3_sub((v3){0}, camera); 541 translate.x = v3_dot(camera, right); 542 translate.y = v3_dot(camera, up); 543 translate.z = v3_dot(camera, normal); 544 545 m4 result; 546 result.c[0] = (v4){{right.x, up.x, normal.x, 0}}; 547 result.c[1] = (v4){{right.y, up.y, normal.y, 0}}; 548 result.c[2] = (v4){{right.z, up.z, normal.z, 0}}; 549 result.c[3] = (v4){{translate.x, translate.y, translate.z, 1}}; 550 return result; 551 } 552 553 /* NOTE(rnp): adapted from "Essential Mathematics for Games and Interactive Applications" (Verth, Bishop) */ 554 function f32 555 obb_raycast(m4 obb_orientation, v3 obb_size, v3 obb_center, ray r) 556 { 557 v3 p = v3_sub(obb_center, r.origin); 558 v3 X = obb_orientation.c[0].xyz; 559 v3 Y = obb_orientation.c[1].xyz; 560 v3 Z = obb_orientation.c[2].xyz; 561 562 /* NOTE(rnp): projects direction vector onto OBB axis */ 563 v3 f; 564 f.x = v3_dot(X, r.direction); 565 f.y = v3_dot(Y, r.direction); 566 f.z = v3_dot(Z, r.direction); 567 568 /* NOTE(rnp): projects relative vector onto OBB axis */ 569 v3 e; 570 e.x = v3_dot(X, p); 571 e.y = v3_dot(Y, p); 572 e.z = v3_dot(Z, p); 573 574 f32 result = 0; 575 f32 t[6] = {0}; 576 for (i32 i = 0; i < 3; i++) { 577 if (f32_cmp(f.E[i], 0)) { 578 if (-e.E[i] - obb_size.E[i] > 0 || -e.E[i] + obb_size.E[i] < 0) 579 result = -1.0f; 580 f.E[i] = F32_EPSILON; 581 } 582 t[i * 2 + 0] = (e.E[i] + obb_size.E[i]) / f.E[i]; 583 t[i * 2 + 1] = (e.E[i] - obb_size.E[i]) / f.E[i]; 584 } 585 586 if (result != -1) { 587 f32 tmin = MAX(MAX(MIN(t[0], t[1]), MIN(t[2], t[3])), MIN(t[4], t[5])); 588 f32 tmax = MIN(MIN(MAX(t[0], t[1]), MAX(t[2], t[3])), MAX(t[4], t[5])); 589 if (tmax >= 0 && tmin <= tmax) { 590 result = tmin > 0 ? tmin : tmax; 591 } else { 592 result = -1; 593 } 594 } 595 596 return result; 597 } 598 599 function f32 600 complex_filter_first_moment(v2 *filter, i32 length, f32 sampling_frequency) 601 { 602 f32 n = 0, d = 0; 603 for (i32 i = 0; i < length; i++) { 604 f32 t = v2_magnitude_squared(filter[i]); 605 n += (f32)i * t; 606 d += t; 607 } 608 f32 result = n / d / sampling_frequency; 609 return result; 610 } 611 612 function f32 613 real_filter_first_moment(f32 *filter, i32 length, f32 sampling_frequency) 614 { 615 f32 n = 0, d = 0; 616 for (i32 i = 0; i < length; i++) { 617 f32 t = filter[i] * filter[i]; 618 n += (f32)i * t; 619 d += t; 620 } 621 f32 result = n / d / sampling_frequency; 622 return result; 623 } 624 625 function f32 626 tukey_window(f32 t, f32 tapering) 627 { 628 f32 r = tapering; 629 f32 result = 1; 630 if (t < r / 2) result = 0.5f * (1 + cos_f32(2 * PI * (t - r / 2) / r)); 631 if (t >= 1 - r / 2) result = 0.5f * (1 + cos_f32(2 * PI * (t - 1 + r / 2) / r)); 632 return result; 633 } 634 635 /* NOTE(rnp): adapted from "Discrete Time Signal Processing" (Oppenheim) */ 636 function f32 * 637 kaiser_low_pass_filter(Arena *arena, f32 cutoff_frequency, f32 sampling_frequency, f32 beta, i32 length) 638 { 639 f32 *result = push_array(arena, f32, length); 640 f32 wc = 2 * PI * cutoff_frequency / sampling_frequency; 641 f32 a = (f32)length / 2.0f; 642 f32 pi_i0_b = PI * (f32)cephes_i0(beta); 643 644 for (i32 n = 0; n < length; n++) { 645 f32 t = (f32)n - a; 646 f32 impulse = !f32_cmp(t, 0) ? sin_f32(wc * t) / t : wc; 647 t = t / a; 648 f32 window = (f32)cephes_i0(beta * sqrt_f32(1 - t * t)) / pi_i0_b; 649 result[n] = impulse * window; 650 } 651 652 return result; 653 } 654 655 function f32 * 656 rf_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency, 657 i32 length, b32 reverse) 658 { 659 f32 *result = push_array(arena, f32, length); 660 for (i32 i = 0; i < length; i++) { 661 i32 index = reverse? length - 1 - i : i; 662 f32 fc = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length); 663 f32 arg = 2 * PI * fc * (f32)i / sampling_frequency; 664 result[index] = sin_f32(arg) * tukey_window((f32)i / (f32)length, 0.2f); 665 } 666 return result; 667 } 668 669 function v2 * 670 baseband_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency, 671 i32 length, b32 reverse, f32 scale) 672 { 673 v2 *result = push_array(arena, v2, length); 674 f32 conjugate = reverse ? -1 : 1; 675 for (i32 i = 0; i < length; i++) { 676 i32 index = reverse? length - 1 - i : i; 677 f32 fc = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length); 678 f32 arg = 2 * PI * fc * (f32)i / sampling_frequency; 679 v2 sample = {{scale * cos_f32(arg), conjugate * scale * sin_f32(arg)}}; 680 result[index] = v2_scale(sample, tukey_window((f32)i / (f32)length, 0.2f)); 681 } 682 return result; 683 } 684 685 function v4 686 hsv_to_rgb(v4 hsv) 687 { 688 /* f(k(n)) = V - V*S*max(0, min(k, min(4 - k, 1))) 689 * k(n) = fmod((n + H * 6), 6) 690 * (R, G, B) = (f(n = 5), f(n = 3), f(n = 1)) 691 */ 692 alignas(16) f32 nval[4] = {5.0f, 3.0f, 1.0f, 0.0f}; 693 f32x4 n = load_f32x4(nval); 694 f32x4 H = dup_f32x4(hsv.x); 695 f32x4 S = dup_f32x4(hsv.y); 696 f32x4 V = dup_f32x4(hsv.z); 697 f32x4 six = dup_f32x4(6); 698 699 f32x4 t = add_f32x4(n, mul_f32x4(six, H)); 700 f32x4 rem = floor_f32x4(div_f32x4(t, six)); 701 f32x4 k = sub_f32x4(t, mul_f32x4(rem, six)); 702 703 t = min_f32x4(sub_f32x4(dup_f32x4(4), k), dup_f32x4(1)); 704 t = max_f32x4(dup_f32x4(0), min_f32x4(k, t)); 705 t = mul_f32x4(t, mul_f32x4(S, V)); 706 707 v4 rgba; 708 store_f32x4(rgba.E, sub_f32x4(V, t)); 709 rgba.a = hsv.a; 710 return rgba; 711 } 712 713 function f32 714 ease_in_out_cubic(f32 t) 715 { 716 f32 result; 717 if (t < 0.5f) { 718 result = 4.0f * t * t * t; 719 } else { 720 t = -2.0f * t + 2.0f; 721 result = 1.0f - t * t * t / 2.0f; 722 } 723 return result; 724 } 725 726 function f32 727 ease_in_out_quartic(f32 t) 728 { 729 f32 result; 730 if (t < 0.5f) { 731 result = 8.0f * t * t * t * t; 732 } else { 733 t = -2.0f * t + 2.0f; 734 result = 1.0f - t * t * t * t / 2.0f; 735 } 736 return result; 737 }