viffy 0.1.5

///! Dumb math modules, use something better like `nalgebra_glm` for
///! better performance - use this only for quick serialization where
///! calculations are not the bottleneck

#[derive(Clone, Copy, Debug, PartialOrd, PartialEq)]
pub struct Vector<const N: usize>(pub [f32; N]);
impl<const N: usize> Default for Vector<N> {
	fn default() -> Self {
		Self(core::array::from_fn(|_| 0.0))
	}
}
impl<const N: usize> Vector<N> {
	pub fn square(self) -> Self {
		Self(core::array::from_fn(|i| self.0[i] * self.0[i]))
	}
	pub fn magnitude(self) -> f32 {
		let mut m = 0.0f32;
		for e in self.0 {
			m += e * e;
		}
		m.sqrt()
	}
	pub fn normalize(self) -> Self {
		let mag = self.magnitude();
		Self(core::array::from_fn(|i| self.0[i] / mag))
	}
	pub fn reduce(self) -> f32 {
		let mut m = 0.0f32;
		for e in self.0 {
			m += e;
		}
		m
	}
	pub fn cross(self, other: Self) -> Self {
		Self(core::array::from_fn(|i|
			self.0[(i + 1) % N] * other.0[(i + 2) % N] - self.0[(i + 2) % N] * other.0[(i + 1) % N]
		))
	}
	pub fn x(self) -> f32 { self.0[0] }
	pub fn y(self) -> f32 { self.0[1] }
	pub fn z(self) -> f32 { self.0[2] }
	pub fn w(self) -> f32 { self.0[3] }
}
impl<const N: usize> std::ops::Index<usize> for Vector<N> {
	type Output = f32;
	fn index(&self, index: usize) -> &Self::Output {
		if index < self.0.len() {
			&self.0[index]
		} else {
			panic!("")
		}
	}
}
impl<const N: usize> std::ops::IndexMut<usize> for Vector<N> {
	fn index_mut<'a>(&'a mut self, index: usize) -> &'a mut f32 {
		if index < self.0.len() {
			&mut self.0[index]
		} else {
			panic!("")
		}
	}
}
impl<const N: usize> std::ops::Add for Vector<N> {
	type Output = Vector<N>;
	fn add(self, rhs: Self) -> Self::Output {
		Self(core::array::from_fn(|i| self.0[i] + rhs.0[i]))
	}
}
impl<const N: usize> std::ops::Sub for Vector<N> {
	type Output = Vector<N>;
	fn sub(self, rhs: Self) -> Self::Output {
		Self(core::array::from_fn(|i| self.0[i] - rhs.0[i]))
	}
}
impl<const N: usize> std::ops::Mul for Vector<N> {
	type Output = Vector<N>;
	fn mul(self, rhs: Self) -> Self::Output {
		Self(core::array::from_fn(|i| self.0[i] * rhs.0[i]))
	}
}
impl<const N: usize> std::ops::Div for Vector<N> {
	type Output = Vector<N>;
	fn div(self, rhs: Self) -> Self::Output {
		Self(core::array::from_fn(|i| self.0[i] / rhs.0[i]))
	}
}

pub type Vector2 = Vector<2>;
pub type Vector3 = Vector<3>;
pub type Vector4 = Vector<4>;

