vec_arith/

lib.rs

// MIT License
//
// Copyright (c) 2019 Tobias Pfeiffer
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#![feature(clamp)]
use std::{
	ops::*,
	iter::FromIterator
};

pub use self::{
	vec2::*,
	vec3::*,
	vec4::*,
	quat32::*,
	quat64::*,
	mat2::*,
	mat3::*,
	mat4::*,
	bezier1::*,
	bezier2::*,
	bezier3::*,
	sphere::*,
	cuboid::*,
	plane::*,
	cylinder::*,
	triangle::*,
	ops::*,
	poly::*,
	misc::*,
	sdf::*
};

pub mod vec2 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Vec2<T>(pub T, pub T);
	
	impl<T> Vec2<T> {
		#[inline]
		pub fn map<F: FnMut(T) -> U, U>(self, mut f: F) -> Vec2<U> {
			Vec2(f(self.0),
				 f(self.1))
		}
	}
	
	impl<T: Copy> Vec2<T> {
		#[inline] pub fn xx(self) -> Vec2<T> { Vec2(self.0, self.0) }
		#[inline] pub fn xy(self) -> Vec2<T> { Vec2(self.0, self.1) }
		#[inline] pub fn yx(self) -> Vec2<T> { Vec2(self.1, self.0) }
		#[inline] pub fn yy(self) -> Vec2<T> { Vec2(self.1, self.1) }
	}
	
	impl<T: Default> Default for Vec2<T> {
		fn default() -> Self {
			Self(T::default(),
				 T::default())
		}
	}
	
	impl<T: Copy> From<T> for Vec2<T> {
		fn from(v: T) -> Self {
			Self(v, v)
		}
	}
	
	impl<T: Copy> From<[T; 2]> for Vec2<T> {
		fn from(v: [T; 2]) -> Self {
			Self(v[0], v[1])
		}
	}
	
	impl<T> From<(T, T)> for Vec2<T> {
		fn from(v: (T, T)) -> Self {
			Self(v.0, v.1)
		}
	}
	
	impl Vec2<f32> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1))
		}
		
		pub fn powi(self, rhs: Vec2<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1))
		}
		
		pub fn len2(self) -> f32 {
			self.dot(self)
		}
		
		pub fn len(self) -> f32 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f32) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance2(self, other: Self) -> f32 {
			(other - self).len2()
		}
		
		pub fn distance(self, other: Self) -> f32 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f32 {
			(self.dot(other) / (self.len() * other.len())).acos()
		}
	}
	
	impl Vec2<f64> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1))
		}
		
		pub fn powi(self, rhs: Vec2<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1))
		}
		
		pub fn len2(self) -> f64 {
			self.dot(self)
		}
		
		pub fn len(self) -> f64 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f64) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance2(self, other: Self) -> f64 {
			(other - self).len2()
		}
		
		pub fn distance(self, other: Self) -> f64 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f64 {
			self.dot(other) / (self.len() * other.len())
		}
	}
	
	impl Vec2<bool> {
		pub fn and(self) -> bool {
			self.0 && self.1
		}
		
		pub fn or(self) -> bool {
			self.0 || self.1
		}
		
		pub fn xor(self) -> bool {
			self.0 ^ self.1
		}
	}
	
	impl<T> Index<usize> for Vec2<T> {
		type Output = T;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Vec2<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				_ => panic!()
			}
		}
	}
	
	impl<T: PartialOrd + Copy> Vec2<T> {
		pub fn min(self, other: Self) -> Self {
			Self(if self.0 < other.0 { self.0 } else { other.0 },
				 if self.1 < other.1 { self.1 } else { other.1 })
		}
		
		pub fn mins(self, other: T) -> Self {
			Self(if self.0 < other { self.0 } else { other },
				 if self.1 < other { self.1 } else { other })
		}
		
		pub fn max(self, other: Self) -> Self {
			Self(if self.0 > other.0 { self.0 } else { other.0 },
				 if self.1 > other.1 { self.1 } else { other.1 })
		}
		
		pub fn maxs(self, other: T) -> Self {
			Self(if self.0 > other { self.0 } else { other },
				 if self.1 > other { self.1 } else { other })
		}
		
		pub fn clamp(self, lower: Self, upper: Self) -> Self {
			Self(if self.0 < lower.0 { lower.0 } else if self.0 > upper.0 { upper.0 } else { self.0 },
				 if self.1 < lower.1 { lower.1 } else if self.1 > upper.1 { upper.1 } else { self.1 })
		}
		
		pub fn clamps(self, lower: T, upper: T) -> Self {
			Self(if self.0 < lower { lower } else if self.0 > upper { upper } else { self.0 },
				 if self.1 < lower { lower } else if self.1 > upper { upper } else { self.1 })
		}
		
		pub fn clampvs(self, lower: Self, upper: T) -> Self {
			Self(if self.0 < lower.0 { lower.0 } else if self.0 > upper { upper } else { self.0 },
				 if self.1 < lower.1 { lower.1 } else if self.1 > upper { upper } else { self.1 })
		}
		
		pub fn clampsv(self, lower: T, upper: Self) -> Self {
			Self(if self.0 < lower { lower } else if self.0 > upper.0 { upper.0 } else { self.0 },
				 if self.1 < lower { lower } else if self.1 > upper.1 { upper.1 } else { self.1 })
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Copy> Vec2<T> {
		pub fn dot(self, other: Self) -> T {
			self.0 * other.0 + self.1 * other.1
		}
		
		pub fn mix(self, a: Self, b: Self) -> Self {
			Self(a.0 + (b.0 - a.0) * self.0,
				 a.1 + (b.1 - a.1) * self.1)
		}
		
		pub fn mixs(self, a: T, b: T) -> Self {
			Self(a + (b - a) * self.0,
				 a + (b - a) * self.1)
		}
	}
	
	impl<T: Add> Add for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn add(self, rhs: Self) -> Self::Output {
			Vec2(self.0 + rhs.0,
				 self.1 + rhs.1)
		}
	}
	
	impl<T: Sub> Sub for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn sub(self, rhs: Self) -> Self::Output {
			Vec2(self.0 - rhs.0,
				 self.1 - rhs.1)
		}
	}
	
	impl<T: Mul> Mul for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn mul(self, rhs: Self) -> Self::Output {
			Vec2(self.0 * rhs.0,
				 self.1 * rhs.1)
		}
	}
	
	impl<T: Mul + Copy> Mul<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn mul(self, rhs: T) -> Self::Output {
			Vec2(self.0 * rhs,
				 self.1 * rhs)
		}
	}
	
	impl<T: Div> Div for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn div(self, rhs: Self) -> Self::Output {
			Vec2(self.0 / rhs.0,
				 self.1 / rhs.1)
		}
	}
	
	impl<T: Div + Copy> Div<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn div(self, rhs: T) -> Self::Output {
			Vec2(self.0 / rhs,
				 self.1 / rhs)
		}
	}
	
	impl<T: Rem> Rem for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn rem(self, rhs: Self) -> Self::Output {
			Vec2(self.0 % rhs.0,
				 self.1 % rhs.1)
		}
	}
	
	impl<T: Rem + Copy> Rem<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn rem(self, rhs: T) -> Self::Output {
			Vec2(self.0 % rhs,
				 self.1 % rhs)
		}
	}
	
	impl<T: Neg> Neg for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn neg(self) -> Self::Output {
			Vec2(-self.0,
				 -self.1)
		}
	}
	
	impl<T: Not> Not for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn not(self) -> Self::Output {
			Vec2(!self.0,
				 !self.1)
		}
	}
	
	impl<T: BitAnd> BitAnd for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitand(self, rhs: Self) -> Self::Output {
			Vec2(self.0 & rhs.0,
				 self.1 & rhs.1)
		}
	}
	
	impl<T: BitAnd + Copy> BitAnd<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitand(self, rhs: T) -> Self::Output {
			Vec2(self.0 & rhs,
				 self.1 & rhs)
		}
	}
	
	impl<T: BitOr> BitOr for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitor(self, rhs: Self) -> Self::Output {
			Vec2(self.0 | rhs.0,
				 self.1 | rhs.1)
		}
	}
	
	impl<T: BitOr + Copy> BitOr<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitor(self, rhs: T) -> Self::Output {
			Vec2(self.0 | rhs,
				 self.1 | rhs)
		}
	}
	
	impl<T: BitXor> BitXor for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitxor(self, rhs: Self) -> Self::Output {
			Vec2(self.0 ^ rhs.0,
				 self.1 ^ rhs.1)
		}
	}
	
	impl<T: BitXor + Copy> BitXor<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn bitxor(self, rhs: T) -> Self::Output {
			Vec2(self.0 ^ rhs,
				 self.1 ^ rhs)
		}
	}
	
	impl<T: Shl> Shl for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn shl(self, rhs: Self) -> Self::Output {
			Vec2(self.0 << rhs.0,
				 self.1 << rhs.1)
		}
	}
	
	impl<T: Shl + Copy> Shl<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn shl(self, rhs: T) -> Self::Output {
			Vec2(self.0 << rhs,
				 self.1 << rhs)
		}
	}
	
	impl<T: Shr> Shr for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn shr(self, rhs: Self) -> Self::Output {
			Vec2(self.0 >> rhs.0,
				 self.1 >> rhs.1)
		}
	}
	
	impl<T: Shr + Copy> Shr<T> for Vec2<T> {
		type Output = Vec2<T::Output>;
		
		fn shr(self, rhs: T) -> Self::Output {
			Vec2(self.0 >> rhs,
				 self.1 >> rhs)
		}
	}
	
	impl<T: AddAssign> AddAssign for Vec2<T> {
		fn add_assign(&mut self, rhs: Self) {
			self.0 += rhs.0;
			self.1 += rhs.1;
		}
	}
	
	impl<T: SubAssign> SubAssign for Vec2<T> {
		fn sub_assign(&mut self, rhs: Self) {
			self.0 -= rhs.0;
			self.1 -= rhs.1;
		}
	}
	
	impl<T: MulAssign> MulAssign for Vec2<T> {
		fn mul_assign(&mut self, rhs: Self) {
			self.0 *= rhs.0;
			self.1 *= rhs.1;
		}
	}
	
	impl<T: MulAssign + Copy> MulAssign<T> for Vec2<T> {
		fn mul_assign(&mut self, rhs: T) {
			self.0 *= rhs;
			self.1 *= rhs;
		}
	}
	
	impl<T: DivAssign> DivAssign for Vec2<T> {
		fn div_assign(&mut self, rhs: Self) {
			self.0 /= rhs.0;
			self.1 /= rhs.1;
		}
	}
	
	impl<T: DivAssign + Copy> DivAssign<T> for Vec2<T> {
		fn div_assign(&mut self, rhs: T) {
			self.0 /= rhs;
			self.1 /= rhs;
		}
	}
	
	impl<T: RemAssign> RemAssign for Vec2<T> {
		fn rem_assign(&mut self, rhs: Self) {
			self.0 %= rhs.0;
			self.1 %= rhs.1;
		}
	}
	
	impl<T: RemAssign + Copy> RemAssign<T> for Vec2<T> {
		fn rem_assign(&mut self, rhs: T) {
			self.0 %= rhs;
			self.1 %= rhs;
		}
	}
	
	impl<T: BitAndAssign> BitAndAssign for Vec2<T> {
		fn bitand_assign(&mut self, rhs: Self) {
			self.0 &= rhs.0;
			self.1 &= rhs.1;
		}
	}
	
	impl<T: BitAndAssign + Copy> BitAndAssign<T> for Vec2<T> {
		fn bitand_assign(&mut self, rhs: T) {
			self.0 &= rhs;
			self.1 &= rhs;
		}
	}
	
	impl<T: BitOrAssign> BitOrAssign for Vec2<T> {
		fn bitor_assign(&mut self, rhs: Self) {
			self.0 |= rhs.0;
			self.1 |= rhs.1;
		}
	}
	
	impl<T: BitOrAssign + Copy> BitOrAssign<T> for Vec2<T> {
		fn bitor_assign(&mut self, rhs: T) {
			self.0 |= rhs;
			self.1 |= rhs;
		}
	}
	
	impl<T: BitXorAssign> BitXorAssign for Vec2<T> {
		fn bitxor_assign(&mut self, rhs: Self) {
			self.0 ^= rhs.0;
			self.1 ^= rhs.1;
		}
	}
	
	impl<T: BitXorAssign + Copy> BitXorAssign<T> for Vec2<T> {
		fn bitxor_assign(&mut self, rhs: T) {
			self.0 ^= rhs;
			self.1 ^= rhs;
		}
	}
	
	impl<T: ShlAssign> ShlAssign for Vec2<T> {
		fn shl_assign(&mut self, rhs: Self) {
			self.0 <<= rhs.0;
			self.1 <<= rhs.1;
		}
	}
	
	impl<T: ShlAssign + Copy> ShlAssign<T> for Vec2<T> {
		fn shl_assign(&mut self, rhs: T) {
			self.0 <<= rhs;
			self.1 <<= rhs;
		}
	}
	
	impl<T: ShrAssign> ShrAssign for Vec2<T> {
		fn shr_assign(&mut self, rhs: Self) {
			self.0 >>= rhs.0;
			self.1 >>= rhs.1;
		}
	}
	
	impl<T: ShrAssign + Copy> ShrAssign<T> for Vec2<T> {
		fn shr_assign(&mut self, rhs: T) {
			self.0 >>= rhs;
			self.1 >>= rhs;
		}
	}
	
	impl<A> FromIterator<A> for Vec2<A> {
		fn from_iter<T: IntoIterator<Item = A>>(iter: T) -> Self {
			let mut iter = iter.into_iter();
			Vec2(iter.next().unwrap(),
				 iter.next().unwrap())
		}
	}
	
	pub struct Iter<T>(Vec2<T>, usize);
	
	impl<T: Clone> Iterator for Iter<T> {
		type Item = T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some((self.0).0.clone()),
				2 => Some((self.0).1.clone()),
				_ => None,
			}
		}
	}
	
	impl<T: Clone> IntoIterator for Vec2<T> {
		type Item = T;
		type IntoIter = Iter<T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			Iter(self, 0)
		}
	}
	
	pub struct IterRef<'a, T>(&'a Vec2<T>, usize);
	
	impl<'a, T> Iterator for IterRef<'a, T> {
		type Item = &'a T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(&(self.0).0),
				2 => Some(&(self.0).1),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a Vec2<T> {
		type Item = &'a T;
		type IntoIter = IterRef<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterRef(self, 0)
		}
	}
	
	pub struct IterMut<'a, T>(&'a mut Vec2<T>, usize);
	
	impl<'a, T> Iterator for IterMut<'a, T> {
		type Item = &'a mut T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(unsafe { ::std::mem::transmute(&mut (self.0).0) }),
				2 => Some(unsafe { ::std::mem::transmute(&mut (self.0).1) }),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a mut Vec2<T> {
		type Item = &'a mut T;
		type IntoIter = IterMut<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterMut(self, 0)
		}
	}
}

