ogl_beamforming

math.c (24736B)

---

 /* See LICENSE for license details. */
 #include "external/cephes.c"
 
 function void
 fill_kronecker_sub_matrix_f16(f16 *out, i32 out_stride, f16 scale, f16 *b, iv2 b_dim)
 {
 	for (i32 i = 0; i < b_dim.y; i++) {
 		for (i32 j = 0; j < b_dim.x; j += 4, b += 4) {
 			out[j + 0] = scale * b[0];
 			out[j + 1] = scale * b[1];
 			out[j + 2] = scale * b[2];
 			out[j + 3] = scale * b[3];
 		}
 		out += out_stride;
 	}
 }
 
 /* NOTE: this won't check for valid space/etc and assumes row major order */
 function void
 kronecker_product_f16(f16 *out, f16 *a, iv2 a_dim, f16 *b, iv2 b_dim)
 {
 	iv2 out_dim = {{a_dim.x * b_dim.x, a_dim.y * b_dim.y}};
 	assert(out_dim.y % 4 == 0);
 	for (i32 i = 0; i < a_dim.y; i++) {
 		f16 *vout = out;
 		for (i32 j = 0; j < a_dim.x; j++, a++) {
 			fill_kronecker_sub_matrix_f16(vout, out_dim.y, *a, b, b_dim);
 			vout += b_dim.y;
 		}
 		out += out_dim.y * b_dim.x;
 	}
 }
 
 /* NOTE/TODO: to support even more hadamard sizes use the Paley construction */
 function f16 *
 make_hadamard_transpose(Arena *a, i32 dim, b32 row_major)
 {
 	read_only local_persist	f16 hadamard_12_12_transpose[] = {
 		1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 		1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1,
 		1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,
 		1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,  1,
 		1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1,
 		1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,
 		1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,
 		1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1,
 		1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1, -1,
 		1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,  1,
 		1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1, -1,
 		1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1, -1,
 	};
 
 	read_only local_persist f16 hadamard_20_20_transpose[] = {
 		1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 		1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1,
 		1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1,
 		1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,
 		1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,
 		1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1,
 		1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1,
 		1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1,
 		1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1,
 		1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,
 		1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1,
 		1,  1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,
 		1, -1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1,
 		1,  1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,
 		1,  1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,
 		1,  1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,
 		1,  1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,
 		1, -1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1,
 		1, -1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1,
 		1,  1, -1, -1,  1,  1, -1, -1, -1, -1,  1, -1,  1, -1,  1,  1,  1,  1, -1, -1,
 	};
 
 
 	f16 *result = 0;
 
 	i32 order          = dim;
 	b32 power_of_2     = IsPowerOfTwo(dim);
 	b32 multiple_of_12 = dim % 12 == 0;
 	b32 multiple_of_20 = dim % 20 == 0;
 	iz elements        = dim * dim;
 
 	i32 base_dim = 0;
 	if (power_of_2) {
 		base_dim  = dim;
 	} else if (multiple_of_20 && IsPowerOfTwo(dim / 20)) {
 		base_dim  = 20;
 		dim      /= 20;
 	} else if (multiple_of_12 && IsPowerOfTwo(dim / 12)) {
 		base_dim  = 12;
 		dim      /= 12;
 	}
 
 	if (power_of_2 && base_dim && arena_capacity(a, f16) >= elements * (1 + (dim != base_dim))) {
 		result = push_array(a, f16, elements);
 
 		Arena tmp = *a;
 		f16 *m = dim == base_dim ? result : push_array(&tmp, f16, elements);
 
 		#define IND(i, j) ((i) * dim + (j))
 		m[0] = 1;
 		for (i32 k = 1; k < dim; k *= 2) {
 			for (i32 i = 0; i < k; i++) {
 				for (i32 j = 0; j < k; j++) {
 					f16 val = m[IND(i, j)];
 					m[IND(i + k, j)]     =  val;
 					m[IND(i, j + k)]     =  val;
 					m[IND(i + k, j + k)] = -val;
 				}
 			}
 		}
 		#undef IND
 
 		f16 *m2 = 0;
 		iv2 m2_dim;
 		switch (base_dim) {
 		case 12:{ m2 = hadamard_12_12_transpose; m2_dim = (iv2){{12, 12}}; }break;
 		case 20:{ m2 = hadamard_20_20_transpose; m2_dim = (iv2){{20, 20}}; }break;
 		}
 		if (m2) kronecker_product_f16(result, m, (iv2){{dim, dim}}, m2, m2_dim);
 	}
 
 	if (row_major) {
 		for (i32 r = 0; r < order; r++)
 			for (i32 c = 0; c < order; c++)
 				swap(result[r * order + c], result[c * order + r]);
 	}
 
 	return result;
 }
 
 function b32
 u128_equal(u128 a, u128 b)
 {
 	b32 result = a.U64[0] == b.U64[0] && a.U64[1] == b.U64[1];
 	return result;
 }
 
 function RangeU64
 subrange_n_from_n_m_count(u64 n, u64 n_count, u64 m)
 {
 	assert(n < n_count);
 
 	u64 per_lane            = m / n_count;
 	u64 leftover            = m - per_lane * n_count;
 	u64 leftovers_before_n  = MIN(leftover, n);
 	u64 base_index          = n * per_lane + leftovers_before_n;
 	u64 one_past_last_index = base_index + per_lane + ((n < leftover) ? 1 : 0);
 
 	RangeU64 result = {base_index, one_past_last_index};
 	return result;
 }
 
 function b32
 iv2_equal(iv2 a, iv2 b)
 {
 	b32 result = a.x == b.x && a.y == b.y;
 	return result;
 }
 
 function b32
 iv3_equal(iv3 a, iv3 b)
 {
 	b32 result = a.x == b.x && a.y == b.y && a.z == b.z;
 	return result;
 }
 
 function i32
 iv3_dimension(iv3 points)
 {
 	i32 result = (points.x > 1) + (points.y > 1) + (points.z > 1);
 	return result;
 }
 
 function v2
 clamp_v2_rect(v2 v, Rect r)
 {
 	v2 result = v;
 	result.x = CLAMP(v.x, r.pos.x, r.pos.x + r.size.x);
 	result.y = CLAMP(v.y, r.pos.y, r.pos.y + r.size.y);
 	return result;
 }
 
 function v2
 v2_from_iv2(iv2 v)
 {
 	v2 result;
 	result.E[0] = (f32)v.E[0];
 	result.E[1] = (f32)v.E[1];
 	return result;
 }
 
 function v2
 v2_abs(v2 a)
 {
 	v2 result;
 	result.x = Abs(a.x);
 	result.y = Abs(a.y);
 	return result;
 }
 
 function v2
 v2_scale(v2 a, f32 scale)
 {
 	v2 result;
 	result.x = a.x * scale;
 	result.y = a.y * scale;
 	return result;
 }
 
 function v2
 v2_add(v2 a, v2 b)
 {
 	v2 result;
 	result.x = a.x + b.x;
 	result.y = a.y + b.y;
 	return result;
 }
 
 function v2
 v2_sub(v2 a, v2 b)
 {
 	v2 result = v2_add(a, v2_scale(b, -1.0f));
 	return result;
 }
 
 function v2
 v2_mul(v2 a, v2 b)
 {
 	v2 result;
 	result.x = a.x * b.x;
 	result.y = a.y * b.y;
 	return result;
 }
 
 function v2
 v2_div(v2 a, v2 b)
 {
 	v2 result;
 	result.x = a.x / b.x;
 	result.y = a.y / b.y;
 	return result;
 }
 
 function v2
 v2_floor(v2 a)
 {
 	v2 result;
 	result.x = (f32)((i32)a.x);
 	result.y = (f32)((i32)a.y);
 	return result;
 }
 
 function f32
 v2_magnitude_squared(v2 a)
 {
 	f32 result = a.x * a.x + a.y * a.y;
 	return result;
 }
 
 function f32
 v2_magnitude(v2 a)
 {
 	f32 result = sqrt_f32(a.x * a.x + a.y * a.y);
 	return result;
 }
 
 function v3
 cross(v3 a, v3 b)
 {
 	v3 result;
 	result.x = a.y * b.z - a.z * b.y;
 	result.y = a.z * b.x - a.x * b.z;
 	result.z = a.x * b.y - a.y * b.x;
 	return result;
 }
 
 function v3
 v3_from_iv3(iv3 v)
 {
 	v3 result;
 	result.E[0] = (f32)v.E[0];
 	result.E[1] = (f32)v.E[1];
 	result.E[2] = (f32)v.E[2];
 	return result;
 }
 
 function v3
 v3_abs(v3 a)
 {
 	v3 result;
 	result.x = Abs(a.x);
 	result.y = Abs(a.y);
 	result.z = Abs(a.z);
 	return result;
 }
 
 function v3
 v3_scale(v3 a, f32 scale)
 {
 	v3 result;
 	result.x = scale * a.x;
 	result.y = scale * a.y;
 	result.z = scale * a.z;
 	return result;
 }
 
 function v3
 v3_add(v3 a, v3 b)
 {
 	v3 result;
 	result.x = a.x + b.x;
 	result.y = a.y + b.y;
 	result.z = a.z + b.z;
 	return result;
 }
 
 function v3
 v3_sub(v3 a, v3 b)
 {
 	v3 result = v3_add(a, v3_scale(b, -1.0f));
 	return result;
 }
 
 function v3
 v3_div(v3 a, v3 b)
 {
 	v3 result;
 	result.x = a.x / b.x;
 	result.y = a.y / b.y;
 	result.z = a.z / b.z;
 	return result;
 }
 
 function f32
 v3_dot(v3 a, v3 b)
 {
 	f32 result = a.x * b.x + a.y * b.y + a.z * b.z;
 	return result;
 }
 
 function f32
 v3_magnitude_squared(v3 a)
 {
 	f32 result = v3_dot(a, a);
 	return result;
 }
 
 function f32
 v3_magnitude(v3 a)
 {
 	f32 result = sqrt_f32(v3_dot(a, a));
 	return result;
 }
 
 function v3
 v3_normalize(v3 a)
 {
 	v3 result = v3_scale(a, 1.0f / v3_magnitude(a));
 	return result;
 }
 
 function v4
 v4_scale(v4 a, f32 scale)
 {
 	v4 result;
 	result.x = scale * a.x;
 	result.y = scale * a.y;
 	result.z = scale * a.z;
 	result.w = scale * a.w;
 	return result;
 }
 
 function v4
 v4_add(v4 a, v4 b)
 {
 	v4 result;
 	result.x = a.x + b.x;
 	result.y = a.y + b.y;
 	result.z = a.z + b.z;
 	result.w = a.w + b.w;
 	return result;
 }
 
 function v4
 v4_sub(v4 a, v4 b)
 {
 	v4 result = v4_add(a, v4_scale(b, -1));
 	return result;
 }
 
 function f32
 v4_dot(v4 a, v4 b)
 {
 	f32 result = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
 	return result;
 }
 
 function v4
 v4_lerp(v4 a, v4 b, f32 t)
 {
 	v4 result = v4_add(a, v4_scale(v4_sub(b, a), t));
 	return result;
 }
 
 function b32
 m4_equal(m4 a, m4 b)
 {
 	b32 result = 1;
 	for EachElement(a.E, it)
 		result &= f32_equal(a.E[it], b.E[it]);
 	return result;
 }
 
 function m4
 m4_identity(void)
 {
 	m4 result;
 	result.c[0] = (v4){{1, 0, 0, 0}};
 	result.c[1] = (v4){{0, 1, 0, 0}};
 	result.c[2] = (v4){{0, 0, 1, 0}};
 	result.c[3] = (v4){{0, 0, 0, 1}};
 	return result;
 }
 
 function v4
 m4_row(m4 a, u32 row)
 {
 	v4 result;
 	result.E[0] = a.c[0].E[row];
 	result.E[1] = a.c[1].E[row];
 	result.E[2] = a.c[2].E[row];
 	result.E[3] = a.c[3].E[row];
 	return result;
 }
 
 function m4
 m4_mul(m4 a, m4 b)
 {
 	m4 result;
 	for (u32 i = 0; i < 4; i++) {
 		for (u32 j = 0; j < 4; j++) {
 			result.c[i].E[j] = v4_dot(m4_row(a, j), b.c[i]);
 		}
 	}
 	return result;
 }
 
 /* NOTE(rnp): based on:
  * https://web.archive.org/web/20131215123403/ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf
  * TODO(rnp): redo with SIMD as given in the link (but need to rewrite for column-major)
  */
 function m4
 m4_inverse(m4 m)
 {
 	m4 result;
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 	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];
 
 	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];
 	determinant = 1.0f / determinant;
 	for(i32 i = 0; i < 16; i++)
 		result.E[i] *= determinant;
 	return result;
 }
 
 function m4
 m4_translation(v3 delta)
 {
 	m4 result;
 	result.c[0] = (v4){{1, 0, 0, 0}};
 	result.c[1] = (v4){{0, 1, 0, 0}};
 	result.c[2] = (v4){{0, 0, 1, 0}};
 	result.c[3] = (v4){{delta.x, delta.y, delta.z, 1}};
 	return result;
 }
 
 function m4
 m4_scale(v3 scale)
 {
 	m4 result;
 	result.c[0] = (v4){{scale.x, 0,       0,       0}};
 	result.c[1] = (v4){{0,       scale.y, 0,       0}};
 	result.c[2] = (v4){{0,       0,       scale.z, 0}};
 	result.c[3] = (v4){{0,       0,       0,       1}};
 	return result;
 }
 
 function m4
 m4_rotation_about_axis(v3 axis, f32 turns)
 {
 	assert(f32_equal(v3_magnitude_squared(axis), 1.0f));
 	f32 sa  = sin_f32(turns * 2 * PI);
 	f32 ca  = cos_f32(turns * 2 * PI);
 	f32 mca = 1.0f - ca;
 
 	f32 x = axis.x, x2 = x * x;
 	f32 y = axis.y, y2 = y * y;
 	f32 z = axis.z, z2 = z * z;
 
 	m4 result;
 	result.c[0] = (v4){{ca + mca * x2,        mca * x * y - sa * z, mca * x * z + sa * y, 0}};
 	result.c[1] = (v4){{mca * x * y + sa * z, ca + mca * y2,        mca * y * z - sa * x, 0}};
 	result.c[2] = (v4){{mca * x * z - sa * y, mca * y * z + sa * x, ca + mca * z2,        0}};
 	result.c[3] = (v4){{0, 0, 0, 1}};
 	return result;
 }
 
 function m4
 m4_rotation_about_y(f32 turns)
 {
 	m4 result = m4_rotation_about_axis((v3){.y = 1.0f}, turns);
 	return result;
 }
 
 function m4
 y_aligned_volume_transform(v3 extent, v3 translation, f32 rotation_turns)
 {
 	m4 T = m4_translation(translation);
 	m4 R = m4_rotation_about_axis((v3){.y = 1.0f}, rotation_turns);
 	m4 S = m4_scale(extent);
 	m4 result = m4_mul(T, m4_mul(R, S));
 	return result;
 }
 
 function v4
 m4_mul_v4(m4 a, v4 v)
 {
 	v4 result;
 	result.x = v4_dot(m4_row(a, 0), v);
 	result.y = v4_dot(m4_row(a, 1), v);
 	result.z = v4_dot(m4_row(a, 2), v);
 	result.w = v4_dot(m4_row(a, 3), v);
 	return result;
 }
 
 function v3
 m4_mul_v3(m4 a, v3 v)
 {
 	v3 result = m4_mul_v4(a, (v4){{v.x, v.y, v.z, 1.0f}}).xyz;
 	return result;
 }
 
 function m4
 orthographic_projection(f32 n, f32 f, f32 t, f32 r)
 {
 	m4 result;
 	f32 a = -2 / (f - n);
 	f32 b = - (f + n) / (f - n);
 	result.c[0] = (v4){{1 / r, 0,     0,  0}};
 	result.c[1] = (v4){{0,     1 / t, 0,  0}};
 	result.c[2] = (v4){{0,     0,     a,  0}};
 	result.c[3] = (v4){{0,     0,     b,  1}};
 	return result;
 }
 
 function m4
 perspective_projection(f32 n, f32 f, f32 fov, f32 aspect)
 {
 	m4 result;
 	f32 t = tan_f32(fov / 2.0f);
 	f32 r = t * aspect;
 	f32 a = -(f + n) / (f - n);
 	f32 b = -2 * f * n / (f - n);
 	result.c[0] = (v4){{1 / r, 0,     0,  0}};
 	result.c[1] = (v4){{0,     1 / t, 0,  0}};
 	result.c[2] = (v4){{0,     0,     a, -1}};
 	result.c[3] = (v4){{0,     0,     b,  0}};
 	return result;
 }
 
 function m4
 camera_look_at(v3 camera, v3 point)
 {
 	v3 orthogonal = {{0, 1.0f, 0}};
 	v3 normal     = v3_normalize(v3_sub(camera, point));
 	v3 right      = cross(orthogonal, normal);
 	v3 up         = cross(normal,     right);
 
 	v3 translate;
 	camera      = v3_sub((v3){0}, camera);
 	translate.x = v3_dot(camera, right);
 	translate.y = v3_dot(camera, up);
 	translate.z = v3_dot(camera, normal);
 
 	m4 result;
 	result.c[0] = (v4){{right.x,     up.x,        normal.x,    0}};
 	result.c[1] = (v4){{right.y,     up.y,        normal.y,    0}};
 	result.c[2] = (v4){{right.z,     up.z,        normal.z,    0}};
 	result.c[3] = (v4){{translate.x, translate.y, translate.z, 1}};
 	return result;
 }
 
 /* NOTE(rnp): adapted from "Essential Mathematics for Games and Interactive Applications" (Verth, Bishop) */
 function f32
 obb_raycast(m4 obb_orientation, v3 obb_size, v3 obb_center, ray r)
 {
 	v3 p = v3_sub(obb_center, r.origin);
 	v3 X = obb_orientation.c[0].xyz;
 	v3 Y = obb_orientation.c[1].xyz;
 	v3 Z = obb_orientation.c[2].xyz;
 
 	/* NOTE(rnp): projects direction vector onto OBB axis */
 	v3 f;
 	f.x = v3_dot(X, r.direction);
 	f.y = v3_dot(Y, r.direction);
 	f.z = v3_dot(Z, r.direction);
 
 	/* NOTE(rnp): projects relative vector onto OBB axis */
 	v3 e;
 	e.x = v3_dot(X, p);
 	e.y = v3_dot(Y, p);
 	e.z = v3_dot(Z, p);
 
 	f32 result = 0;
 	f32 t[6] = {0};
 	for (i32 i = 0; i < 3; i++) {
 		if (f32_equal(f.E[i], 0)) {
 			if (-e.E[i] - obb_size.E[i] > 0 || -e.E[i] + obb_size.E[i] < 0)
 				result = -1.0f;
 			f.E[i] = F32_EPSILON;
 		}
 		t[i * 2 + 0] = (e.E[i] + obb_size.E[i]) / f.E[i];
 		t[i * 2 + 1] = (e.E[i] - obb_size.E[i]) / f.E[i];
 	}
 
 	if (result != -1) {
 		f32 tmin = MAX(MAX(MIN(t[0], t[1]), MIN(t[2], t[3])), MIN(t[4], t[5]));
 		f32 tmax = MIN(MIN(MAX(t[0], t[1]), MAX(t[2], t[3])), MAX(t[4], t[5]));
 		if (tmax >= 0 && tmin <= tmax) {
 			result = tmin > 0 ? tmin : tmax;
 		} else {
 			result = -1;
 		}
 	}
 
 	return result;
 }
 
 function f32
 complex_filter_first_moment(v2 *filter, i32 length, f32 sampling_frequency)
 {
 	f32 n = 0, d = 0;
 	for (i32 i = 0; i < length; i++) {
 		f32 t = v2_magnitude_squared(filter[i]);
 		n += (f32)i * t;
 		d += t;
 	}
 	f32 result = n / d / sampling_frequency;
 	return result;
 }
 
 function f32
 real_filter_first_moment(f32 *filter, i32 length, f32 sampling_frequency)
 {
 	f32 n = 0, d = 0;
 	for (i32 i = 0; i < length; i++) {
 		f32 t = filter[i] * filter[i];
 		n += (f32)i * t;
 		d += t;
 	}
 	f32 result = n / d / sampling_frequency;
 	return result;
 }
 
 function f32
 tukey_window(f32 t, f32 tapering)
 {
 	f32 r = tapering;
 	f32 result = 1;
 	if (t < r / 2)      result = 0.5f * (1 + cos_f32(2 * PI * (t - r / 2)     / r));
 	if (t >= 1 - r / 2) result = 0.5f * (1 + cos_f32(2 * PI * (t - 1 + r / 2) / r));
 	return result;
 }
 
 /* NOTE(rnp): adapted from "Discrete Time Signal Processing" (Oppenheim) */
 function f32 *
 kaiser_low_pass_filter(Arena *arena, f32 cutoff_frequency, f32 sampling_frequency, f32 beta, i32 length)
 {
 	f32 *result = push_array(arena, f32, length);
 	f32 wc      = 2 * PI * cutoff_frequency / sampling_frequency;
 	f32 a       = (f32)length / 2.0f;
 	f32 pi_i0_b = PI * (f32)cephes_i0(beta);
 
 	for (i32 n = 0; n < length; n++) {
 		f32 t       = (f32)n - a;
 		f32 impulse = !f32_equal(t, 0) ? sin_f32(wc * t) / t : wc;
 		t           = t / a;
 		f32 window  = (f32)cephes_i0(beta * sqrt_f32(1 - t * t)) / pi_i0_b;
 		result[n]   = impulse * window;
 	}
 
 	return result;
 }
 
 function f32 *
 rf_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency,
          i32 length, b32 reverse)
 {
 	f32 *result = push_array(arena, f32, length);
 	for (i32 i = 0; i < length; i++) {
 		i32 index = reverse? length - 1 - i : i;
 		f32 fc    = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length);
 		f32 arg   = 2 * PI * fc * (f32)i / sampling_frequency;
 		result[index] = sin_f32(arg) * tukey_window((f32)i / (f32)length, 0.2f);
 	}
 	return result;
 }
 
 function v2 *
 baseband_chirp(Arena *arena, f32 min_frequency, f32 max_frequency, f32 sampling_frequency,
                i32 length, b32 reverse, f32 scale)
 {
 	v2 *result    = push_array(arena, v2, length);
 	f32 conjugate = reverse ? -1 : 1;
 	for (i32 i = 0; i < length; i++) {
 		i32 index = reverse? length - 1 - i : i;
 		f32 fc    = min_frequency + (f32)i * (max_frequency - min_frequency) / (2 * (f32)length);
 		f32 arg   = 2 * PI * fc * (f32)i / sampling_frequency;
 		v2 sample = {{scale * cos_f32(arg), conjugate * scale * sin_f32(arg)}};
 		result[index] = v2_scale(sample, tukey_window((f32)i / (f32)length, 0.2f));
 	}
 	return result;
 }
 
 function iv3
 das_output_dimension(iv3 points)
 {
 	iv3 result;
 	result.x = Max(points.x, 1);
 	result.y = Max(points.y, 1);
 	result.z = Max(points.z, 1);
 
 	switch (iv3_dimension(result)) {
 	case 1:{
 		if (result.y > 1) result.x = result.y;
 		if (result.z > 1) result.x = result.z;
 		result.y = result.z = 1;
 	}break;
 
 	case 2:{
 		if (result.x > 1) {
 			if (result.z > 1) result.y = result.z;
 		} else {
 			result.x = result.z;
 		}
 		result.z = 1;
 	}break;
 
 	case 3:{}break;
 
 	InvalidDefaultCase;
 	}
 
 	return result;
 }
 
 function m4
 das_transform_1d(v3 p1, v3 p2)
 {
 	v3 extent = v3_sub(p2, p1);
 	m4 result = {
 		.c[0] = (v4){{extent.x, extent.y, extent.z, 0.0f}},
 		.c[1] = (v4){{0.0f, 0.0f, 0.0f, 0.0f}},
 		.c[2] = (v4){{0.0f, 0.0f, 0.0f, 0.0f}},
 		.c[3] = (v4){{p1.x, p1.y, p1.z, 1.0f}},
 	};
 	return result;
 }
 
 function m4
 das_transform_2d_with_normal(v3 normal, v2 min_coordinate, v2 max_coordinate, f32 offset)
 {
 	v3 U = {{0, 1.0f, 0}};
 	if (f32_equal(v3_dot(U, normal), 1.0f))
 		U = (v3){{1.0f, 0, 0}};
 
 	v3 N = normal;
 	v3 V = cross(U, N);
 
 	v3 min = v3_add(v3_scale(U, min_coordinate.x), v3_scale(V, min_coordinate.y));
 	v3 max = v3_add(v3_scale(U, max_coordinate.x), v3_scale(V, max_coordinate.y));
 
 	v3 extent = v3_sub(max, min);
 	U = v3_scale(U, v3_dot(U, extent));
 	V = v3_scale(V, v3_dot(V, extent));
 
 	v3 t = v3_add(v3_scale(N, offset), min);
 
 	m4 result;
 	result.c[0] = (v4){{U.x,  U.y,  U.z,  0.0f}};
 	result.c[1] = (v4){{V.x,  V.y,  V.z,  0.0f}};
 	result.c[2] = (v4){{N.x,  N.y,  N.z,  0.0f}};
 	result.c[3] = (v4){{t.x,  t.y,  t.z,  1.0f}};
 
 	return result;
 }
 
 function m4
 das_transform_2d_xz(v2 min_coordinate, v2 max_coordinate, f32 y_off)
 {
 	m4 result = das_transform_2d_with_normal((v3){.y = 1.0f}, min_coordinate, max_coordinate, y_off);
 	return result;
 }
 
 function m4
 das_transform_2d_yz(v2 min_coordinate, v2 max_coordinate, f32 x_off)
 {
 	m4 result = das_transform_2d_with_normal((v3){.x = 1.0f}, min_coordinate, max_coordinate, x_off);
 	return result;
 }
 
 function m4
 das_transform_2d_xy(v2 min_coordinate, v2 max_coordinate, f32 z_off)
 {
 	m4 result = das_transform_2d_with_normal((v3){.z = 1.0f}, min_coordinate, max_coordinate, z_off);
 	return result;
 }
 
 function m4
 das_transform_3d(v3 min_coordinate, v3 max_coordinate)
 {
 	v3 extent = v3_sub(max_coordinate, min_coordinate);
 	m4 result;
 	result.c[0] = (v4){{extent.x,         0.0f,             0.0f,             0.0f}};
 	result.c[1] = (v4){{0.0f,             extent.y,         0.0f,             0.0f}};
 	result.c[2] = (v4){{0.0f,             0.0f,             extent.z,         0.0f}};
 	result.c[3] = (v4){{min_coordinate.x, min_coordinate.y, min_coordinate.z, 1.0f}};
 	return result;
 }
 
 function m4
 das_transform(v3 min_coordinate, v3 max_coordinate, iv3 *points)
 {
 	m4 result;
 
 	*points = das_output_dimension(*points);
 
 	switch (iv3_dimension(*points)) {
 	case 1:{result = das_transform_1d(      min_coordinate,     max_coordinate);    }break;
 	case 2:{result = das_transform_2d_xz(XY(min_coordinate), XY(max_coordinate), 0);}break;
 	case 3:{result = das_transform_3d(      min_coordinate,     max_coordinate);    }break;
 	}
 
 	return result;
 }
 
 function v2
 plane_uv(v3 point, v3 U, v3 V)
 {
 	v2 result;
 	result.x = v3_dot(U, point) / v3_dot(U, U);
 	result.y = v3_dot(V, point) / v3_dot(V, V);
 	return result;
 }
 
 function v4
 hsv_to_rgb(v4 hsv)
 {
 	/* f(k(n))   = V - V*S*max(0, min(k, min(4 - k, 1)))
 	 * k(n)      = fmod((n + H * 6), 6)
 	 * (R, G, B) = (f(n = 5), f(n = 3), f(n = 1))
 	 */
 	alignas(16) f32 nval[4] = {5.0f, 3.0f, 1.0f, 0.0f};
 	f32x4 n   = load_f32x4(nval);
 	f32x4 H   = dup_f32x4(hsv.x);
 	f32x4 S   = dup_f32x4(hsv.y);
 	f32x4 V   = dup_f32x4(hsv.z);
 	f32x4 six = dup_f32x4(6);
 
 	f32x4 t   = add_f32x4(n, mul_f32x4(six, H));
 	f32x4 rem = floor_f32x4(div_f32x4(t, six));
 	f32x4 k   = sub_f32x4(t, mul_f32x4(rem, six));
 
 	t = min_f32x4(sub_f32x4(dup_f32x4(4), k), dup_f32x4(1));
 	t = max_f32x4(dup_f32x4(0), min_f32x4(k, t));
 	t = mul_f32x4(t, mul_f32x4(S, V));
 
 	v4 rgba;
 	store_f32x4(rgba.E, sub_f32x4(V, t));
 	rgba.a = hsv.a;
 	return rgba;
 }
 
 function f32
 ease_in_out_cubic(f32 t)
 {
 	f32 result;
 	if (t < 0.5f) {
 		result = 4.0f * t * t * t;
 	} else {
 		t      = -2.0f * t + 2.0f;
 		result =  1.0f - t * t * t / 2.0f;
 	}
 	return result;
 }
 
 function f32
 ease_in_out_quartic(f32 t)
 {
 	f32 result;
 	if (t < 0.5f) {
 		result = 8.0f * t * t * t * t;
 	} else {
 		t      = -2.0f * t + 2.0f;
 		result =  1.0f - t * t * t * t / 2.0f;
 	}
 	return result;
 }