#[derive(Default, Clone, Copy, Debug, PartialOrd, PartialEq)]
pub struct Matrix4x4(pub [[f32; 4]; 4]);
impl Matrix4x4 {
	pub fn size(&self) -> usize {
		4
	}
	pub fn identity() -> Self {
		Self([
			[1.0, 0.0, 0.0, 0.0],
			[0.0, 1.0, 0.0, 0.0],
			[0.0, 0.0, 1.0, 0.0],
			[0.0, 0.0, 0.0, 1.0],
		])
	}
	pub fn rotate(self, a: Vector3, r: f32) -> Self {
		let c1 = 1.0 - r.cos();
		self * Self([
			[a.x() * a.x() * c1 + r.cos(),			a.x() * a.y() * c1 - a.z() * r.sin(),	a.x() * a.z() * c1 + a.y() * r.sin(), 0.0],
			[a.x() * a.y() * c1 + a.z() * r.sin(),	a.y() * a.y() * c1 + r.cos(),			a.y() * a.z() * c1 - a.x() * r.sin(), 0.0],
			[a.x() * a.z() * c1 - a.y() * r.sin(),	a.y() * a.z() * c1 + a.x() * r.sin(),	a.z() * a.z() * c1 + r.cos(), 0.0],
			[0.0, 0.0, 0.0, 1.0],
		])
	}
	pub fn scale(self, a: Vector3) -> Self {
		self * Self([
			[a.x(), 0.0, 0.0, 0.0],
			[0.0, a.y(), 0.0, 0.0],
			[0.0, 0.0, a.z(), 0.0],
			[0.0, 0.0, 0.0, 1.0],
		])
	}
	pub fn translate(self, a: Vector3) -> Self {
		self * Self([
			[1.0, 0.0, 0.0, a.x()],
			[0.0, 1.0, 0.0, a.y()],
			[0.0, 0.0, 1.0, a.z()],
			[0.0, 0.0, 0.0, 1.0],
		])
	}
	pub fn transpose(self) -> Self {
		let mut m = Self::default();
		for i in 0..4 {
			for j in 0..4 {
				m.0[j][i] = self.0[i][j];
			}
		}
		m
	}
	pub fn ortho(left: f32, right: f32, bottom: f32, top: f32, near: f32, far: f32) -> Self {
		let tx = -(right + left) / (right - left);
		let ty = -(top + bottom) / (top - bottom);
		let tz = -(far + near) / (far - near);
		Self([
			[2.0 / (right - left), 0.0, 0.0, tx],
			[0.0, 2.0 / (top - bottom), 0.0, ty],
			[0.0, 0.0, -2.0 / (far - near), tz],
			[0.0, 0.0, 0.0, 1.0],
		])
	}
	pub fn perspective(fov: f32, aspx: f32, near: f32, far: f32) -> Self {
		let f = 1.0 / fov.tan();
		Self([
			[f / aspx, 0.0, 0.0, 0.0],
			[0.0, f, 0.0, 0.0],
			[0.0, 0.0, (far + near) / (near - far), (2.0 * far * near) / (near - far)],
			[0.0, 0.0, -1.0, 0.0],
		])
	}
	// First tuple element is L (lower triangular), second is U (upper triangular)
	pub fn lu_decomposition(self) -> (Self, Self) {
		let mut s = self;
		let mut r = Self::identity();
		for j in 0..4 {
			for i in (j + 1)..4 { // pivot (i)
				let pivot = Self::f32_div(s.0[i][j], s.0[j][j]);
				let row_a = s.0[j];
				let row_b = s.0[i];
				r.0[i][j] = pivot;
				for k in j..4 {
					s.0[i][k] = row_b[k] - pivot * row_a[k];
				}
			}
		}
		(r, s)
	}
	pub fn lu_doolittle_decomposition(self) -> Self {
		let mut lu = Matrix4x4::default();
		for i in 0..4 {
			for j in i..4 {
				let mut sum = 0.0;
				for k in 0..i {
					sum += lu.0[i][k] * lu.0[k][j];
				}
				lu.0[i][j] = self.0[i][j] - sum;
			}
			for j in (i + 1)..4 {
				let mut sum = 0.0;
				for k in 0..i {
					sum += lu.0[j][k] * lu.0[k][i];
				}
				lu.0[j][i] = Self::f32_div(1.0, lu.0[i][i]) * (self.0[j][i] - sum);
			}
		}
		lu
	}
	fn f32_div(x: f32, y: f32) -> f32 {
		[x / y, 0.0][(y == 0.0) as usize]
	}
	// Solves a 3D vector using LU decomposition, solves the equation
	// of the form "Ax = b"
	pub fn solve_lu_vector3(self, b: Vector4) -> Vector4 {
		let lu = self.lu_doolittle_decomposition();
		let mut y = Vector4::default();
		for i in 0..4 {
			let mut sum = 0.0;
			for j in 0..i {
				sum += lu.0[i][j] * y[j];
			}
			y[i] = b[i] - sum;
		}
		let mut x = Vector4::default();
		for i in (0..4).into_iter().rev() {
			let mut sum = 0.0;
			for j in (i + 1)..4 {
				sum += lu.0[i][j] * x[j];
			}
			x[i] = Self::f32_div(1.0, lu.0[i][i]) * (y[i] - sum);
		}
		x
	}

	pub fn look_at(eye: Vector3, center: Vector3, up: Vector3) -> Self {
		let big_f = center - eye;
		let f = big_f.normalize();
		let up_d = up.normalize();
		let s = f.cross(up_d);
		let u = s.normalize().cross(f);
		Matrix4x4([
			[s[0], s[1], s[2], 0.0f32],
			[u[0], u[1], u[2], 0.0f32],
			[-f[0], -f[1], -f[2], 0.0f32],
			[0.0f32, 0.0f32, 0.0f32, 1.0f32],
		])
	}
}
impl std::ops::Add for Matrix4x4 {
	type Output = Matrix4x4;
	fn add(self, rhs: Self) -> Self::Output {
		let mut m = Self::default();
		for i in 0..4 {
			for j in 0..4 {
				m.0[i][j] = self.0[i][j] + rhs.0[i][j];
			}
		}
		m
	}
}
impl std::ops::Sub for Matrix4x4 {
	type Output = Matrix4x4;
	fn sub(self, rhs: Self) -> Self::Output {
		let mut m = Self::default();
		for i in 0..4 {
			for j in 0..4 {
				m.0[i][j] = self.0[i][j] - rhs.0[i][j];
			}
		}
		m
	}
}
impl std::ops::Mul for Matrix4x4 {
	type Output = Matrix4x4;
	fn mul(self, rhs: Self) -> Self::Output {
		let mut m = Self::default();
		for i in 0..4 {
			for k in 0..4 {
				for j in 0..4 {
					m.0[i][j] += self.0[i][k] * rhs.0[k][j];
				}
			}
		}
		m
	}
}

/// Deterministic exp(x)
pub fn exp(x: f32) -> f32 {
	if x > -176.059402 {
		let a: f32 = ((1 << 22) as f32) / 0.69314718055994530942;
		let b: i32 = 127 * (1 << 23);
		let r = (a * x) as i32;
		let s = b + r;
		let t = b - r;
		f32::from_bits(s as u32) / f32::from_bits(t as u32)
	} else {
		0.0
	}
}

/// Deterministic floor
fn floor(x: f32) -> f32 {
	let i = x as i32;
	i as f32 - (i as f32 > x) as i32 as f32
}

/// Deterministic sin
fn sin(mut x: f32) -> f32 {
	use core::f32::consts::PI as PI;
	let px = (floor(x / PI) as u32 & 0x1) << 31;
	let qx = floor(1.0 + x / PI);
	x = x + PI * (1.0 - qx);
	let p = 0.036783 * x * x * x * x
		- 0.23111 * x * x * x
		+ 0.048506 * x * x
		+ 0.98811 * x;
	unsafe {
		let mut q = f32::to_bits(p);
		q |= px;
		f32::from_bits(q)
	}
}