pub mod vec3 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Vec3<T>(pub T, pub T, pub T);
	
	impl<T> Vec3<T> {
		#[inline]
		pub fn map<F: FnMut(T) -> U, U>(self, mut f: F) -> Vec3<U> {
			Vec3(f(self.0),
				 f(self.1),
				 f(self.2))
		}
	}
	
	impl<T: Copy> Vec3<T> {
		#[inline] pub fn xx(self) -> Vec2<T> { Vec2(self.0, self.0) }
		#[inline] pub fn xy(self) -> Vec2<T> { Vec2(self.0, self.1) }
		#[inline] pub fn xz(self) -> Vec2<T> { Vec2(self.0, self.2) }
		#[inline] pub fn yx(self) -> Vec2<T> { Vec2(self.1, self.0) }
		#[inline] pub fn yy(self) -> Vec2<T> { Vec2(self.1, self.1) }
		#[inline] pub fn yz(self) -> Vec2<T> { Vec2(self.1, self.2) }
		#[inline] pub fn zx(self) -> Vec2<T> { Vec2(self.2, self.0) }
		#[inline] pub fn zy(self) -> Vec2<T> { Vec2(self.2, self.1) }
		#[inline] pub fn zz(self) -> Vec2<T> { Vec2(self.2, self.2) }
		#[inline] pub fn xxx(self) -> Vec3<T> { Vec3(self.0, self.0, self.0) }
		#[inline] pub fn xxy(self) -> Vec3<T> { Vec3(self.0, self.0, self.1) }
		#[inline] pub fn xxz(self) -> Vec3<T> { Vec3(self.0, self.0, self.2) }
		#[inline] pub fn xyx(self) -> Vec3<T> { Vec3(self.0, self.1, self.0) }
		#[inline] pub fn xyy(self) -> Vec3<T> { Vec3(self.0, self.1, self.1) }
		#[inline] pub fn xyz(self) -> Vec3<T> { Vec3(self.0, self.1, self.2) }
		#[inline] pub fn xzx(self) -> Vec3<T> { Vec3(self.0, self.2, self.0) }
		#[inline] pub fn xzy(self) -> Vec3<T> { Vec3(self.0, self.2, self.1) }
		#[inline] pub fn xzz(self) -> Vec3<T> { Vec3(self.0, self.2, self.2) }
		#[inline] pub fn yxx(self) -> Vec3<T> { Vec3(self.1, self.0, self.0) }
		#[inline] pub fn yxy(self) -> Vec3<T> { Vec3(self.1, self.0, self.1) }
		#[inline] pub fn yxz(self) -> Vec3<T> { Vec3(self.1, self.0, self.2) }
		#[inline] pub fn yyx(self) -> Vec3<T> { Vec3(self.1, self.1, self.0) }
		#[inline] pub fn yyy(self) -> Vec3<T> { Vec3(self.1, self.1, self.1) }
		#[inline] pub fn yyz(self) -> Vec3<T> { Vec3(self.1, self.1, self.2) }
		#[inline] pub fn yzx(self) -> Vec3<T> { Vec3(self.1, self.2, self.0) }
		#[inline] pub fn yzy(self) -> Vec3<T> { Vec3(self.1, self.2, self.1) }
		#[inline] pub fn yzz(self) -> Vec3<T> { Vec3(self.1, self.2, self.2) }
		#[inline] pub fn zxx(self) -> Vec3<T> { Vec3(self.2, self.0, self.0) }
		#[inline] pub fn zxy(self) -> Vec3<T> { Vec3(self.2, self.0, self.1) }
		#[inline] pub fn zxz(self) -> Vec3<T> { Vec3(self.2, self.0, self.2) }
		#[inline] pub fn zyx(self) -> Vec3<T> { Vec3(self.2, self.1, self.0) }
		#[inline] pub fn zyy(self) -> Vec3<T> { Vec3(self.2, self.1, self.1) }
		#[inline] pub fn zyz(self) -> Vec3<T> { Vec3(self.2, self.1, self.2) }
		#[inline] pub fn zzx(self) -> Vec3<T> { Vec3(self.2, self.2, self.0) }
		#[inline] pub fn zzy(self) -> Vec3<T> { Vec3(self.2, self.2, self.1) }
		#[inline] pub fn zzz(self) -> Vec3<T> { Vec3(self.2, self.2, self.2) }
	}
	
	impl<T: Default> Default for Vec3<T> {
		fn default() -> Self {
			Self(T::default(),
				 T::default(),
				 T::default())
		}
	}
	
	impl<T: Copy> From<T> for Vec3<T> {
		fn from(v: T) -> Self {
			Vec3(v, v, v)
		}
	}
	
	impl<T: Copy> From<[T; 3]> for Vec3<T> {
		fn from(v: [T; 3]) -> Self {
			Vec3(v[0], v[1], v[2])
		}
	}
	
	impl<T> From<(T, T, T)> for Vec3<T> {
		fn from(v: (T, T, T)) -> Self {
			Vec3(v.0, v.1, v.2)
		}
	}
	
	impl<T> From<(Vec2<T>, T)> for Vec3<T> {
		fn from(v: (Vec2<T>, T)) -> Self {
			Vec3((v.0).0, (v.0).1, v.1)
		}
	}
	
	impl<T> From<(T, Vec2<T>)> for Vec3<T> {
		fn from(v: (T, Vec2<T>)) -> Self {
			Vec3(v.0, (v.1).0, (v.1).1)
		}
	}
	
	impl Vec3<f32> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs(),
				 self.2.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1),
				 self.2.powf(rhs.2))
		}
		
		pub fn powi(self, rhs: Vec3<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1),
				 self.2.powi(rhs.2))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt(),
				 self.2.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt(),
				 self.2.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1),
				 self.2.log(rhs.2))
		}
		
		pub fn len2(self) -> f32 {
			self.dot(self)
		}
		
		pub fn len(self) -> f32 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f32) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance2(self, other: Self) -> f32 {
			(other - self).len2()
		}
		
		pub fn distance(self, other: Self) -> f32 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f32 {
			self.dot(other) / (self.len() * other.len())
		}
	}
	
	impl Vec3<f64> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs(),
				 self.2.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1),
				 self.2.powf(rhs.2))
		}
		
		pub fn powi(self, rhs: Vec3<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1),
				 self.2.powi(rhs.2))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt(),
				 self.2.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt(),
				 self.2.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1),
				 self.2.log(rhs.2))
		}
		
		pub fn len2(self) -> f64 {
			self.dot(self)
		}
		
		pub fn len(self) -> f64 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f64) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance(self, other: Self) -> f64 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f64 {
			(self.dot(other) / (self.len() * other.len())).acos()
		}
	}
	
	impl Vec3<bool> {
		pub fn and(self) -> bool {
			self.0 && self.1 && self.2
		}
		
		pub fn or(self) -> bool {
			self.0 || self.1 || self.2
		}
		
		pub fn xor(self) -> bool {
			self.0 ^ self.1 ^ self.2
		}
	}
	
	impl<T> Index<usize> for Vec3<T> {
		type Output = T;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				2 => &self.2,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Vec3<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				2 => &mut self.2,
				_ => panic!()
			}
		}
	}
	
	impl<T: PartialOrd + Copy> Vec3<T> {
		pub fn min(self, other: Self) -> Self {
			Self(if self.0 < other.0 { self.0 } else { other.0 },
				 if self.1 < other.1 { self.1 } else { other.1 },
				 if self.2 < other.2 { self.2 } else { other.2 })
		}
		
		pub fn mins(self, other: T) -> Self {
			Self(if self.0 < other { self.0 } else { other },
				 if self.1 < other { self.1 } else { other },
				 if self.2 < other { self.2 } else { other })
		}
		
		pub fn max(self, other: Self) -> Self {
			Self(if self.0 > other.0 { self.0 } else { other.0 },
				 if self.1 > other.1 { self.1 } else { other.1 },
				 if self.2 > other.2 { self.2 } else { other.2 })
		}
		
		pub fn maxs(self, other: T) -> Self {
			Self(if self.0 > other { self.0 } else { other },
				 if self.1 > other { self.1 } else { other },
				 if self.2 > other { self.2 } else { other })
		}
		
		pub fn clamp(self, lower: Self, upper: Self) -> Self {
			Self(if self.0 < lower.0 { lower.0 } else if self.0 > upper.0 { upper.0 } else { self.0 },
				 if self.1 < lower.1 { lower.1 } else if self.1 > upper.1 { upper.1 } else { self.1 },
				 if self.2 < lower.2 { lower.2 } else if self.2 > upper.2 { upper.2 } else { self.2 })
		}
		
		pub fn clamps(self, lower: T, upper: T) -> Self {
			Self(if self.0 < lower { lower } else if self.0 > upper { upper } else { self.0 },
				 if self.1 < lower { lower } else if self.1 > upper { upper } else { self.1 },
				 if self.2 < lower { lower } else if self.2 > upper { upper } else { self.2 })
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Copy> Vec3<T> {
		pub fn dot(self, other: Self) -> T {
			self.0 * other.0 + self.1 * other.1 + self.2 * other.2
		}
		
		pub fn cross(self, other: Self) -> Self {
			Self(self.1 * other.2 - self.2 * other.1,
				 self.2 * other.0 - self.0 * other.2,
				 self.0 * other.1 - self.1 * other.0)
		}
		
		pub fn mix(self, a: Self, b: Self) -> Self {
			Self(a.0 + (b.0 - a.0) * self.0,
				 a.1 + (b.1 - a.1) * self.1,
				 a.2 + (b.2 - a.2) * self.2)
		}
	}
	
	impl<T: Add> Add for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn add(self, rhs: Self) -> Self::Output {
			Vec3(self.0 + rhs.0,
				 self.1 + rhs.1,
				 self.2 + rhs.2)
		}
	}
	
	impl<T: Sub> Sub for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn sub(self, rhs: Self) -> Self::Output {
			Vec3(self.0 - rhs.0,
				 self.1 - rhs.1,
				 self.2 - rhs.2)
		}
	}
	
	impl<T: Mul> Mul for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn mul(self, rhs: Self) -> Self::Output {
			Vec3(self.0 * rhs.0,
				 self.1 * rhs.1,
				 self.2 * rhs.2)
		}
	}
	
	impl<T: Mul + Copy> Mul<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn mul(self, rhs: T) -> Self::Output {
			Vec3(self.0 * rhs,
				 self.1 * rhs,
				 self.2 * rhs)
		}
	}
	
	impl<T: Div> Div for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn div(self, rhs: Self) -> Self::Output {
			Vec3(self.0 / rhs.0,
				 self.1 / rhs.1,
				 self.2 / rhs.2)
		}
	}
	
	impl<T: Div + Copy> Div<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn div(self, rhs: T) -> Self::Output {
			Vec3(self.0 / rhs,
				 self.1 / rhs,
				 self.2 / rhs)
		}
	}
	
	impl<T: Rem> Rem for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn rem(self, rhs: Self) -> Self::Output {
			Vec3(self.0 % rhs.0,
				 self.1 % rhs.1,
				 self.2 % rhs.2)
		}
	}
	
	impl<T: Rem + Copy> Rem<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn rem(self, rhs: T) -> Self::Output {
			Vec3(self.0 % rhs,
				 self.1 % rhs,
				 self.2 % rhs)
		}
	}
	
	impl<T: Neg> Neg for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn neg(self) -> Self::Output {
			Vec3(-self.0,
				 -self.1,
				 -self.2)
		}
	}
	
	impl<T: Not> Not for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn not(self) -> Self::Output {
			Vec3(!self.0,
				 !self.1,
				 !self.2)
		}
	}
	
	impl<T: BitAnd> BitAnd for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitand(self, rhs: Self) -> Self::Output {
			Vec3(self.0 & rhs.0,
				 self.1 & rhs.1,
				 self.2 & rhs.2)
		}
	}
	
	impl<T: BitAnd + Copy> BitAnd<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitand(self, rhs: T) -> Self::Output {
			Vec3(self.0 & rhs,
				 self.1 & rhs,
				 self.2 & rhs)
		}
	}
	
	impl<T: BitOr> BitOr for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitor(self, rhs: Self) -> Self::Output {
			Vec3(self.0 | rhs.0,
				 self.1 | rhs.1,
				 self.2 | rhs.2)
		}
	}
	
	impl<T: BitOr + Copy> BitOr<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitor(self, rhs: T) -> Self::Output {
			Vec3(self.0 | rhs,
				 self.1 | rhs,
				 self.2 | rhs)
		}
	}
	
	impl<T: BitXor> BitXor for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitxor(self, rhs: Self) -> Self::Output {
			Vec3(self.0 ^ rhs.0,
				 self.1 ^ rhs.1,
				 self.2 ^ rhs.2)
		}
	}
	
	impl<T: BitXor + Copy> BitXor<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn bitxor(self, rhs: T) -> Self::Output {
			Vec3(self.0 ^ rhs,
				 self.1 ^ rhs,
				 self.2 ^ rhs)
		}
	}
	
	impl<T: Shl> Shl for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn shl(self, rhs: Self) -> Self::Output {
			Vec3(self.0 << rhs.0,
				 self.1 << rhs.1,
				 self.2 << rhs.2)
		}
	}
	
	impl<T: Shl + Copy> Shl<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn shl(self, rhs: T) -> Self::Output {
			Vec3(self.0 << rhs,
				 self.1 << rhs,
				 self.2 << rhs)
		}
	}
	
	impl<T: Shr> Shr for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn shr(self, rhs: Self) -> Self::Output {
			Vec3(self.0 >> rhs.0,
				 self.1 >> rhs.1,
				 self.2 >> rhs.2)
		}
	}
	
	impl<T: Shr + Copy> Shr<T> for Vec3<T> {
		type Output = Vec3<T::Output>;
		
		fn shr(self, rhs: T) -> Self::Output {
			Vec3(self.0 >> rhs,
				 self.1 >> rhs,
				 self.2 >> rhs)
		}
	}
	
	impl<T: AddAssign> AddAssign for Vec3<T> {
		fn add_assign(&mut self, rhs: Self) {
			self.0 += rhs.0;
			self.1 += rhs.1;
			self.2 += rhs.2;
		}
	}
	
	impl<T: SubAssign> SubAssign for Vec3<T> {
		fn sub_assign(&mut self, rhs: Self) {
			self.0 -= rhs.0;
			self.1 -= rhs.1;
			self.2 -= rhs.2;
		}
	}
	
	impl<T: MulAssign> MulAssign for Vec3<T> {
		fn mul_assign(&mut self, rhs: Self) {
			self.0 *= rhs.0;
			self.1 *= rhs.1;
			self.2 *= rhs.2;
		}
	}
	
	impl<T: MulAssign + Copy> MulAssign<T> for Vec3<T> {
		fn mul_assign(&mut self, rhs: T) {
			self.0 *= rhs;
			self.1 *= rhs;
			self.2 *= rhs;
		}
	}
	
	impl<T: DivAssign> DivAssign for Vec3<T> {
		fn div_assign(&mut self, rhs: Self) {
			self.0 /= rhs.0;
			self.1 /= rhs.1;
			self.2 /= rhs.2;
		}
	}
	
	impl<T: DivAssign + Copy> DivAssign<T> for Vec3<T> {
		fn div_assign(&mut self, rhs: T) {
			self.0 /= rhs;
			self.1 /= rhs;
			self.2 /= rhs;
		}
	}
	
	impl<T: RemAssign> RemAssign for Vec3<T> {
		fn rem_assign(&mut self, rhs: Self) {
			self.0 %= rhs.0;
			self.1 %= rhs.1;
			self.2 %= rhs.2;
		}
	}
	
	impl<T: RemAssign + Copy> RemAssign<T> for Vec3<T> {
		fn rem_assign(&mut self, rhs: T) {
			self.0 %= rhs;
			self.1 %= rhs;
			self.2 %= rhs;
		}
	}
	
	impl<T: BitAndAssign> BitAndAssign for Vec3<T> {
		fn bitand_assign(&mut self, rhs: Self) {
			self.0 &= rhs.0;
			self.1 &= rhs.1;
			self.2 &= rhs.2;
		}
	}
	
	impl<T: BitAndAssign + Copy> BitAndAssign<T> for Vec3<T> {
		fn bitand_assign(&mut self, rhs: T) {
			self.0 &= rhs;
			self.1 &= rhs;
			self.2 &= rhs;
		}
	}
	
	impl<T: BitOrAssign> BitOrAssign for Vec3<T> {
		fn bitor_assign(&mut self, rhs: Self) {
			self.0 |= rhs.0;
			self.1 |= rhs.1;
			self.2 |= rhs.2;
		}
	}
	
	impl<T: BitOrAssign + Copy> BitOrAssign<T> for Vec3<T> {
		fn bitor_assign(&mut self, rhs: T) {
			self.0 |= rhs;
			self.1 |= rhs;
			self.2 |= rhs;
		}
	}
	
	impl<T: BitXorAssign> BitXorAssign for Vec3<T> {
		fn bitxor_assign(&mut self, rhs: Self) {
			self.0 ^= rhs.0;
			self.1 ^= rhs.1;
			self.2 ^= rhs.2;
		}
	}
	
	impl<T: BitXorAssign + Copy> BitXorAssign<T> for Vec3<T> {
		fn bitxor_assign(&mut self, rhs: T) {
			self.0 ^= rhs;
			self.1 ^= rhs;
			self.2 ^= rhs;
		}
	}
	
	impl<T: ShlAssign> ShlAssign for Vec3<T> {
		fn shl_assign(&mut self, rhs: Self) {
			self.0 <<= rhs.0;
			self.1 <<= rhs.1;
			self.2 <<= rhs.2;
		}
	}
	
	impl<T: ShlAssign + Copy> ShlAssign<T> for Vec3<T> {
		fn shl_assign(&mut self, rhs: T) {
			self.0 <<= rhs;
			self.1 <<= rhs;
			self.2 <<= rhs;
		}
	}
	
	impl<T: ShrAssign> ShrAssign for Vec3<T> {
		fn shr_assign(&mut self, rhs: Self) {
			self.0 >>= rhs.0;
			self.1 >>= rhs.1;
			self.2 >>= rhs.2;
		}
	}
	
	impl<T: ShrAssign + Copy> ShrAssign<T> for Vec3<T> {
		fn shr_assign(&mut self, rhs: T) {
			self.0 >>= rhs;
			self.1 >>= rhs;
			self.2 >>= rhs;
		}
	}
	
	impl<A> FromIterator<A> for Vec3<A> {
		fn from_iter<T: IntoIterator<Item = A>>(iter: T) -> Self {
			let mut iter = iter.into_iter();
			Vec3(iter.next().unwrap(),
				 iter.next().unwrap(),
				 iter.next().unwrap())
		}
	}
	
	pub struct Iter<T>(Vec3<T>, usize);
	
	impl<T: Clone> Iterator for Iter<T> {
		type Item = T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some((self.0).0.clone()),
				2 => Some((self.0).1.clone()),
				3 => Some((self.0).2.clone()),
				_ => None,
			}
		}
	}
	
	impl<T: Clone> IntoIterator for Vec3<T> {
		type Item = T;
		type IntoIter = Iter<T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			Iter(self, 0)
		}
	}
	
	pub struct IterRef<'a, T>(&'a Vec3<T>, usize);
	
	impl<'a, T> Iterator for IterRef<'a, T> {
		type Item = &'a T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(&(self.0).0),
				2 => Some(&(self.0).1),
				3 => Some(&(self.0).2),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a Vec3<T> {
		type Item = &'a T;
		type IntoIter = IterRef<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterRef(self, 0)
		}
	}
	
	pub struct IterMut<'a, T>(&'a mut Vec3<T>, usize);
	
	impl<'a, T> Iterator for IterMut<'a, T> {
		type Item = &'a mut T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(unsafe { ::std::mem::transmute(&mut (self.0).0) }),
				2 => Some(unsafe { ::std::mem::transmute(&mut (self.0).1) }),
				3 => Some(unsafe { ::std::mem::transmute(&mut (self.0).2) }),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a mut Vec3<T> {
		type Item = &'a mut T;
		type IntoIter = IterMut<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterMut(self, 0)
		}
	}
}

pub mod vec4 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Vec4<T>(pub T, pub T, pub T, pub T);
	
	impl<T> Vec4<T> {
		#[inline]
		pub fn map<F: FnMut(T) -> U, U>(self, mut f: F) -> Vec4<U> {
			Vec4(f(self.0),
				 f(self.1),
				 f(self.2),
				 f(self.3))
		}
	}
	
	impl<T: Default> Default for Vec4<T> {
		fn default() -> Self {
			Self(T::default(),
				 T::default(),
				 T::default(),
				 T::default())
		}
	}
	
	impl<T> Vec4<T> {
		#[inline]
		pub fn xyz(self) -> Vec3<T> {
			Vec3(self.0, self.1, self.2)
		}
	}
	
	impl<T: Copy> From<T> for Vec4<T> {
		fn from(v: T) -> Self {
			Self(v, v, v, v)
		}
	}
	
	impl<T: Copy> From<[T; 4]> for Vec4<T> {
		fn from(v: [T; 4]) -> Self {
			Self(v[0], v[1], v[2], v[3])
		}
	}
	
	impl<T> From<(T, T, T, T)> for Vec4<T> {
		fn from(v: (T, T, T, T)) -> Self {
			Self(v.0, v.1, v.2, v.3)
		}
	}
	
	impl Vec4<f32> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs(),
				 self.2.abs(),
				 self.3.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1),
				 self.2.powf(rhs.2),
				 self.3.powf(rhs.3))
		}
		
		pub fn powi(self, rhs: Vec4<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1),
				 self.2.powi(rhs.2),
				 self.3.powi(rhs.3))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt(),
				 self.2.sqrt(),
				 self.3.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt(),
				 self.2.cbrt(),
				 self.3.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1),
				 self.2.log(rhs.2),
				 self.3.log(rhs.3))
		}
		
		pub fn len2(self) -> f32 {
			self.dot(self)
		}
		
		pub fn len(self) -> f32 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f32) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance2(self, other: Self) -> f32 {
			(other - self).len2()
		}
		
		pub fn distance(self, other: Self) -> f32 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f32 {
			self.dot(other) / (self.len() * other.len())
		}
	}
	
	impl Vec4<f64> {
		pub fn abs(self) -> Self {
			Self(self.0.abs(),
				 self.1.abs(),
				 self.2.abs(),
				 self.3.abs())
		}
		
		pub fn powf(self, rhs: Self) -> Self {
			Self(self.0.powf(rhs.0),
				 self.1.powf(rhs.1),
				 self.2.powf(rhs.2),
				 self.3.powf(rhs.3))
		}
		
		pub fn powi(self, rhs: Vec4<i32>) -> Self {
			Self(self.0.powi(rhs.0),
				 self.1.powi(rhs.1),
				 self.2.powi(rhs.2),
				 self.3.powi(rhs.3))
		}
		
		pub fn sqrt(self) -> Self {
			Self(self.0.sqrt(),
				 self.1.sqrt(),
				 self.2.sqrt(),
				 self.3.sqrt())
		}
		
		pub fn cbrt(self) -> Self {
			Self(self.0.cbrt(),
				 self.1.cbrt(),
				 self.2.cbrt(),
				 self.3.cbrt())
		}
		
		pub fn log(self, rhs: Self) -> Self {
			Self(self.0.log(rhs.0),
				 self.1.log(rhs.1),
				 self.2.log(rhs.2),
				 self.2.log(rhs.3))
		}
		
		pub fn len2(self) -> f64 {
			self.dot(self)
		}
		
		pub fn len(self) -> f64 {
			self.dot(self).sqrt()
		}
		
		pub fn normalize(self, len: f64) -> Self {
			self * (len / self.len())
		}
		
		pub fn distance2(self, other: Self) -> f64 {
			(other - self).len2()
		}
		
		pub fn distance(self, other: Self) -> f64 {
			(other - self).len()
		}
		
		pub fn angle(self, other: Self) -> f64 {
			(self.dot(other) / (self.len() * other.len())).acos()
		}
	}
	
	impl Vec4<bool> {
		pub fn and(self) -> bool {
			self.0 && self.1 && self.2 && self.3
		}
		
		pub fn or(self) -> bool {
			self.0 || self.1 || self.2 || self.3
		}
		
		pub fn xor(self) -> bool {
			self.0 ^ self.1 ^ self.2 ^ self.3
		}
	}
	
	impl<T> Index<usize> for Vec4<T> {
		type Output = T;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				2 => &self.2,
				3 => &self.3,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Vec4<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				2 => &mut self.2,
				3 => &mut self.3,
				_ => panic!()
			}
		}
	}
	
	impl<T: PartialOrd + Copy> Vec4<T> {
		pub fn eq(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 == other.0,
				 self.1 == other.1,
				 self.2 == other.2,
				 self.3 == other.3)
		}
		
		pub fn eqs(self, other: T) -> Vec4<bool> {
			Vec4(self.0 == other,
				 self.1 == other,
				 self.2 == other,
				 self.3 == other)
		}
		
		
		pub fn ne(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 != other.0,
				 self.1 != other.1,
				 self.2 != other.2,
				 self.3 != other.3)
		}
		
		pub fn nes(self, other: T) -> Vec4<bool> {
			Vec4(self.0 != other,
				 self.1 != other,
				 self.2 != other,
				 self.3 != other)
		}
		
		pub fn gt(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 > other.0,
				 self.1 > other.1,
				 self.2 > other.2,
				 self.3 > other.3)
		}
		
		pub fn gts(self, other: T) -> Vec4<bool> {
			Vec4(self.0 > other,
				 self.1 > other,
				 self.2 > other,
				 self.3 > other)
		}
		
		
		pub fn lt(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 < other.0,
				 self.1 < other.1,
				 self.2 < other.2,
				 self.3 < other.3)
		}
		
		pub fn lts(self, other: T) -> Vec4<bool> {
			Vec4(self.0 < other,
				 self.1 < other,
				 self.2 < other,
				 self.3 < other)
		}
		
		pub fn ge(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 >= other.0,
				 self.1 >= other.1,
				 self.2 >= other.2,
				 self.3 >= other.3)
		}
		
		pub fn ges(self, other: T) -> Vec4<bool> {
			Vec4(self.0 >= other,
				 self.1 >= other,
				 self.2 >= other,
				 self.3 >= other)
		}
		
		pub fn le(self, other: Self) -> Vec4<bool> {
			Vec4(self.0 <= other.0,
				 self.1 <= other.1,
				 self.2 <= other.2,
				 self.3 <= other.3)
		}
		
		pub fn les(self, other: T) -> Vec4<bool> {
			Vec4(self.0 <= other,
				 self.1 <= other,
				 self.2 <= other,
				 self.3 <= other)
		}
		
		pub fn min(self, other: Self) -> Self {
			Self(if self.0 < other.0 { self.0 } else { other.0 },
				 if self.1 < other.1 { self.1 } else { other.1 },
				 if self.2 < other.2 { self.2 } else { other.2 },
				 if self.3 < other.3 { self.3 } else { other.3 })
		}
		
		pub fn mins(self, other: T) -> Self {
			Self(if self.0 < other { self.0 } else { other },
				 if self.1 < other { self.1 } else { other },
				 if self.2 < other { self.2 } else { other },
				 if self.3 < other { self.3 } else { other })
		}
		
		pub fn max(self, other: Self) -> Self {
			Self(if self.0 > other.0 { self.0 } else { other.0 },
				 if self.1 > other.1 { self.1 } else { other.1 },
				 if self.2 > other.2 { self.2 } else { other.2 },
				 if self.3 > other.3 { self.3 } else { other.3 })
		}
		
		pub fn maxs(self, other: T) -> Self {
			Self(if self.0 > other { self.0 } else { other },
				 if self.1 > other { self.1 } else { other },
				 if self.2 > other { self.2 } else { other },
				 if self.3 > other { self.3 } else { other })
		}
		
		pub fn clamp(self, lower: Self, upper: Self) -> Self {
			Self(if self.0 < lower.0 { lower.0 } else if self.0 > upper.0 { upper.0 } else { self.0 },
				 if self.1 < lower.1 { lower.1 } else if self.1 > upper.1 { upper.1 } else { self.1 },
				 if self.2 < lower.2 { lower.2 } else if self.2 > upper.2 { upper.2 } else { self.2 },
				 if self.3 < lower.3 { lower.3 } else if self.3 > upper.3 { upper.3 } else { self.3 })
		}
		
		pub fn clamps(self, lower: T, upper: T) -> Self {
			Self(if self.0 < lower { lower } else if self.0 > upper { upper } else { self.0 },
				 if self.1 < lower { lower } else if self.1 > upper { upper } else { self.1 },
				 if self.2 < lower { lower } else if self.2 > upper { upper } else { self.2 },
				 if self.3 < lower { lower } else if self.3 > upper { upper } else { self.3 })
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Copy> Vec4<T> {
		pub fn dot(self, other: Self) -> T {
			self.0 * other.0 + self.1 * other.1 + self.2 * other.2 + self.3 * other.3
		}
		
		pub fn mix(self, a: Self, b: Self) -> Self {
			Self(a.0 + (b.0 - a.0) * self.0,
				 a.1 + (b.1 - a.1) * self.1,
				 a.2 + (b.2 - a.2) * self.2,
				 a.3 + (b.3 - a.3) * self.3)
		}
	}
	
	impl<T: Add> Add for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn add(self, rhs: Self) -> Self::Output {
			Vec4(self.0 + rhs.0,
				 self.1 + rhs.1,
				 self.2 + rhs.2,
				 self.3 + rhs.3)
		}
	}
	
	impl<T: Sub> Sub for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn sub(self, rhs: Self) -> Self::Output {
			Vec4(self.0 - rhs.0,
				 self.1 - rhs.1,
				 self.2 - rhs.2,
				 self.3 - rhs.3)
		}
	}
	
	impl<T: Mul> Mul for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn mul(self, rhs: Self) -> Self::Output {
			Vec4(self.0 * rhs.0,
				 self.1 * rhs.1,
				 self.2 * rhs.2,
				 self.3 * rhs.3)
		}
	}
	
	impl<T: Mul + Copy> Mul<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn mul(self, rhs: T) -> Self::Output {
			Vec4(self.0 * rhs,
				 self.1 * rhs,
				 self.2 * rhs,
				 self.3 * rhs)
		}
	}
	
	impl<T: Div> Div for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn div(self, rhs: Self) -> Self::Output {
			Vec4(self.0 / rhs.0,
				 self.1 / rhs.1,
				 self.2 / rhs.2,
				 self.3 / rhs.3)
		}
	}
	
	impl<T: Div + Copy> Div<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn div(self, rhs: T) -> Self::Output {
			Vec4(self.0 / rhs,
				 self.1 / rhs,
				 self.2 / rhs,
				 self.3 / rhs)
		}
	}
	
	impl<T: Rem> Rem for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn rem(self, rhs: Self) -> Self::Output {
			Vec4(self.0 % rhs.0,
				 self.1 % rhs.1,
				 self.2 % rhs.2,
				 self.3 % rhs.3)
		}
	}
	
	impl<T: Rem + Copy> Rem<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn rem(self, rhs: T) -> Self::Output {
			Vec4(self.0 % rhs,
				 self.1 % rhs,
				 self.2 % rhs,
				 self.3 % rhs)
		}
	}
	
	impl<T: Neg> Neg for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn neg(self) -> Self::Output {
			Vec4(-self.0,
				 -self.1,
				 -self.2,
				 -self.3)
		}
	}
	
	impl<T: Not> Not for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn not(self) -> Self::Output {
			Vec4(!self.0,
				 !self.1,
				 !self.2,
				 !self.3)
		}
	}
	
	impl<T: BitAnd> BitAnd for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitand(self, rhs: Self) -> Self::Output {
			Vec4(self.0 & rhs.0,
				 self.1 & rhs.1,
				 self.2 & rhs.2,
				 self.3 & rhs.3)
		}
	}
	
	impl<T: BitAnd + Copy> BitAnd<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitand(self, rhs: T) -> Self::Output {
			Vec4(self.0 & rhs,
				 self.1 & rhs,
				 self.2 & rhs,
				 self.3 & rhs)
		}
	}
	
	impl<T: BitOr> BitOr for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitor(self, rhs: Self) -> Self::Output {
			Vec4(self.0 | rhs.0,
				 self.1 | rhs.1,
				 self.2 | rhs.2,
				 self.3 | rhs.3)
		}
	}
	
	impl<T: BitOr + Copy> BitOr<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitor(self, rhs: T) -> Self::Output {
			Vec4(self.0 | rhs,
				 self.1 | rhs,
				 self.2 | rhs,
				 self.3 | rhs)
		}
	}
	
	impl<T: BitXor> BitXor for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitxor(self, rhs: Self) -> Self::Output {
			Vec4(self.0 ^ rhs.0,
				 self.1 ^ rhs.1,
				 self.2 ^ rhs.2,
				 self.3 ^ rhs.3)
		}
	}
	
	impl<T: BitXor + Copy> BitXor<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn bitxor(self, rhs: T) -> Self::Output {
			Vec4(self.0 ^ rhs,
				 self.1 ^ rhs,
				 self.2 ^ rhs,
				 self.3 ^ rhs)
		}
	}
	
	impl<T: Shl> Shl for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn shl(self, rhs: Self) -> Self::Output {
			Vec4(self.0 << rhs.0,
				 self.1 << rhs.1,
				 self.2 << rhs.2,
				 self.3 << rhs.3)
		}
	}
	
	impl<T: Shl + Copy> Shl<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn shl(self, rhs: T) -> Self::Output {
			Vec4(self.0 << rhs,
				 self.1 << rhs,
				 self.2 << rhs,
				 self.3 << rhs)
		}
	}
	
	impl<T: Shr> Shr for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn shr(self, rhs: Self) -> Self::Output {
			Vec4(self.0 >> rhs.0,
				 self.1 >> rhs.1,
				 self.2 >> rhs.2,
				 self.3 >> rhs.3)
		}
	}
	
	impl<T: Shr + Copy> Shr<T> for Vec4<T> {
		type Output = Vec4<T::Output>;
		
		fn shr(self, rhs: T) -> Self::Output {
			Vec4(self.0 >> rhs,
				 self.1 >> rhs,
				 self.2 >> rhs,
				 self.3 >> rhs)
		}
	}
	
	impl<T: AddAssign> AddAssign for Vec4<T> {
		fn add_assign(&mut self, rhs: Self) {
			self.0 += rhs.0;
			self.1 += rhs.1;
			self.2 += rhs.2;
			self.3 += rhs.3;
		}
	}
	
	impl<T: SubAssign> SubAssign for Vec4<T> {
		fn sub_assign(&mut self, rhs: Self) {
			self.0 -= rhs.0;
			self.1 -= rhs.1;
			self.2 -= rhs.2;
			self.3 -= rhs.3;
		}
	}
	
	impl<T: MulAssign> MulAssign for Vec4<T> {
		fn mul_assign(&mut self, rhs: Self) {
			self.0 *= rhs.0;
			self.1 *= rhs.1;
			self.2 *= rhs.2;
			self.3 *= rhs.3;
		}
	}
	
	impl<T: MulAssign + Copy> MulAssign<T> for Vec4<T> {
		fn mul_assign(&mut self, rhs: T) {
			self.0 *= rhs;
			self.1 *= rhs;
			self.2 *= rhs;
			self.3 *= rhs;
		}
	}
	
	impl<T: DivAssign> DivAssign for Vec4<T> {
		fn div_assign(&mut self, rhs: Self) {
			self.0 /= rhs.0;
			self.1 /= rhs.1;
			self.2 /= rhs.2;
			self.3 /= rhs.3;
		}
	}
	
	impl<T: DivAssign + Copy> DivAssign<T> for Vec4<T> {
		fn div_assign(&mut self, rhs: T) {
			self.0 /= rhs;
			self.1 /= rhs;
			self.2 /= rhs;
			self.3 /= rhs;
		}
	}
	
	impl<T: RemAssign> RemAssign for Vec4<T> {
		fn rem_assign(&mut self, rhs: Self) {
			self.0 %= rhs.0;
			self.1 %= rhs.1;
			self.2 %= rhs.2;
			self.3 %= rhs.3;
		}
	}
	
	impl<T: RemAssign + Copy> RemAssign<T> for Vec4<T> {
		fn rem_assign(&mut self, rhs: T) {
			self.0 %= rhs;
			self.1 %= rhs;
			self.2 %= rhs;
			self.3 %= rhs;
		}
	}
	
	impl<T: BitAndAssign> BitAndAssign for Vec4<T> {
		fn bitand_assign(&mut self, rhs: Self) {
			self.0 &= rhs.0;
			self.1 &= rhs.1;
			self.2 &= rhs.2;
			self.3 &= rhs.3;
		}
	}
	
	impl<T: BitAndAssign + Copy> BitAndAssign<T> for Vec4<T> {
		fn bitand_assign(&mut self, rhs: T) {
			self.0 &= rhs;
			self.1 &= rhs;
			self.2 &= rhs;
			self.3 &= rhs;
		}
	}
	
	impl<T: BitOrAssign> BitOrAssign for Vec4<T> {
		fn bitor_assign(&mut self, rhs: Self) {
			self.0 |= rhs.0;
			self.1 |= rhs.1;
			self.2 |= rhs.2;
			self.3 |= rhs.3;
		}
	}
	
	impl<T: BitOrAssign + Copy> BitOrAssign<T> for Vec4<T> {
		fn bitor_assign(&mut self, rhs: T) {
			self.0 |= rhs;
			self.1 |= rhs;
			self.2 |= rhs;
			self.3 |= rhs;
		}
	}
	
	impl<T: BitXorAssign> BitXorAssign for Vec4<T> {
		fn bitxor_assign(&mut self, rhs: Self) {
			self.0 ^= rhs.0;
			self.1 ^= rhs.1;
			self.2 ^= rhs.2;
			self.3 ^= rhs.3;
		}
	}
	
	impl<T: BitXorAssign + Copy> BitXorAssign<T> for Vec4<T> {
		fn bitxor_assign(&mut self, rhs: T) {
			self.0 ^= rhs;
			self.1 ^= rhs;
			self.2 ^= rhs;
			self.3 ^= rhs;
		}
	}
	
	impl<T: ShlAssign> ShlAssign for Vec4<T> {
		fn shl_assign(&mut self, rhs: Self) {
			self.0 <<= rhs.0;
			self.1 <<= rhs.1;
			self.2 <<= rhs.2;
			self.3 <<= rhs.3;
		}
	}
	
	impl<T: ShlAssign + Copy> ShlAssign<T> for Vec4<T> {
		fn shl_assign(&mut self, rhs: T) {
			self.0 <<= rhs;
			self.1 <<= rhs;
			self.2 <<= rhs;
			self.3 <<= rhs;
		}
	}
	
	impl<T: ShrAssign> ShrAssign for Vec4<T> {
		fn shr_assign(&mut self, rhs: Self) {
			self.0 >>= rhs.0;
			self.1 >>= rhs.1;
			self.2 >>= rhs.2;
			self.3 >>= rhs.3;
		}
	}
	
	impl<T: ShrAssign + Copy> ShrAssign<T> for Vec4<T> {
		fn shr_assign(&mut self, rhs: T) {
			self.0 >>= rhs;
			self.1 >>= rhs;
			self.2 >>= rhs;
			self.3 >>= rhs;
		}
	}
	
	impl<A> FromIterator<A> for Vec4<A> {
		fn from_iter<T: IntoIterator<Item = A>>(iter: T) -> Self {
			let mut iter = iter.into_iter();
			Vec4(iter.next().unwrap(),
				 iter.next().unwrap(),
				 iter.next().unwrap(),
				 iter.next().unwrap())
		}
	}
	
	pub struct Iter<T>(Vec4<T>, usize);
	
	impl<T: Clone> Iterator for Iter<T> {
		type Item = T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some((self.0).0.clone()),
				2 => Some((self.0).1.clone()),
				3 => Some((self.0).2.clone()),
				4 => Some((self.0).3.clone()),
				_ => None,
			}
		}
	}
	
	impl<T: Clone> IntoIterator for Vec4<T> {
		type Item = T;
		type IntoIter = Iter<T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			Iter(self, 0)
		}
	}
	
	pub struct IterRef<'a, T>(&'a Vec4<T>, usize);
	
	impl<'a, T> Iterator for IterRef<'a, T> {
		type Item = &'a T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(&(self.0).0),
				2 => Some(&(self.0).1),
				3 => Some(&(self.0).2),
				4 => Some(&(self.0).3),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a Vec4<T> {
		type Item = &'a T;
		type IntoIter = IterRef<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterRef(self, 0)
		}
	}
	
	pub struct IterMut<'a, T>(&'a mut Vec4<T>, usize);
	
	impl<'a, T> Iterator for IterMut<'a, T> {
		type Item = &'a mut T;
		
		fn next(&mut self) -> Option<Self::Item> {
			self.1 += 1;
			match self.1 {
				1 => Some(unsafe { ::std::mem::transmute(&mut (self.0).0) }),
				2 => Some(unsafe { ::std::mem::transmute(&mut (self.0).1) }),
				3 => Some(unsafe { ::std::mem::transmute(&mut (self.0).2) }),
				4 => Some(unsafe { ::std::mem::transmute(&mut (self.0).3) }),
				_ => None,
			}
		}
	}
	
	impl<'a, T> IntoIterator for &'a mut Vec4<T> {
		type Item = &'a mut T;
		type IntoIter = IterMut<'a, T>;
		
		#[inline]
		fn into_iter(self) -> Self::IntoIter {
			IterMut(self, 0)
		}
	}
}

pub mod quat32 {
	use super::*;
	
	pub type Quat32 = Vec4<f32>;
	
	impl Quat32 {
		pub fn from_rotation_x(angle: f32) -> Self {
			Self((angle * 0.5f32).sin(), 0f32, 0f32, (angle * 0.5f32).cos())
		}
		
		pub fn from_rotation_y(angle: f32) -> Self {
			Self(0f32, (angle * 0.5f32).sin(), 0f32, (angle * 0.5f32).cos())
		}
		
		pub fn from_rotation_z(angle: f32) -> Self {
			Self(0f32, 0f32, (angle * 0.5f32).sin(), (angle * 0.5f32).cos())
		}
		
		pub fn from_axis_angle(axis: Vec3<f32>, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let inv_len = 1f32 / axis.len();
			Self(axis.0 * inv_len * sin,
				 axis.1 * inv_len * sin,
				 axis.2 * inv_len * sin,
				 (angle * 0.5f32).cos())
		}
		
		pub fn to_axis_angle(self) -> (Vec3<f32>, f32) {
			let acos = self.3.acos();
			let inv_sqrt = 1f32 / (1f32 - self.3 * self.3).sqrt();
			(Vec3(self.0 * inv_sqrt, self.1 * inv_sqrt, self.3 * inv_sqrt), acos + acos)
		}
		
		pub fn rotate_x_local(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.0 * cos + self.3 * sin,
				 self.1 * cos - self.2 * sin,
				 self.2 * cos + self.1 * sin,
				 self.3 * cos - self.0 * sin)
		}
		
		pub fn rotate_y_local(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.0 * cos + self.2 * sin,
				 self.1 * cos + self.3 * sin,
				 self.2 * cos - self.0 * sin,
				 self.3 * cos - self.1 * sin)
		}
		
		pub fn rotate_z_local(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.0 * cos - self.1 * sin,
				 self.1 * cos + self.0 * sin,
				 self.2 * cos + self.3 * sin,
				 self.3 * cos - self.2 * sin)
		}
		
		pub fn rotate_x(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.3 * sin  + self.0 * cos,
				 self.1 * cos  + self.2 * sin,
				 self.2 * cos  - self.1 * sin,
				 self.3 * cos  - self.0 * sin)
		}
		
		pub fn rotate_y(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.0 * cos  - self.2 * sin,
				 self.3 * sin  + self.1 * cos,
				 self.0 * sin  + self.2 * cos,
				 self.3 * cos  - self.1 * sin)
		}
		
		pub fn rotate_z(self, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let cos = (angle * 0.5f32).cos();
			Self(self.0 * cos  + self.1 * sin,
				 self.1 * cos  - self.0 * sin,
				 self.3 * sin  + self.2 * cos,
				 self.3 * cos  - self.2 * sin)
		}
		
		pub fn rotate_axis(self, axis: Vec3<f32>, angle: f32) -> Self {
			let sin = (angle * 0.5f32).sin();
			let inv_len = 1f32 / axis.len();
			let rx = axis.0 * inv_len * sin;
			let ry = axis.1 * inv_len * sin;
			let rz = axis.2 * inv_len * sin;
			let rw = (angle * 0.5f32).cos();
			Self(self.3 * rx + self.0 * rw + self.1 * rz - self.2 * ry,
				 self.3 * ry + self.0 * rz + self.1 * rw - self.2 * rx,
				 self.3 * rz + self.0 * ry + self.1 * rx - self.2 * rw,
				 self.3 * rw + self.0 * rx + self.1 * ry - self.2 * rz)
		}
	}
}

pub mod quat64 {
	use super::*;
	
	pub type Quat64 = Vec4<f64>;
	
	impl Quat64 {
		pub fn from_rotation_x(angle: f64) -> Self {
			Self((angle * 0.5f64).sin(), 0f64, 0f64, (angle * 0.5f64).cos())
		}
		
		pub fn from_rotation_y(angle: f64) -> Self {
			Self(0f64, (angle * 0.5f64).sin(), 0f64, (angle * 0.5f64).cos())
		}
		
		pub fn from_rotation_z(angle: f64) -> Self {
			Self(0f64, 0f64, (angle * 0.5f64).sin(), (angle * 0.5f64).cos())
		}
		
		pub fn from_axis_angle(axis: Vec3<f64>, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let inv_len = 1f64 / axis.len();
			Self(axis.0 * inv_len * sin,
				 axis.1 * inv_len * sin,
				 axis.2 * inv_len * sin,
				 (angle * 0.5f64).cos())
		}
		
		pub fn to_axis_angle(self) -> (Vec3<f64>, f64) {
			let acos = self.3.acos();
			let inv_sqrt = 1f64 / (1f64 - self.3 * self.3).sqrt();
			(Vec3(self.0 * inv_sqrt, self.1 * inv_sqrt, self.3 * inv_sqrt), acos + acos)
		}
		
		pub fn rotate_x_local(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.0 * cos + self.3 * sin,
				 self.1 * cos + self.2 * sin,
				 self.2 * cos - self.1 * sin,
				 self.3 * cos - self.0 * sin)
		}
		
		pub fn rotate_y_local(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.0 * cos + self.2 * sin,
				 self.1 * cos + self.3 * sin,
				 self.2 * cos - self.0 * sin,
				 self.3 * cos - self.1 * sin)
		}
		
		pub fn rotate_z_local(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.0 * cos + self.1 * sin,
				 self.1 * cos + self.0 * sin,
				 self.2 * cos - self.3 * sin,
				 self.3 * cos - self.2 * sin)
		}
		
		pub fn rotate_x(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.3 * sin  + self.0 * cos,
				 self.1 * cos  + self.2 * sin,
				 self.2 * cos  - self.1 * sin,
				 self.3 * cos  - self.0 * sin)
		}
		
		pub fn rotate_y(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.0 * cos  - self.2 * sin,
				 self.3 * sin  + self.1 * cos,
				 self.0 * sin  + self.2 * cos,
				 self.3 * cos  - self.1 * sin)
		}
		
		pub fn rotate_z(self, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let cos = (angle * 0.5f64).cos();
			Self(self.0 * cos  + self.1 * sin,
				 self.1 * cos  + self.0 * sin,
				 self.3 * sin  - self.2 * cos,
				 self.3 * cos  - self.2 * sin)
		}
		
		pub fn rotate_axis(self, axis: Vec3<f64>, angle: f64) -> Self {
			let sin = (angle * 0.5f64).sin();
			let inv_len = 1f64 / axis.len();
			let rx = axis.0 * inv_len * sin;
			let ry = axis.1 * inv_len * sin;
			let rz = axis.2 * inv_len * sin;
			let rw = (angle * 0.5f64).cos();
			Self(self.3 * rx + self.0 * rw + self.1 * rz - self.2 * ry,
				 self.3 * ry + self.0 * rz + self.1 * rw - self.2 * rx,
				 self.3 * rz + self.0 * ry + self.1 * rx - self.2 * rw,
				 self.3 * rw + self.0 * rx + self.1 * ry - self.2 * rz)
		}
	}
}

pub mod mat2 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Mat2<T>(pub Vec2<T>, pub Vec2<T>);
	
	impl Default for Mat2<f32> {
		fn default() -> Self {
			Self(Vec2(1f32, 0f32),
				 Vec2(0f32, 1f32))
		}
	}
	
	impl Default for Mat2<f64> {
		fn default() -> Self {
			Self(Vec2(1f64, 0f64),
				 Vec2(0f64, 1f64))
		}
	}
	
	impl<T: Copy> From<[T; 4]> for Mat2<T> {
		fn from(v: [T; 4]) -> Self {
			Self(Vec2(v[0], v[1]),
				 Vec2(v[2], v[3]))
		}
	}
	
	impl<T: Copy> From<[[T; 2]; 2]> for Mat2<T> {
		fn from(v: [[T; 2]; 2]) -> Self {
			Self(Vec2::from(v[0]),
				 Vec2::from(v[1]))
		}
	}
	
	impl<T> Index<usize> for Mat2<T> {
		type Output = Vec2<T>;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Mat2<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				_ => panic!()
			}
		}
	}
}

pub mod mat3 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Mat3<T>(pub Vec3<T>, pub Vec3<T>, pub Vec3<T>);
	
	impl Default for Mat3<f32> {
		fn default() -> Self {
			Self(Vec3(1f32, 0f32, 0f32),
				 Vec3(0f32, 1f32, 0f32),
				 Vec3(0f32, 0f32, 1f32))
		}
	}
	
	impl Default for Mat3<f64> {
		fn default() -> Self {
			Self(Vec3(1f64, 0f64, 0f64),
				 Vec3(0f64, 1f64, 0f64),
				 Vec3(0f64, 0f64, 1f64))
		}
	}
	
	impl<T: Copy> From<[T; 9]> for Mat3<T> {
		fn from(v: [T; 9]) -> Self {
			Self(Vec3(v[0], v[1], v[2]),
				 Vec3(v[3], v[4], v[5]),
				 Vec3(v[6], v[7], v[8]))
		}
	}
	
	impl<T: Copy> From<[[T; 3]; 3]> for Mat3<T> {
		fn from(v: [[T; 3]; 3]) -> Self {
			Self(Vec3::from(v[0]),
				 Vec3::from(v[1]),
				 Vec3::from(v[2]))
		}
	}
	
	impl<T> Index<usize> for Mat3<T> {
		type Output = Vec3<T>;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				2 => &self.2,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Mat3<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				2 => &mut self.2,
				_ => panic!()
			}
		}
	}
}

pub mod mat4 {
	use super::*;
	
	#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
	pub struct Mat4<T>(pub Vec4<T>, pub Vec4<T>, pub Vec4<T>, pub Vec4<T>);
	
	impl Default for Mat4<f32> {
		fn default() -> Self {
			Self(Vec4(1f32, 0f32, 0f32, 0f32),
				 Vec4(0f32, 1f32, 0f32, 0f32),
				 Vec4(0f32, 0f32, 1f32, 0f32),
				 Vec4(0f32, 0f32, 0f32, 1f32))
		}
	}
	
	impl Default for Mat4<f64> {
		fn default() -> Self {
			Self(Vec4(1f64, 0f64, 0f64, 0f64),
				 Vec4(0f64, 1f64, 0f64, 0f64),
				 Vec4(0f64, 0f64, 1f64, 0f64),
				 Vec4(0f64, 0f64, 0f64, 1f64))
		}
	}
	
	impl<T: Copy> From<[T; 16]> for Mat4<T> {
		fn from(v: [T; 16]) -> Self {
			Self(Vec4(v[0], v[1], v[2], v[3]),
				 Vec4(v[4], v[5], v[6], v[7]),
				 Vec4(v[8], v[9], v[10], v[11]),
				 Vec4(v[12], v[13], v[14], v[15]))
		}
	}
	
	impl<T: Copy> From<[[T; 4]; 4]> for Mat4<T> {
		fn from(v: [[T; 4]; 4]) -> Self {
			Self(Vec4::from(v[0]),
				 Vec4::from(v[1]),
				 Vec4::from(v[2]),
				 Vec4::from(v[3]))
		}
	}
	
	impl<T> Index<usize> for Mat4<T> {
		type Output = Vec4<T>;
		
		fn index(&self, index: usize) -> &Self::Output {
			match index {
				0 => &self.0,
				1 => &self.1,
				2 => &self.2,
				3 => &self.3,
				_ => panic!()
			}
		}
	}
	
	impl<T> IndexMut<usize> for Mat4<T> {
		fn index_mut(&mut self, index: usize) -> &mut Self::Output {
			match index {
				0 => &mut self.0,
				1 => &mut self.1,
				2 => &mut self.2,
				3 => &mut self.3,
				_ => panic!()
			}
		}
	}
	
	impl Mat4<f32> {
		pub fn from_translation(trans: Vec3<f32>) -> Self {
			let mut mat = Self::default();
			(mat.3).0 = trans.0;
			(mat.3).1 = trans.1;
			(mat.3).2 = trans.2;
			mat
		}
		
		pub fn from_rotation(quat: Quat32) -> Self {
			let mut mat = Self::default();
			let w2 = quat.3 * quat.3;
			let x2 = quat.0 * quat.0;
			let y2 = quat.1 * quat.1;
			let z2 = quat.2 * quat.2;
			let zw = quat.2 * quat.3;
			let xy = quat.0 * quat.1;
			let xz = quat.0 * quat.2;
			let yw = quat.1 * quat.3;
			let yz = quat.1 * quat.2;
			let xw = quat.0 * quat.3;
			(mat.0).0 = w2 + x2 - z2 - y2;
			(mat.0).1 = xy + zw + zw + xy;
			(mat.0).2 = xz - yw + xz - yw;
			(mat.1).0 = -zw + xy - zw + xy;
			(mat.1).1 = y2 - z2 + w2 - x2;
			(mat.1).2 = yz + yz + xw + xw;
			(mat.2).0 = yw + xz + xz + yw;
			(mat.2).1 = yz + yz - xw - xw;
			(mat.2).2 = z2 - y2 - x2 + w2;
			mat
		}
		
		pub fn from_rotation_axis_angle(axis: Vec3<f32>, angle: f32) -> Self {
			let axis = axis.normalize(1f32);
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f32 - cos;
			let xyc = axis.0 * axis.1 * c;
			let xzc = axis.0 * axis.2 * c;
			let yzc = axis.1 * axis.2 * c;
			let xs = axis.0 * sin;
			let ys = axis.1 * sin;
			let zs = axis.2 * sin;
			Self(Vec4(cos + axis.0 * axis.0 * c, xyc + zs, xzc - ys, 0f32),
				 Vec4(xyc - zs, cos + axis.1 * axis.1 * c, yzc + xs, 0f32),
				 Vec4(xzc + ys, yzc - xs, cos + axis.2 * axis.2 * c, 0f32),
				 Vec4(0f32, 0f32, 0f32, 1f32))
		}
		
		pub fn from_scale(scale: Vec3<f32>) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = scale.0;
			(mat.1).1 = scale.1;
			(mat.2).2 = scale.2;
			mat
		}
		
		pub fn from_transform(t: Vec3<f32>, q: Quat32, s: Vec3<f32>) -> Self {
			let mut mat = Self::default();
			let dqx = q.0 + q.0;
			let dqy = q.1 + q.1;
			let dqz = q.2 + q.2;
			let q00 = dqx * q.0;
			let q11 = dqy * q.1;
			let q22 = dqz * q.2;
			let q01 = dqx * q.1;
			let q02 = dqx * q.2;
			let q03 = dqx * q.3;
			let q12 = dqy * q.2;
			let q13 = dqy * q.3;
			let q23 = dqz * q.3;
			(mat.0).0 = s.0 - (q11 + q22) * s.0;
			(mat.0).1 = (q01 + q23) * s.0;
			(mat.0).2 = (q02 - q13) * s.0;
			(mat.0).3 = 0f32;
			(mat.1).0 = (q01 - q23) * s.1;
			(mat.1).1 = s.1 - (q22 + q00) * s.1;
			(mat.1).2 = (q12 + q03) * s.1;
			(mat.1).3 = 0f32;
			(mat.2).0 = (q02 + q13) * s.2;
			(mat.2).1 = (q12 - q03) * s.2;
			(mat.2).2 = s.2 - (q11 + q00) * s.2;
			(mat.2).3 = 0f32;
			(mat.3).0 = t.0;
			(mat.3).1 = t.1;
			(mat.3).2 = t.2;
			(mat.3).3 = 1f32;
			mat
		}
		
		pub fn from_perspective(y_fov: f32, aspect: f32, z_near: f32, z_far: f32, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			let h = (y_fov * 0.5f32).tan();
			(mat.0).0 = 1f32 / (h * aspect);
			(mat.1).1 = 1f32 / h;
			if z_far > 0f32 && (z_far == ::std::f32::INFINITY || z_far == ::std::f32::NEG_INFINITY) {
				(mat.2).2 = 1E-6f32 - 1f32;
				(mat.3).2 = (1E-6f32 - if z_zero_to_one { 1f32 } else { 2f32 }) * z_near;
			} else if z_near > 0f32 && (z_near == ::std::f32::INFINITY || z_near == ::std::f32::NEG_INFINITY) {
				(mat.2).2 = (if z_zero_to_one { 0f32 } else { 1f32 }) - 1E-6f32;
				(mat.3).2 = ((if z_zero_to_one { 1f32 } else { 2f32 }) - 1E-6f32) * z_far;
			} else {
				(mat.2).2 = (if z_zero_to_one { z_far } else { z_far + z_near }) / (z_near - z_far);
				(mat.3).2 = (if z_zero_to_one { z_far } else { z_far + z_far }) * z_near / (z_near - z_far);
			}
			(mat.2).3 = -1f32;
			(mat.3).3 = 0f32;
			mat
		}
		
		pub fn from_ortho(left: f32, right: f32, bottom: f32, top: f32, z_near: f32, z_far: f32, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = 2f32 / (right - left);
			(mat.1).1 = 2f32 / (top - bottom);
			(mat.2).2 = (if z_zero_to_one { 1f32 } else { 2f32 }) / (z_near - z_far);
			(mat.3).0 = (right + left) / (left - right);
			(mat.3).1 = (top + bottom) / (bottom - top);
			(mat.3).2 = (if z_zero_to_one { z_near } else { z_far + z_near }) / (z_near - z_far);
			mat
		}
		
		pub fn from_ortho_symmetric(width: f32, height: f32, z_near: f32, z_far: f32, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = 2f32 / width;
			(mat.1).1 = 2f32 / height;
			(mat.2).2 = (if z_zero_to_one { 1f32 } else { 2f32 }) / (z_near - z_far);
			(mat.3).2 = (if z_zero_to_one { z_near } else { z_far + z_near }) / (z_near - z_far);
			mat
		}
		
		pub fn rotate(self, axis: Vec3<f32>, angle: f32) -> Self {
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f32 - cos;
			let m00 = axis.0 * axis.0 * c + cos;
			let m01 = axis.0 * axis.1 * c + axis.2 * sin;
			let m02 = axis.0 * axis.2 * c - axis.1 * sin;
			let m10 = axis.0 * axis.1 * c - axis.2 * sin;
			let m11 = axis.1 * axis.1 * c + cos;
			let m12 = axis.1 * axis.2 * c + axis.0 * sin;
			let m20 = axis.0 * axis.2 * c + axis.1 * sin;
			let m21 = axis.1 * axis.2 * c - axis.0 * sin;
			let m22 = axis.2 * axis.2 * c + cos;
			
			Self(Vec4((self.0).0 * m00 + (self.1).0 * m01 + (self.2).0 * m02,
					  (self.0).1 * m00 + (self.1).1 * m01 + (self.2).1 * m02,
					  (self.0).2 * m00 + (self.1).2 * m01 + (self.2).2 * m02,
					  (self.0).3 * m00 + (self.1).3 * m01 + (self.2).3 * m02),
				 Vec4((self.0).0 * m10 + (self.1).0 * m11 + (self.2).0 * m12,
					  (self.0).1 * m10 + (self.1).1 * m11 + (self.2).1 * m12,
					  (self.0).2 * m10 + (self.1).2 * m11 + (self.2).2 * m12,
					  (self.0).3 * m10 + (self.1).3 * m11 + (self.2).3 * m12),
				 Vec4((self.0).0 * m20 + (self.1).0 * m21 + (self.2).0 * m22,
					  (self.0).1 * m20 + (self.1).1 * m21 + (self.2).1 * m22,
					  (self.0).2 * m20 + (self.1).2 * m21 + (self.2).2 * m22,
					  (self.0).3 * m20 + (self.1).3 * m21 + (self.2).3 * m22),
				 self.3)
		}
		
		pub fn rotate_local(self, axis: Vec3<f32>, angle: f32) -> Self {
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f32 - cos;
			let m00 = axis.0 * axis.0 * c + cos;
			let m01 = axis.0 * axis.1 * c + axis.2 * sin;
			let m02 = axis.0 * axis.2 * c - axis.1 * sin;
			let m10 = axis.0 * axis.1 * c - axis.2 * sin;
			let m11 = axis.1 * axis.1 * c + cos;
			let m12 = axis.1 * axis.2 * c + axis.0 * sin;
			let m20 = axis.0 * axis.2 * c + axis.1 * sin;
			let m21 = axis.1 * axis.2 * c - axis.0 * sin;
			let m22 = axis.2 * axis.2 * c + cos;
			
			Self(Vec4(m00 * (self.0).0 + m10 * (self.0).1 + m20 * (self.0).2,
					  m01 * (self.0).0 + m11 * (self.0).1 + m21 * (self.0).2,
					  m02 * (self.0).0 + m12 * (self.0).1 + m22 * (self.0).2,
					  (self.0).3),
				 Vec4(m00 * (self.1).0 + m10 * (self.1).1 + m20 * (self.1).2,
					  m01 * (self.1).0 + m11 * (self.1).1 + m21 * (self.1).2,
					  m02 * (self.1).0 + m12 * (self.1).1 + m22 * (self.1).2,
					  (self.1).3),
				 Vec4(m00 * (self.2).0 + m10 * (self.2).1 + m20 * (self.2).2,
					  m01 * (self.2).0 + m11 * (self.2).1 + m21 * (self.2).2,
					  m02 * (self.2).0 + m12 * (self.2).1 + m22 * (self.2).2,
					  (self.2).3),
				 Vec4(m00 * (self.3).0 + m10 * (self.3).1 + m20 * (self.3).2,
					  m01 * (self.3).0 + m11 * (self.3).1 + m21 * (self.3).2,
					  m02 * (self.3).0 + m12 * (self.3).1 + m22 * (self.3).2,
					  (self.3).3))
		}
		
		pub fn get_scale(&self) -> Vec3<f32> {
			Vec3(((self.0).0 * (self.0).0 + (self.0).1 * (self.0).1 + (self.0).2 * (self.0).2).sqrt(),
				 ((self.1).0 * (self.1).0 + (self.1).1 * (self.1).1 + (self.1).2 * (self.1).2).sqrt(),
				 ((self.2).0 * (self.2).0 + (self.2).1 * (self.2).1 + (self.2).2 * (self.2).2).sqrt())
		}
		
		pub fn get_rotation(&self) -> Quat32 {
			let tr = (self.0).0 + (self.1).1 + (self.2).2;
			if tr >= 0f32 {
				let t = (tr + 1f32).sqrt();
				let w = 0.5f32 * t;
				let t = 0.5f32 / t;
				Vec4(((self.1).2 - (self.2).1) * t,
					 ((self.2).0 - (self.0).2) * t,
					 ((self.0).1 - (self.1).0) * t,
					 w)
			} else if (self.0).0 >= (self.1).1 && (self.0).0 >= (self.2).2 {
				let t = ((self.0).0 - ((self.1).1 + (self.2).2) + 1f32).sqrt();
				let x = 0.5f32 * t;
				let t = 0.5f32 / t;
				Vec4(x,
					 ((self.1).0 + (self.0).1) * t,
					 ((self.0).2 + (self.2).0) * t,
					 ((self.1).2 - (self.2).1) * t)
			} else if (self.1).1 > (self.2).2 {
				let t = ((self.1).1 - ((self.2).2 + (self.0).0) + 1f32).sqrt();
				let y = 0.5f32 * t;
				let t = 0.5f32 / t;
				Vec4(((self.1).0 + (self.0).1) * t,
					 y,
					 ((self.2).1 + (self.1).2) * t,
					 ((self.2).0 - (self.0).2) * t)
			} else {
				let t = ((self.2).2 - ((self.0).0 + (self.1).1) + 1f32).sqrt();
				let z = 0.5f32 * t;
				let t = 0.5f32 / t;
				Vec4(((self.0).2 + (self.2).0) * t,
					 ((self.2).1 + (self.1).2) * t,
					 z,
					 ((self.0).1 - (self.1).0) * t)
			}
		}
		
		pub fn get_rotation_axis_angle(self) -> (Vec3<f32>, f32) {
			if ((self.1).0 - (self.0).1).abs() < 1E-4f32
				&& ((self.2).0 - (self.0).2).abs() < 1E-4f32
				&& ((self.2).1 - (self.1).2).abs() < 1E-4f32 {
				let xx = ((self.0).0 + 1f32) / 2f32;
				let yy = ((self.1).1 + 1f32) / 2f32;
				let zz = ((self.2).2 + 1f32) / 2f32;
				let xy = ((self.1).0 + (self.0).1) / 4f32;
				let xz = ((self.2).0 + (self.0).2) / 4f32;
				let yz = ((self.2).1 + (self.1).2) / 4f32;
				(if xx > yy && xx > zz {
					let x = xx.sqrt();
					Vec3(x, xy / x, xz / x)
				} else if yy > zz {
					let y = yy.sqrt();
					Vec3(xy / y, y, yz / y)
				} else {
					let z = zz.sqrt();
					Vec3(xz / z, yz / z, z)
				}, ::std::f32::consts::PI)
			} else {
				let s = (((self.1).2 - (self.2).1) * ((self.1).2 - (self.2).1)
					+ ((self.2).0 - (self.0).2) * ((self.2).0 - (self.0).2)
					+ ((self.0).1 - (self.1).0) * ((self.0).1 - (self.1).0)).sqrt();
				(Vec3(((self.1).2 - (self.2).1) / s,
					  ((self.1).2 - (self.2).1) / s,
					  ((self.0).1 - (self.1).0) / s),
				 (((self.0).0 + (self.1).1 + (self.2).2 - 1f32) / 2f32).acos())
			}
		}
		
		pub fn normalize3x3(self) -> Self {
			let inv_x_len = 1f32 / ((self.0).0 * (self.0).0 + (self.0).1 * (self.0).1 + (self.0).2 * (self.0).2).sqrt();
			let inv_y_len = 1f32 / ((self.1).0 * (self.1).0 + (self.1).1 * (self.1).1 + (self.1).2 * (self.1).2).sqrt();
			let inv_z_len = 1f32 / ((self.2).0 * (self.2).0 + (self.2).1 * (self.2).1 + (self.2).2 * (self.2).2).sqrt();
			let mut mat = Self(self.0 * inv_x_len,
							   self.1 * inv_y_len,
							   self.2 * inv_z_len,
							   self.3);
			(mat.0).3 = (self.0).3;
			(mat.1).3 = (self.1).3;
			(mat.2).3 = (self.2).3;
			mat
		}
	}
	
	impl Mat4<f64> {
		pub fn from_translation(trans: Vec3<f64>) -> Self {
			let mut mat = Self::default();
			(mat.3).0 = trans.0;
			(mat.3).1 = trans.1;
			(mat.3).2 = trans.2;
			mat
		}
		
		pub fn from_rotation(quat: Quat64) -> Self {
			let mut mat = Self::default();
			let w2 = quat.3 * quat.3;
			let x2 = quat.0 * quat.0;
			let y2 = quat.1 * quat.1;
			let z2 = quat.2 * quat.2;
			let zw = quat.2 * quat.3;
			let xy = quat.0 * quat.1;
			let xz = quat.0 * quat.2;
			let yw = quat.1 * quat.3;
			let yz = quat.1 * quat.2;
			let xw = quat.0 * quat.3;
			(mat.0).0 = w2 + x2 - z2 - y2;
			(mat.0).1 = xy + zw + zw + xy;
			(mat.0).2 = xz - yw + xz - yw;
			(mat.1).0 = -zw + xy - zw + xy;
			(mat.1).1 = y2 - z2 + w2 - x2;
			(mat.1).2 = yz + yz + xw + xw;
			(mat.2).0 = yw + xz + xz + yw;
			(mat.2).1 = yz + yz - xw - xw;
			(mat.2).2 = z2 - y2 - x2 + w2;
			mat
		}
		
		pub fn from_rotation_axis_angle(axis: Vec3<f64>, angle: f64) -> Self {
			let axis = axis.normalize(1f64);
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f64 - cos;
			let xyc = axis.0 * axis.1 * c;
			let xzc = axis.0 * axis.2 * c;
			let yzc = axis.1 * axis.2 * c;
			let xs = axis.0 * sin;
			let ys = axis.1 * sin;
			let zs = axis.2 * sin;
			Self(Vec4(cos + axis.0 * axis.0 * c, xyc + zs, xzc - ys, 0f64),
				 Vec4(xyc - zs, cos + axis.1 * axis.1 * c, yzc + xs, 0f64),
				 Vec4(xzc + ys, yzc - xs, cos + axis.2 * axis.2 * c, 0f64),
				 Vec4(0f64, 0f64, 0f64, 1f64))
		}
		
		pub fn from_scale(scale: Vec3<f64>) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = scale.0;
			(mat.1).1 = scale.1;
			(mat.2).2 = scale.2;
			mat
		}
		
		pub fn from_transform(t: Vec3<f64>, q: Quat64, s: Vec3<f64>) -> Self {
			let mut mat = Self::default();
			let dqx = q.0 + q.0;
			let dqy = q.1 + q.1;
			let dqz = q.2 + q.2;
			let q00 = dqx * q.0;
			let q11 = dqy * q.1;
			let q22 = dqz * q.2;
			let q01 = dqx * q.1;
			let q02 = dqx * q.2;
			let q03 = dqx * q.3;
			let q12 = dqy * q.2;
			let q13 = dqy * q.3;
			let q23 = dqz * q.3;
			(mat.0).0 = s.0 - (q11 + q22) * s.0;
			(mat.0).1 = (q01 + q23) * s.0;
			(mat.0).2 = (q02 - q13) * s.0;
			(mat.0).3 = 0f64;
			(mat.1).0 = (q01 - q23) * s.1;
			(mat.1).1 = s.1 - (q22 + q00) * s.1;
			(mat.1).2 = (q12 + q03) * s.1;
			(mat.1).3 = 0f64;
			(mat.2).0 = (q02 + q13) * s.2;
			(mat.2).1 = (q12 - q03) * s.2;
			(mat.2).2 = s.2 - (q11 + q00) * s.2;
			(mat.2).3 = 0f64;
			(mat.3).0 = t.0;
			(mat.3).1 = t.1;
			(mat.3).2 = t.2;
			(mat.3).3 = 1f64;
			mat
		}
		
		pub fn from_perspective(y_fov: f64, aspect: f64, z_near: f64, z_far: f64, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			let h = (y_fov * 0.5f64).tan();
			(mat.0).0 = 1f64 / (h * aspect);
			(mat.1).1 = 1f64 / h;
			if z_far > 0f64 && (z_far == ::std::f64::INFINITY || z_far == ::std::f64::NEG_INFINITY) {
				(mat.2).2 = 1E-6f64 - 1f64;
				(mat.3).2 = (1E-6f64 - if z_zero_to_one { 1f64 } else { 2f64 }) * z_near;
			} else if z_near > 0f64 && (z_near == ::std::f64::INFINITY || z_near == ::std::f64::NEG_INFINITY) {
				(mat.2).2 = (if z_zero_to_one { 0f64 } else { 1f64 }) - 1E-6f64;
				(mat.3).2 = ((if z_zero_to_one { 1f64 } else { 2f64 }) - 1E-6f64) * z_far;
			} else {
				(mat.2).2 = (if z_zero_to_one { z_far } else { z_far + z_near }) / (z_near - z_far);
				(mat.3).2 = (if z_zero_to_one { z_far } else { z_far + z_far }) * z_near / (z_near - z_far);
			}
			(mat.2).3 = -1f64;
			(mat.3).3 = 0f64;
			mat
		}
		
		pub fn from_ortho(left: f64, right: f64, bottom: f64, top: f64, z_near: f64, z_far: f64, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = 2f64 / (right - left);
			(mat.1).1 = 2f64 / (top - bottom);
			(mat.2).2 = (if z_zero_to_one { 1f64 } else { 2f64 }) / (z_near - z_far);
			(mat.3).0 = (right + left) / (left - right);
			(mat.3).1 = (top + bottom) / (bottom - top);
			(mat.3).2 = (if z_zero_to_one { z_near } else { z_far + z_near }) / (z_near - z_far);
			mat
		}
		
		pub fn from_ortho_symmetric(width: f64, height: f64, z_near: f64, z_far: f64, z_zero_to_one: bool) -> Self {
			let mut mat = Self::default();
			(mat.0).0 = 2f64 / width;
			(mat.1).1 = 2f64 / height;
			(mat.2).2 = (if z_zero_to_one { 1f64 } else { 2f64 }) / (z_near - z_far);
			(mat.3).2 = (if z_zero_to_one { z_near } else { z_far + z_near }) / (z_near - z_far);
			mat
		}
		
		pub fn rotate(self, axis: Vec3<f64>, angle: f64) -> Self {
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f64 - cos;
			let m00 = axis.0 * axis.0 * c + cos;
			let m01 = axis.0 * axis.1 * c + axis.2 * sin;
			let m02 = axis.0 * axis.2 * c - axis.1 * sin;
			let m10 = axis.0 * axis.1 * c - axis.2 * sin;
			let m11 = axis.1 * axis.1 * c + cos;
			let m12 = axis.1 * axis.2 * c + axis.0 * sin;
			let m20 = axis.0 * axis.2 * c + axis.1 * sin;
			let m21 = axis.1 * axis.2 * c - axis.0 * sin;
			let m22 = axis.2 * axis.2 * c + cos;
			
			Self(Vec4((self.0).0 * m00 + (self.1).0 * m01 + (self.2).0 * m02,
					  (self.0).1 * m00 + (self.1).1 * m01 + (self.2).1 * m02,
					  (self.0).2 * m00 + (self.1).2 * m01 + (self.2).2 * m02,
					  (self.0).3 * m00 + (self.1).3 * m01 + (self.2).3 * m02),
				 Vec4((self.0).0 * m10 + (self.1).0 * m11 + (self.2).0 * m12,
					  (self.0).1 * m10 + (self.1).1 * m11 + (self.2).1 * m12,
					  (self.0).2 * m10 + (self.1).2 * m11 + (self.2).2 * m12,
					  (self.0).3 * m10 + (self.1).3 * m11 + (self.2).3 * m12),
				 Vec4((self.0).0 * m20 + (self.1).0 * m21 + (self.2).0 * m22,
					  (self.0).1 * m20 + (self.1).1 * m21 + (self.2).1 * m22,
					  (self.0).2 * m20 + (self.1).2 * m21 + (self.2).2 * m22,
					  (self.0).3 * m20 + (self.1).3 * m21 + (self.2).3 * m22),
				 self.3)
		}
		
		pub fn rotate_local(self, axis: Vec3<f64>, angle: f64) -> Self {
			let sin = angle.sin();
			let cos = angle.cos();
			let c = 1f64 - cos;
			let m00 = axis.0 * axis.0 * c + cos;
			let m01 = axis.0 * axis.1 * c + axis.2 * sin;
			let m02 = axis.0 * axis.2 * c - axis.1 * sin;
			let m10 = axis.0 * axis.1 * c - axis.2 * sin;
			let m11 = axis.1 * axis.1 * c + cos;
			let m12 = axis.1 * axis.2 * c + axis.0 * sin;
			let m20 = axis.0 * axis.2 * c + axis.1 * sin;
			let m21 = axis.1 * axis.2 * c - axis.0 * sin;
			let m22 = axis.2 * axis.2 * c + cos;
			
			Self(Vec4(m00 * (self.0).0 + m10 * (self.0).1 + m20 * (self.0).2,
					  m01 * (self.0).0 + m11 * (self.0).1 + m21 * (self.0).2,
					  m02 * (self.0).0 + m12 * (self.0).1 + m22 * (self.0).2,
					  (self.0).3),
				 Vec4(m00 * (self.1).0 + m10 * (self.1).1 + m20 * (self.1).2,
					  m01 * (self.1).0 + m11 * (self.1).1 + m21 * (self.1).2,
					  m02 * (self.1).0 + m12 * (self.1).1 + m22 * (self.1).2,
					  (self.1).3),
				 Vec4(m00 * (self.2).0 + m10 * (self.2).1 + m20 * (self.2).2,
					  m01 * (self.2).0 + m11 * (self.2).1 + m21 * (self.2).2,
					  m02 * (self.2).0 + m12 * (self.2).1 + m22 * (self.2).2,
					  (self.2).3),
				 Vec4(m00 * (self.3).0 + m10 * (self.3).1 + m20 * (self.3).2,
					  m01 * (self.3).0 + m11 * (self.3).1 + m21 * (self.3).2,
					  m02 * (self.3).0 + m12 * (self.3).1 + m22 * (self.3).2,
					  (self.3).3))
		}
		
		pub fn get_scale(&self) -> Vec3<f64> {
			Vec3(((self.0).0 * (self.0).0 + (self.0).1 * (self.0).1 + (self.0).2 * (self.0).2).sqrt(),
				 ((self.1).0 * (self.1).0 + (self.1).1 * (self.1).1 + (self.1).2 * (self.1).2).sqrt(),
				 ((self.2).0 * (self.2).0 + (self.2).1 * (self.2).1 + (self.2).2 * (self.2).2).sqrt())
		}
		
		pub fn get_rotation(&self) -> Quat64 {
			let tr = (self.0).0 + (self.1).1 + (self.2).2;
			if tr >= 0f64 {
				let t = (tr + 1f64).sqrt();
				let w = 0.5f64 * t;
				let t = 0.5f64 / t;
				Vec4(((self.1).2 - (self.2).1) * t,
					 ((self.2).0 - (self.0).2) * t,
					 ((self.0).1 - (self.1).0) * t,
					 w)
			} else if (self.0).0 >= (self.1).1 && (self.0).0 >= (self.2).2 {
				let t = ((self.0).0 - ((self.1).1 + (self.2).2) + 1f64).sqrt();
				let x = 0.5f64 * t;
				let t = 0.5f64 / t;
				Vec4(x,
					 ((self.1).0 + (self.0).1) * t,
					 ((self.0).2 + (self.2).0) * t,
					 ((self.1).2 - (self.2).1) * t)
			} else if (self.1).1 > (self.2).2 {
				let t = ((self.1).1 - ((self.2).2 + (self.0).0) + 1f64).sqrt();
				let y = 0.5f64 * t;
				let t = 0.5f64 / t;
				Vec4(((self.1).0 + (self.0).1) * t,
					 y,
					 ((self.2).1 + (self.1).2) * t,
					 ((self.2).0 - (self.0).2) * t)
			} else {
				let t = ((self.2).2 - ((self.0).0 + (self.1).1) + 1f64).sqrt();
				let z = 0.5f64 * t;
				let t = 0.5f64 / t;
				Vec4(((self.0).2 + (self.2).0) * t,
					 ((self.2).1 + (self.1).2) * t,
					 z,
					 ((self.0).1 - (self.1).0) * t)
			}
		}
		
		pub fn get_rotation_axis_angle(self) -> (Vec3<f64>, f64) {
			if ((self.1).0 - (self.0).1).abs() < 1E-4f64
				&& ((self.2).0 - (self.0).2).abs() < 1E-4f64
				&& ((self.2).1 - (self.1).2).abs() < 1E-4f64 {
				let xx = ((self.0).0 + 1f64) / 2f64;
				let yy = ((self.1).1 + 1f64) / 2f64;
				let zz = ((self.2).2 + 1f64) / 2f64;
				let xy = ((self.1).0 + (self.0).1) / 4f64;
				let xz = ((self.2).0 + (self.0).2) / 4f64;
				let yz = ((self.2).1 + (self.1).2) / 4f64;
				(if xx > yy && xx > zz {
					let x = xx.sqrt();
					Vec3(x, xy / x, xz / x)
				} else if yy > zz {
					let y = yy.sqrt();
					Vec3(xy / y, y, yz / y)
				} else {
					let z = zz.sqrt();
					Vec3(xz / z, yz / z, z)
				}, ::std::f64::consts::PI)
			} else {
				let s = (((self.1).2 - (self.2).1) * ((self.1).2 - (self.2).1)
					+ ((self.2).0 - (self.0).2) * ((self.2).0 - (self.0).2)
					+ ((self.0).1 - (self.1).0) * ((self.0).1 - (self.1).0)).sqrt();
				(Vec3(((self.1).2 - (self.2).1) / s,
					  ((self.1).2 - (self.2).1) / s,
					  ((self.0).1 - (self.1).0) / s),
				 (((self.0).0 + (self.1).1 + (self.2).2 - 1f64) / 2f64).acos())
			}
		}
		
		pub fn normalize3x3(self) -> Self {
			let inv_x_len = 1f64 / ((self.0).0 * (self.0).0 + (self.0).1 * (self.0).1 + (self.0).2 * (self.0).2).sqrt();
			let inv_y_len = 1f64 / ((self.1).0 * (self.1).0 + (self.1).1 * (self.1).1 + (self.1).2 * (self.1).2).sqrt();
			let inv_z_len = 1f64 / ((self.2).0 * (self.2).0 + (self.2).1 * (self.2).1 + (self.2).2 * (self.2).2).sqrt();
			let mut mat = Self(self.0 * inv_x_len,
							   self.1 * inv_y_len,
							   self.2 * inv_z_len,
							   self.3);
			(mat.0).3 = (self.0).3;
			(mat.1).3 = (self.1).3;
			(mat.2).3 = (self.2).3;
			mat
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + MulAssign + DivAssign + Copy> Mat4<T> {
		pub fn mul_component_wise(self, rhs: Self) -> Self {
			Self(self.0 * rhs.0,
				 self.1 * rhs.1,
				 self.2 * rhs.2,
				 self.3 * rhs.3)
		}
		
		pub fn div_component_wise(self, rhs: Self) -> Self {
			Self(self.0 / rhs.0,
				 self.1 / rhs.1,
				 self.2 / rhs.2,
				 self.3 / rhs.3)
		}
		
		pub fn mul_assign_component_wise(&mut self, rhs: Self) {
			self.0 *= rhs.0;
			self.1 *= rhs.1;
			self.2 *= rhs.2;
			self.3 *= rhs.3;
		}
		
		pub fn div_assign_component_wise(&mut self, rhs: Self) {
			self.0 /= rhs.0;
			self.1 /= rhs.1;
			self.2 /= rhs.2;
			self.3 /= rhs.3;
		}
		
		pub fn det(&self) -> T {
			((self.0).0 * (self.1).1 - (self.0).1 * (self.1).0) * ((self.2).2 * (self.3).3 - (self.2).3 * (self.3).2)
				+ ((self.0).2 * (self.1).0 - (self.0).0 * (self.1).2) * ((self.2).1 * (self.3).3 - (self.2).3 * (self.3).1)
				+ ((self.0).0 * (self.1).3 - (self.0).3 * (self.1).0) * ((self.2).1 * (self.3).2 - (self.2).2 * (self.3).1)
				+ ((self.0).1 * (self.1).2 - (self.0).2 * (self.1).1) * ((self.2).0 * (self.3).3 - (self.2).3 * (self.3).0)
				+ ((self.0).3 * (self.1).1 - (self.0).1 * (self.1).3) * ((self.2).0 * (self.3).2 - (self.2).2 * (self.3).0)
				+ ((self.0).2 * (self.1).3 - (self.0).3 * (self.1).2) * ((self.2).0 * (self.3).1 - (self.2).1 * (self.3).0)
		}
		
		pub fn transform(&self, vec: Vec4<T>) -> Vec4<T> {
			Vec4((self.0).0 * vec.0 + (self.1).0 * vec.1 + (self.2).0 * vec.2 + (self.3).0 * vec.3,
				 (self.0).1 * vec.0 + (self.1).1 * vec.1 + (self.2).1 * vec.2 + (self.3).1 * vec.3,
				 (self.0).2 * vec.0 + (self.1).2 * vec.1 + (self.2).2 * vec.2 + (self.3).2 * vec.3,
				 (self.0).3 * vec.0 + (self.1).3 * vec.1 + (self.2).3 * vec.2 + (self.3).3 * vec.3)
		}
		
		pub fn transform_pos(&self, vec: Vec3<T>) -> Vec3<T> {
			Vec3((self.0).0 * vec.0 + (self.1).0 * vec.1 + (self.2).0 * vec.2 + (self.3).0,
				 (self.0).1 * vec.0 + (self.1).1 * vec.1 + (self.2).1 * vec.2 + (self.3).1,
				 (self.0).2 * vec.0 + (self.1).2 * vec.1 + (self.2).2 * vec.2 + (self.3).2)
		}
		
		pub fn transform_dir(&self, vec: Vec3<T>) -> Vec3<T> {
			Vec3((self.0).0 * vec.0 + (self.1).0 * vec.1 + (self.2).0 * vec.2,
				 (self.0).1 * vec.0 + (self.1).1 * vec.1 + (self.2).1 * vec.2,
				 (self.0).2 * vec.0 + (self.1).2 * vec.1 + (self.2).2 * vec.2)
		}
		
		pub fn scale(self, scale: Vec3<T>) -> Self {
			Self(self.0 * scale.0,
				 self.1 * scale.1,
				 self.2 * scale.2,
				 self.3)
		}
		
		pub fn scale_local(self, scale: Vec3<T>) -> Self {
			Self(Vec4(scale.0 * (self.0).0, scale.1 * (self.0).1, scale.2 * (self.0).2, (self.0).3),
				 Vec4(scale.0 * (self.1).0, scale.1 * (self.1).1, scale.2 * (self.1).2, (self.1).3),
				 Vec4(scale.0 * (self.2).0, scale.1 * (self.2).1, scale.2 * (self.2).2, (self.2).3),
				 Vec4(scale.0 * (self.3).0, scale.1 * (self.3).1, scale.2 * (self.3).2, (self.3).3))
		}
		
		pub fn translate(self, t: Vec3<T>) -> Self {
			Self(self.0,
				 self.1,
				 self.3,
				 Vec4((self.0).0 * t.0 + (self.1).0 * t.1 + (self.2).0 * t.2 + (self.3).0,
					  (self.0).1 * t.0 + (self.1).1 * t.1 + (self.2).1 * t.2 + (self.3).1,
					  (self.0).2 * t.0 + (self.1).2 * t.1 + (self.2).2 * t.2 + (self.3).2,
					  (self.0).3 * t.0 + (self.1).3 * t.1 + (self.2).3 * t.2 + (self.3).3))
		}
		
		pub fn translate_local(self, t: Vec3<T>) -> Self {
			Self(Vec4((self.0).0 + t.0 * (self.0).3,
					  (self.0).1 + t.1 * (self.0).3,
					  (self.0).2 + t.2 * (self.0).3,
					  (self.0).3),
				 Vec4((self.1).0 + t.0 * (self.1).3,
					  (self.1).1 + t.1 * (self.1).3,
					  (self.1).2 + t.2 * (self.1).3,
					  (self.1).3),
				 Vec4((self.2).0 + t.0 * (self.2).3,
					  (self.2).1 + t.1 * (self.2).3,
					  (self.2).2 + t.2 * (self.2).3,
					  (self.2).3),
				 Vec4((self.3).0 + t.0 * (self.3).3,
					  (self.3).1 + t.1 * (self.3).3,
					  (self.3).2 + t.2 * (self.3).3,
					  (self.3).3))
		}
		
		pub fn get_translation(&self) -> Vec3<T> {
			Vec3((self.0).0, (self.1).0, (self.2).0)
		}
	}
	
	impl<T: Add<Output = U>, U> Add for Mat4<T> {
		type Output = Mat4<U>;
		
		fn add(self, rhs: Self) -> Self::Output {
			Mat4(self.0 + rhs.0,
				 self.1 + rhs.1,
				 self.2 + rhs.2,
				 self.3 + rhs.3)
		}
	}
	
	impl<T: Sub<Output = U>, U> Sub for Mat4<T> {
		type Output = Mat4<U>;
		
		fn sub(self, rhs: Self) -> Self::Output {
			Mat4(self.0 - rhs.0,
				 self.1 - rhs.1,
				 self.2 - rhs.2,
				 self.3 - rhs.3)
		}
	}
	
	impl<T: Add<Output = T> + Mul<Output = T> + Copy> Mul for Mat4<T> {
		type Output = Self;
		
		fn mul(self, rhs: Self) -> Self::Output {
			Mat4(Vec4((self.0).0 * (rhs.0).0 + (self.1).0 * (rhs.0).1 + (self.2).0 * (rhs.0).2 + (self.3).0 * (rhs.0).3,
					  (self.0).1 * (rhs.0).0 + (self.1).1 * (rhs.0).1 + (self.2).1 * (rhs.0).2 + (self.3).1 * (rhs.0).3,
					  (self.0).2 * (rhs.0).0 + (self.1).2 * (rhs.0).1 + (self.2).2 * (rhs.0).2 + (self.3).2 * (rhs.0).3,
					  (self.0).3 * (rhs.0).0 + (self.1).3 * (rhs.0).1 + (self.2).3 * (rhs.0).2 + (self.3).3 * (rhs.0).3),
				 Vec4((self.0).0 * (rhs.1).0 + (self.1).0 * (rhs.1).1 + (self.2).0 * (rhs.1).2 + (self.3).0 * (rhs.1).3,
					  (self.0).1 * (rhs.1).0 + (self.1).1 * (rhs.1).1 + (self.2).1 * (rhs.1).2 + (self.3).1 * (rhs.1).3,
					  (self.0).2 * (rhs.1).0 + (self.1).2 * (rhs.1).1 + (self.2).2 * (rhs.1).2 + (self.3).2 * (rhs.1).3,
					  (self.0).3 * (rhs.1).0 + (self.1).3 * (rhs.1).1 + (self.2).3 * (rhs.1).2 + (self.3).3 * (rhs.1).3),
				 Vec4((self.0).0 * (rhs.2).0 + (self.1).0 * (rhs.2).1 + (self.2).0 * (rhs.2).2 + (self.3).0 * (rhs.2).3,
					  (self.0).1 * (rhs.2).0 + (self.1).1 * (rhs.2).1 + (self.2).1 * (rhs.2).2 + (self.3).1 * (rhs.2).3,
					  (self.0).2 * (rhs.2).0 + (self.1).2 * (rhs.2).1 + (self.2).2 * (rhs.2).2 + (self.3).2 * (rhs.2).3,
					  (self.0).3 * (rhs.2).0 + (self.1).3 * (rhs.2).1 + (self.2).3 * (rhs.2).2 + (self.3).3 * (rhs.2).3),
				 Vec4((self.0).0 * (rhs.3).0 + (self.1).0 * (rhs.3).1 + (self.2).0 * (rhs.3).2 + (self.3).0 * (rhs.3).3,
					  (self.0).1 * (rhs.3).0 + (self.1).1 * (rhs.3).1 + (self.2).1 * (rhs.3).2 + (self.3).1 * (rhs.3).3,
					  (self.0).2 * (rhs.3).0 + (self.1).2 * (rhs.3).1 + (self.2).2 * (rhs.3).2 + (self.3).2 * (rhs.3).3,
					  (self.0).3 * (rhs.3).0 + (self.1).3 * (rhs.3).1 + (self.2).3 * (rhs.3).2 + (self.3).3 * (rhs.3).3))
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + Neg<Output = T> + Copy> Div for Mat4<T> {
		type Output = Self;
		
		fn div(self, rhs: Self) -> Self::Output {
			self * (!rhs)
		}
	}
	
	impl<T: Neg<Output = U>, U> Neg for Mat4<T> {
		type Output = Mat4<U>;
		
		fn neg(self) -> Self::Output {
			Mat4(-self.0, -self.1, -self.2, -self.3)
		}
	}
	
	impl<T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + Neg<Output = T> + Copy> Not for Mat4<T> {
		type Output = Self;
		
		fn not(self) -> Self {
			let a = (self.0).0 * (self.1).1 - (self.0).1 * (self.1).0;
			let b = (self.0).0 * (self.1).2 - (self.0).2 * (self.1).0;
			let c = (self.0).0 * (self.1).3 - (self.0).3 * (self.1).0;
			let d = (self.0).1 * (self.1).2 - (self.0).2 * (self.1).1;
			let e = (self.0).1 * (self.1).3 - (self.0).3 * (self.1).1;
			let f = (self.0).2 * (self.1).3 - (self.0).3 * (self.1).2;
			let g = (self.2).0 * (self.3).1 - (self.2).1 * (self.3).0;
			let h = (self.2).0 * (self.3).2 - (self.2).2 * (self.3).0;
			let i = (self.2).0 * (self.3).3 - (self.2).3 * (self.3).0;
			let j = (self.2).1 * (self.3).2 - (self.2).2 * (self.3).1;
			let k = (self.2).1 * (self.3).3 - (self.2).3 * (self.3).1;
			let l = (self.2).2 * (self.3).3 - (self.2).3 * (self.3).2;
			let det = a * l - b * k + c * j + d * i - e * h + f * g;
			Self(Vec4(( (self.1).1 * l - (self.1).2 * k + (self.1).3 * j) / det,
					  (-(self.0).1 * l + (self.0).2 * k - (self.0).3 * j) / det,
					  ( (self.3).1 * f - (self.3).2 * e + (self.3).3 * d) / det,
					  (-(self.2).1 * f + (self.2).2 * e - (self.2).3 * d) / det),
				 Vec4((-(self.1).0 * l + (self.1).2 * i - (self.1).3 * h) / det,
					  ( (self.0).0 * l - (self.0).2 * i + (self.0).3 * h) / det,
					  (-(self.3).0 * f + (self.3).2 * c - (self.3).3 * b) / det,
					  ( (self.2).0 * f - (self.2).2 * c + (self.2).3 * b) / det),
				 Vec4(( (self.1).0 * k - (self.1).1 * i + (self.1).3 * g) / det,
					  (-(self.0).0 * k + (self.0).1 * i - (self.0).3 * g) / det,
					  ( (self.3).0 * e - (self.3).1 * c + (self.3).3 * a) / det,
					  (-(self.2).0 * e + (self.2).1 * c - (self.2).3 * a) / det),
				 Vec4((-(self.1).0 * j + (self.1).1 * h - (self.1).2 * g) / det,
					  ( (self.0).0 * j - (self.0).1 * h + (self.0).2 * g) / det,
					  (-(self.3).0 * d + (self.3).1 * b - (self.3).2 * a) / det,
					  ( (self.2).0 * d - (self.2).1 * b + (self.2).2 * a) / det))
		}
	}
	
	impl<T: AddAssign> AddAssign for Mat4<T> {
		fn add_assign(&mut self, rhs: Self) {
			self.0 += rhs.0;
			self.1 += rhs.1;
			self.2 += rhs.2;
			self.3 += rhs.3;
		}
	}
	
	impl<T: SubAssign> SubAssign for Mat4<T> {
		fn sub_assign(&mut self, rhs: Self) {
			self.0 -= rhs.0;
			self.1 -= rhs.1;
			self.2 -= rhs.2;
			self.3 -= rhs.3;
		}
	}
}

pub mod bezier1 {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Bezier1(pub Vec2<f32>, pub Vec2<f32>);
	
	impl Bezier1 {
		/// https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
		pub fn intersections(self, l0: Vec2<f32>, l1: Vec2<f32>) -> u32 {
			let Bezier1(p0, p1) = self;
			let d0 = p1 - p0;
			let d1 = l1 - l0;
			let t0 = ((p0.0 - l0.0) * d1.1 - (p0.1 - l0.1) * d1.0)
				* ((p1.0 - l0.0) * d1.1 - (p1.1 - l0.1) * d1.0);
			let t1 = ((l0.0 - p0.0) * d0.1 - (l0.1 - p0.1) * d0.0)
				* ((l1.0 - p0.0) * d0.1 - (l1.1 - p0.1) * d0.0);
			match t0 >= 0.0 && t0 <= 1.0 && t1 >= 0.0 && t1 <= 1.0 {
				true  => 0,
				false => 1
			}
		}
		
		pub fn distance2(self, p: Vec2<f32>) -> f32 {
			let Bezier1(p0, p1) = self;
			let a = p - p0;
			let b = p1 - p0;
			let e = a.dot(b);
			a.len2() - e * e / b.len2()
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn intersections() {
			let b = Bezier1(Vec2(0.0, 0.0), Vec2(1.0, 0.0));
			assert_eq!(0, b.intersections(Vec2(0.0, 1.0), Vec2(1.0, 1.0)));
			assert_eq!(1, b.intersections(Vec2(0.0, -0.5), Vec2(1.0, 0.5)));
		}
		
		#[test]
		fn distance() {
			let b = Bezier1(Vec2(0.0, 0.0), Vec2(1.0, 0.0));
			assert_eq!(1.0, b.distance2(Vec2(0.5, 1.0)));
		}
	}
}

pub mod bezier2 {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Bezier2(pub Vec2<f32>, pub Vec2<f32>, pub Vec2<f32>);
	
	impl Bezier2 {
		pub fn b(self, t: f32) -> Vec2<f32> {
			debug_assert!(t >= 0.0 && t <= 1.0);
			let Bezier2(p0, p1, p2) = self;
			let tinv = 1.0 - t;
			p0 * (tinv * tinv)
				+ p1 * 2.0 * t * tinv
				+ p2 * (t * t)
		}
		
		/// https://github.com/w8r/bezier-intersect/blob/master/src/quadratic.js
		pub fn intersections(self, p: Vec2<f32>, r: Vec2<f32>) -> usize {
			let Bezier2(p0, p1, p2) = self;
			
			/*
			// calculate coefficients
			let c0 = r.1 * (p0.0 - 2.0 * p1.0 + p2.0) - p.1 * (p0.0 - 2.0 * p1.0 + p2.0)
						+ p.0 * (p0.1 - 2.0 * p1.1 + p2.1) - r.0 * (p0.1 - 2.0 * p1.1 + p2.1);
			let c1 = r.1 * (2.0 * p1.0 - 2.0 * p0.0)  - p.1 * (2.0 * p1.0 - 2.0 * p0.0)
						+ p.0 * (2.0 * p1.1 - 2.0 * p0.1)  - r.0 * (2.0 * p1.1 - 2.0 * p0.1);
			let c2 = r.1 * p0.0 - p.1 * p0.0
				+ p.0 * p0.1 + r.0 * p0.1
				+ p.0 * (p.1 - r.1)
				+ p.1 * (r.0 - p.0);
			*/
			
			// calculate coefficients
			let c0 = p0 + p1 * -2.0 + p2;
			let c1 = p0 * -2.0 + p1 * 2.0;
			let c2 = p0;
			
			// Convert line to normal form: ax + by + c = 0
			let n = Vec2(p.1 - r.1, r.0 - p.0); // Find normal to line: negative inverse of original line's slope
			let cl = p.0 * r.1 - r.0 * p.1;      // Determine new c coefficient
			
			let min = p.0.min(r.0);
			let max = p.0.max(r.0);
			
			// solve quadratic equation
			quadratic_roots_iter(
				n.dot(c0),
				n.dot(c1),
				n.dot(c2) + cl)
				.filter(|t| *t >= 0.0 && *t <= 1.0) // check if point is on bezier
				.map(|t| self.b(t).0)                       // calculate point on bezier
				.filter(|b| *b >= min && *b <= max) // check if point is on line segment
				.count()                                           // count intersections
		}
		
		/// http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
		/// alternative: https://www.academia.edu/34616746/An_algorithm_for_computing_the_shortest_distance_between_a_point_and_quadratic_Bezier_curve
		pub fn distance2(self, p: Vec2<f32>) -> f32 {
			let Bezier2(p0, p1, p2) = self;
			
			let a = p1 - p0;
			let b = p2 - p1 - a;
			let c = p0 - p;
			
			cubic_roots_iter(
				b.dot(b),
				3.0 * a.dot(b),
				2.0 * a.dot(a) + c.dot(b),
				c.dot(a))
				.filter(|t| *t >= 0.0 && *t <= 1.0)
				.map(|t| p.distance2(self.b(t)))
				.fold(p.distance2(p0), f32::min)
				.min(p.distance2(p2))
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn intersections() {
			let b = Bezier2(Vec2(0.0, 0.0), Vec2(1.0, 0.0), Vec2(2.0, 0.0));
			assert_eq!(0, b.intersections(Vec2(0.0, 1.0), Vec2(1.0, 1.0)));
			assert_eq!(1, b.intersections(Vec2(0.0, -0.5), Vec2(1.0, 0.5)));
		}
		
		#[test]
		fn distance() {
			let b = Bezier2(Vec2(0.0, 0.0), Vec2(1.0, 0.0), Vec2(2.0, 0.0));
			assert_eq!(1.0, b.distance2(Vec2(0.5, 1.0)));
		}
	}
}

pub mod bezier3 {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Bezier3(pub Vec2<f32>, pub Vec2<f32>, pub Vec2<f32>, pub Vec2<f32>);
	
	/// https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Properties
	impl From<Bezier2> for Bezier3 {
		fn from(Bezier2(p0, p1, p2): Bezier2) -> Self {
			const FRAC_2_3: f32 = 2.0 / 3.0;
			Self(p0,
				 p0 * (1f32 - FRAC_2_3) + p1 * FRAC_2_3,
				 p2 * (1f32 - FRAC_2_3) + p1 * FRAC_2_3,
				 p2)
		}
	}
	
	impl Bezier3 {
		pub fn b(self, t: f32) -> Vec2<f32> {
			debug_assert!(t >= 0.0 && t <= 1.0);
			let Bezier3(p0, p1, p2, p3) = self;
			let t2 = t * t;
			let t3 = t * t2;
			let tinv = 1.0 - t;
			let tinv2 = tinv * tinv;
			let tinv3 = tinv * tinv2;
			p0 * tinv3
				+ p1 * 3.0 * tinv2 * t
				+ p2 * 3.0 * tinv2 * t2
				+ p3 * t3
		}
		
		/// https://github.com/w8r/bezier-intersect/blob/master/src/cubic.js
		pub fn intersections(self, p: Vec2<f32>, r: Vec2<f32>) -> usize {
			let Bezier3(p0, p1, p2, p3) = self;
			
			let min = p.min(r);
			let max = p.max(r);
			
			// calculate coefficients
			let c0 = p0 * -1.0 + p1 *  3.0 + p2 * -3.0 + p3;
			let c1 = p0 *  3.0 + p1 * -6.0 + p2 *  3.0;
			let c2 = p0 * -3.0 + p1 *  3.0;
			let c3 = p0;
			
			// Convert line to normal form: ax + by + c = 0
			let n = Vec2(p.1 - r.1, r.0 - p.0); // Find normal to line: negative inverse of original line's slope
			let cl = p.0 * r.1 - r.0 * p.1;      // Determine new c coefficient
			
			// ?Rotate each cubic coefficient using line for new coordinate system?
			// Find roots of rotated cubic
			cubic_roots_iter(
				n.dot(c0),
				n.dot(c1),
				n.dot(c2),
				n.dot(c3) + cl)
				.filter(|t| *t >= 0.0 && *t <= 1.0)            // check if t in [0;1] i.e. intersects Bezier
				.map(|t| self.b(t))                                    // calculate point on Bezier
				.filter(|b| b.0 >= min.0 && b.1 >= min.1
					&& b.0 <= max.0 && b.1 <= max.1)                          // check if point is on line segment
				.count()                                                      // count intersections
		}
		
		/// alternative: https://hal.inria.fr/file/index/docid/518379/filename/Xiao-DiaoChen2007c.pdf
		pub fn distance2(self, p: Vec2<f32>) -> f32 {
			const STEPS: usize = 1000;
			let mut distance = self.3.distance2(p);
			for i in 0..STEPS {
				distance = distance.min(self.b(1.0 / STEPS as f32 * i as f32).distance2(p));
			}
			distance
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn intersections() {
			let b = Bezier3(Vec2(0.0, 0.0), Vec2(1.0, 0.0), Vec2(2.0, 0.0), Vec2(3.0, 0.0));
			assert_eq!(0, b.intersections(Vec2(0.0, 1.0), Vec2(1.0, 1.0)));
			assert_eq!(1, b.intersections(Vec2(0.0, -0.5), Vec2(1.0, 0.5)));
		}
		
		#[test]
		fn distance() {
			let b = Bezier3(Vec2(0.0, 0.0), Vec2(1.0, 0.0), Vec2(2.0, 0.0), Vec2(3.0, 0.0));
			assert!((b.distance2(Vec2(0.5, 1.0)) - 1.0).abs() < 1e-4);
		}
	}
}

pub mod sphere {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Sphere {
		pub radius: f32
	}
	
	impl Distance<Vec3<f32>> for Sphere {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			(*point).len() - self.radius
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			assert_eq!(1f32, Sphere { radius: 0.5 }.distance(&Vec3(0.0, 0.0, 1.5)))
		}
	}
}

pub mod cuboid {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Cuboid(Vec3<f32>);
	
	impl Distance<Vec3<f32>> for Cuboid {
		fn distance(&self, p: &Vec3<f32>) -> f32 {
			let d = (*p).abs() - self.0;
			d.maxs(0f32).len() + d.0.max(d.1).max(d.2).min(0f32)
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			assert_eq!(1f32, Cuboid(Vec3(0.5, 0.5, 0.5)).distance(&Vec3(0.0, 0.0, 1.5)))
		}
	}
}

pub mod plane {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Plane(Vec4<f32>);
	
	impl Distance<Vec3<f32>> for Plane {
		fn distance(&self, p: &Vec3<f32>) -> f32 {
			(*p).dot(self.0.xyz()) + (self.0).3
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			// TODO add test
		}
	}
}

pub mod cylinder {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Cylinder(Vec2<f32>);
	
	impl Distance<Vec3<f32>> for Cylinder {
		fn distance(&self, p: &Vec3<f32>) -> f32 {
			let d = Vec2((*p).xz().len(), p.1).abs() - self.0;
			d.0.max(d.1).min(0.0) + d.maxs(0.0).len()
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			assert_eq!(1f32, Cylinder(Vec2(0.5, 0.5)).distance(&Vec3(0.0, 0.0, 1.5)))
		}
	}
}

pub mod triangle {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Triangle {
		pub a: Vec3<f32>,
		pub b: Vec3<f32>,
		pub c: Vec3<f32>
	}
	
	impl Distance2<Vec3<f32>> for Triangle {
		fn distance2(&self, p: &Vec3<f32>) -> f32 {
			let Triangle { a, b, c } = *self;
			let p = *p;
			let ba = b - a;
			let pa = p - a;
			let cb = c - b;
			let pb = p - b;
			let ac = a - c;
			let pc = p - c;
			let nor = ba.cross(ac);
			
			if ba.cross(nor).dot(pa).sign() +
				cb.cross(nor).dot(pb).sign() +
				ac.cross(nor).dot(pc).sign() < 2 {
				(ba * (ba.dot(pa) / ba.dot(ba)).clamp(0f32, 1f32) - pa).len2()
					.min((cb * (cb.dot(pb) / cb.dot(cb)).clamp(0f32, 1f32) - pb).len2())
					.min((ac * (ac.dot(pc) / ac.dot(ac)).clamp(0f32, 1f32) - pc).len2())
			} else {
				nor.dot(pa) * nor.dot(pa) / nor.dot(nor)
			}
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			assert_eq!(1f32, Triangle {
				a: Vec3(1.0, 0.0, 0.0),
				b: Vec3(0.0, 0.0, 0.0),
				c: Vec3(0.0, 1.0, 0.0)
			}.distance(&Vec3(0.0, 0.0, 1.0)))
		}
	}
}

pub mod quad {
	use super::*;
	
	#[derive(Copy, Clone, Debug)]
	pub struct Quad {
		pub a: Vec3<f32>,
		pub b: Vec3<f32>,
		pub c: Vec3<f32>,
		pub d: Vec3<f32>
	}
	
	impl Distance2<Vec3<f32>> for Quad {
		fn distance2(&self, p: &Vec3<f32>) -> f32 {
			let Quad { a, b, c, d } = *self;
			let p = *p;
			let ba = b - a;
			let pa = p - a;
			let cb = c - b;
			let pb = p - b;
			let dc = d - c;
			let pc = p - c;
			let ad = a - d;
			let pd = p - d;
			let nor = ba.cross(ad);
			
			if ba.cross(nor).dot(pa).sign() +
				cb.cross(nor).dot(pb).sign() +
				dc.cross(nor).dot(pc).sign() +
				ad.cross(nor).dot(pd).sign() < 3 {
				(ba * (ba.dot(pa) / ba.dot(ba)).clamp(0f32, 1f32) - pa).len2()
					.min((cb * (cb.dot(pb) / cb.dot(cb)).clamp(0f32, 1f32) - pb).len2())
					.min((dc * (dc.dot(pc) / dc.dot(dc)).clamp(0f32, 1f32) - pc).len2())
					.min((ad * (ad.dot(pd) / ad.dot(ad)).clamp(0f32, 1f32) - pd).len2())
			} else {
				nor.dot(pa) * nor.dot(pa) / nor.dot(nor)
			}
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn distance() {
			assert_eq!(1f32, Quad {
				a: Vec3(0.0, 0.0, 0.0),
				b: Vec3(1.0, 0.0, 0.0),
				c: Vec3(1.0, 1.0, 0.0),
				d: Vec3(0.0, 1.0, 0.0)
			}.distance(&Vec3(0.0, 0.0, 1.0)))
		}
	}
}

pub mod ops {
	use super::*;
	
	pub trait Distance<T> {
		fn distance(&self, _: &T) -> f32;
	}
	
	pub trait Distance2<T> {
		fn distance2(&self, _: &T) -> f32;
	}
	
	impl<T: Distance2<U>, U> Distance<U> for T {
		#[inline]
		fn distance(&self, v: &U) -> f32 {
			self.distance2(v).sqrt()
		}
	}
	
	pub trait Intersections<T> {
		fn intersections(&self, _: &T) -> u32;
		
		#[inline]
		fn intersects(&self, v: &T) -> bool {
			self.intersections(v) > 0
		}
	}
}

pub mod sdf {
	use super::*;
	
	#[inline]
	pub fn union(d1: f32, d2: f32) -> f32 {
		d1.min(d2)
	}
	
	#[inline]
	pub fn subtraction(d1: f32, d2: f32) -> f32 {
		(-d1).max(d2)
	}
	
	#[inline]
	pub fn intersection(d1: f32, d2: f32) -> f32 {
		d1.max(d2)
	}
	
	pub fn smooth_union(d1: f32, d2: f32, k: f32) -> f32 {
		let h = (0.5f32 + 0.5f32 * (d2 - d1) / k).clamp(0f32, 1f32);
		mix(d2, d1, h) - k * h * (1f32 - h)
	}
	
	pub fn smooth_subtraction(d1: f32, d2: f32, k: f32) -> f32 {
		let h = (0.5f32 - 0.5f32 * (d2 + d1) / k).clamp(0f32, 1f32);
		mix(d2, -d1, h) + k * h * (1f32 - h)
	}
	
	pub fn smooth_intersection(d1: f32, d2: f32, k: f32) -> f32 {
		let h = (0.5f32 - 0.5f32 * (d2 - d1) / k).clamp(0f32, 1f32);
		mix(d2, d1, h) + k * h * (1f32 - h)
	}
	
	pub struct Union<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Union<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			union(self.0.distance(point), self.1.distance(point))
		}
	}
	
	pub struct Subtraction<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Subtraction<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			subtraction(self.0.distance(point), self.1.distance(point))
		}
	}
	
	pub struct Intersection<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Intersection<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			intersection(self.0.distance(point), self.1.distance(point))
		}
	}
	
	pub struct SmoothUnion<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U, f32);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for SmoothUnion<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			smooth_union(self.0.distance(point), self.1.distance(point), self.2)
		}
	}
	
	pub struct SmoothSubtraction<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U, f32);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for SmoothSubtraction<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			smooth_subtraction(self.0.distance(point), self.1.distance(point), self.2)
		}
	}
	
	pub struct SmoothIntersection<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>>(T, U, f32);
	
	impl<T: Distance<Vec3<f32>>, U: Distance<Vec3<f32>>> Distance<Vec3<f32>> for SmoothIntersection<T, U> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			smooth_intersection(self.0.distance(point), self.1.distance(point), self.2)
		}
	}
	
	pub struct Round<T: Distance<Vec3<f32>>>(T, f32);
	
	impl<T: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Round<T> {
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			self.0.distance(point) - self.1
		}
	}
	
	pub struct Translation<T: Distance<Vec3<f32>>>(T, Vec3<f32>);
	
	impl<T: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Translation<T> {
		#[inline]
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			self.0.distance(&(*point + self.1))
		}
	}
	
	pub struct Rotation<T: Distance<Vec3<f32>>>(T, Quat32);
	
	impl<T: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Rotation<T> {
		#[inline]
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			self.0.distance(&Mat4::<f32>::from_rotation(self.1).transform_pos(*point))
		}
	}
	
	pub struct Scale<T: Distance<Vec3<f32>>>(T, f32);
	
	impl<T: Distance<Vec3<f32>>> Distance<Vec3<f32>> for Scale<T> {
		#[inline]
		fn distance(&self, point: &Vec3<f32>) -> f32 {
			self.0.distance(&(*point / self.1)) * self.1
		}
	}
}

pub mod poly {
	const EPS: f32 = 1E-6;
	const FRAC_1_3: f32 = 1.0 / 3.0;
	
	pub fn quadratic_roots(c0: f32, c1: f32, c2: f32) -> ([f32; 2], usize) {
		if c0.abs() < EPS { return ([-c2 / c1, 0.0], 1) }
		let discriminant = c1 * c1 - (4.0 * c0 * c2);
		([
			 (-c1 + discriminant.sqrt()) / (2.0 * c0),
			 (-c1 - discriminant.sqrt()) / (2.0 * c0)
		 ], if discriminant > EPS {
			2
		} else if discriminant < EPS {
			0
		} else {
			1
		})
	}
	
	#[inline]
	pub fn quadratic_roots_iter(c0: f32, c1: f32, c2: f32) -> QuadraticRootsIter {
		let (roots, count) = quadratic_roots(c0, c1, c2);
		QuadraticRootsIter {
			roots,
			count,
			index: 0
		}
	}
	
	pub struct QuadraticRootsIter {
		roots: [f32; 2],
		count: usize,
		index: usize
	}
	
	impl Iterator for QuadraticRootsIter {
		type Item = f32;
		
		fn next(&mut self) -> Option<Self::Item> {
			if self.index < self.count {
				let v = Some(self.roots[self.index]);
				self.index += 1;
				v
			} else {
				None
			}
		}
	}
	
	/// https://github.com/nical/lyon/blob/master/geom/src/utils.rs
	pub fn cubic_roots(c0: f32, mut c1: f32, mut c2: f32, mut c3: f32) -> ([f32; 3], usize) {
		// check if it's a quadratic equation to avoid div by 0
		if c0.abs() < EPS {
			let (roots, count) = quadratic_roots(c1, c2, c3);
			return ([
						roots[0],
						roots[1],
						0.0
					], count)
		}
		
		c1 /= c0;
		c2 /= c0;
		c3 /= c0;
		
		let d0 = (3.0 * c2 - c1 * c1) / 9.0;
		let d1 = (9.0 * c1 * c2 - 27.0 * c3 - 2.0 * c1 * c1 * c1) / 54.0;
		let d = d0 * d0 * d0 + d1 * d1;
		
		if d > EPS {
			let p = d1 + d.sqrt();
			let m = d1 - d.sqrt();
			let s = p.signum() * p.abs().powf(FRAC_1_3);
			let t = m.signum() * m.abs().powf(FRAC_1_3);
			
			if (s - t).abs() < EPS && (s + t) > EPS {
				([
					 -c1 * FRAC_1_3 + s + t,
					 -c1 * FRAC_1_3 - (s + t) * 0.5,
					 0.0
				 ], 2)
			} else {
				([
					 -c1 * FRAC_1_3 + s + t,
					 0.0,
					 0.0
				 ], 1)
			}
		} else {
			let theta = (d1 / (-d0 * d0 * d0).sqrt()).acos();
			let d0 = 2.0 * (-d0).sqrt();
			([
				 d0 * (theta * FRAC_1_3).cos() - c1 * FRAC_1_3,
				 d0 * ((theta + 2.0 * ::std::f32::consts::PI) * FRAC_1_3).cos() - c1 * FRAC_1_3,
				 d0 * ((theta + 4.0 * ::std::f32::consts::PI) * FRAC_1_3).cos() - c1 * FRAC_1_3,
			 ], 3)
		}
	}
	
	#[inline]
	pub fn cubic_roots_iter(c0: f32, c1: f32, c2: f32, c3: f32) -> CubicRootsIter {
		let (roots, count) = cubic_roots(c0, c1, c2, c3);
		CubicRootsIter {
			roots,
			count,
			index: 0
		}
	}
	
	pub struct CubicRootsIter {
		roots: [f32; 3],
		count: usize,
		index: usize
	}
	
	impl Iterator for CubicRootsIter {
		type Item = f32;
		
		fn next(&mut self) -> Option<Self::Item> {
			if self.index < self.count {
				let v = Some(self.roots[self.index]);
				self.index += 1;
				v
			} else {
				None
			}
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn quadratic() {
			assert_eq!(([0.5, -3.0], 2), quadratic_roots(2.0, 5.0, -3.0));
		}
		
		#[test]
		fn cubic() {
			let (roots, count) = cubic_roots(2.0, -4.0, -22.0, 24.0);
			assert_eq!(3, count);
			assert!((roots[0] - 4.0).abs() < EPS);
			assert!((roots[1] + 3.0).abs() < EPS);
			assert!((roots[2] - 1.0).abs() < EPS);
		}
	}
}

pub mod misc {
	#[inline]
	pub fn mix(x: f32, y: f32, a: f32) -> f32 {
		x * (1f32 - a) + y * a
	}
	
	pub trait Sign {
		fn sign(self) -> i32;
	}
	
	impl Sign for f32 {
		fn sign(self) -> i32 {
			if self < 0f32 {
				-1
			} else if self > 0f32 {
				1
			} else {
				0
			}
		}
	}
	
	impl Sign for f64 {
		fn sign(self) -> i32 {
			if self < 0f64 {
				-1
			} else if self > 0f64 {
				1
			} else {
				0
			}
		}
	}
	
	#[cfg(test)]
	mod tests {
		use super::*;
		
		#[test]
		fn test_mix() {
			assert_eq!(0.5, mix(0.0, 1.0, 0.5));
			assert_eq!(0.0, mix(-1.0, 1.0, 0.5));
		}
		
		#[test]
		fn sign() {
			assert_eq!(0, 0.0.sign());
			assert_eq!(-1, -1.5.sign());
			assert_eq!(1, 1.5.sign());
		}
	}
}