hyperbolic_drawing 1.0.0

#![allow(mixed_script_confusables)]
//#![warn(missing_docs)]



//! Library for manipulating and displaying elements of 2D hyperbolic space (ℍ²).
//! 
//! All implementation is done in terms of the [hyperboloid model],
//! without requirements for normalization.
//! 
//! The component along the axis of the hyperboloid is named `t`
//! (in reference to [spacetime algebra]).
//! However, contrary to popular convention, we use γ₀² = −1 and γ₁² = γ₂² = 1.
//! As such, elements on the hyperboloid are normalized to p² = −1.
//! 
//! [hyperboloid model]: https://en.wikipedia.org/wiki/Hyperboloid_model
//! [spacetime algebra]: https://en.wikipedia.org/wiki/Spacetime_algebra



use core::{fmt, ops};



#[cfg(feature = "plotters")]
pub mod plotters;
#[cfg(any(feature = "svg", feature = "svg_fmt"))]
pub mod svg_common;
#[cfg(feature = "svg")]
pub mod svg;
#[cfg(feature = "svg_fmt")]
pub mod svg_fmt;
//#[cfg(feature = "svg")]
//use svg::



/// Common trait for operations available on float types.
/// 
/// Additional items may be added at any time,
/// so it is advised to not provide additional implementations in a library.
pub trait Scalar: Sized + Copy + fmt::Debug + PartialOrd
	+ ops::Neg<Output = Self>
	+ ops::Add<Output = Self> + ops::Sub<Output = Self>
	+ ops::Mul<Output = Self> + ops::Div<Output = Self>
	+ ops::AddAssign + ops::SubAssign
	+ ops::MulAssign + ops::DivAssign
{
	/// The additive identity (0).
	const ZERO: Self;
	/// The multiplicative identity (1).
	const ONE: Self;
	/// The multiplicative involutive element (-1).
	const NEG_ONE: Self;
	const HALF: Self;
	const EPSILON: Self;
	const FEW_EPSILON: Self;
	/// A sufficiently large constant (the float implementations use `1e6`).
	/// 
	/// This is used to determine when a circle is, visually, approximately a line.
	const VERY_LARGE: Self;
	fn abs(self) -> Self;
	/// The maximum of two scalars.
	fn max(self, other: Self) -> Self;
	/// The reciprocal of a scalar.
	fn recip(self) -> Self;
	/// The square root of a scalar.
	fn sqrt(self) -> Self;
	/// The exponential of a scalar (x ↦ eˣ).
	fn exp(self) -> Self;
	/// The natural logarithm of a scalar.
	fn ln(self) -> Self;
	/// The sine and cosine of a scalar.
	fn sin_cos(self) -> (Self, Self);
	/// The arctangent of a scalar.
	fn atan(self) -> Self;
	/// The angle from the origin to a point (other, self).
	fn atan2(self, other: Self) -> Self;
	/// The hyperbolic sine and cosine of a scalar.
	fn sinh_cosh(self) -> (Self, Self);
	/// The inverse hyperbolic sine of a scalar.
	fn asinh(self) -> Self;
	/// The inverse hyperbolic cosine of a scalar.
	fn acosh(self) -> Self;
	/// Whether the scalar is not zero, infinite, subnormal, or NaN.
	/// 
	/// This is primarily relevant for floating point implementations.
	fn is_normal(self) -> bool;
	/// Whether the sign of the number is positive. The value returned for zero is irrelevant.
	fn is_sign_positive(self) -> bool;
}

macro_rules! scalar_impl {
	($f:ident) => { impl Scalar for $f {
		const ZERO: Self = 0.0;
		const ONE: Self = 1.0;
		const NEG_ONE: Self = -1.0;
		const HALF: Self = 0.5;
		const EPSILON: Self = Self::EPSILON;
		const FEW_EPSILON: Self = 2.0 * Self::EPSILON;
		const VERY_LARGE: Self = 1e6;
		fn abs(self) -> Self { self.abs() }
		fn max(self, other: Self) -> Self { self.max(other) }
		fn recip(self) -> Self { self.recip() }
		fn sqrt(self) -> Self { self.sqrt() }
		fn exp(self) -> Self { self.exp() }
		fn ln(self) -> Self { self.ln() }
		fn sin_cos(self) -> (Self, Self) { self.sin_cos() }
		fn atan(self) -> Self { self.atan() }
		fn atan2(self, other: Self) -> Self { self.atan2(other) }
		fn sinh_cosh(self) -> (Self, Self) { (self.sinh(), self.cosh()) }
		fn asinh(self) -> Self { self.asinh() }
		fn acosh(self) -> Self { self.acosh() }
		fn is_normal(self) -> bool { self.is_normal() }
		fn is_sign_positive(self) -> bool { self.is_sign_positive() }
	}};
}
//scalar_impl!(f16);
scalar_impl!(f32);
scalar_impl!(f64);
//scalar_impl!(f128);



/// An element of the Clifford algebra 𝒢(2,1) under the homogenous coordinates interpretation.
/// 
/// # Implementation rationale
/// 
/// We require implementers to specify not only the individual component values
/// but also constant `bool`s (e.g. [`Self::HAS_K`]) for whether these components will ever be present.
/// This is so that we can ensure components are optimized at compile time
/// even in cases where we could otherwise have e.g. `(f64::NAN * 0.0).is_nan()` preventing us from doing so.
/// 
/// Additionally, setting [`Self::UNITIZED`] provides faster calculation options in a few places.
pub trait Multivector: Sized + Copy {
	/// Scalar type for this multivector. The type of each component.
	type K: Scalar;
	/// Whether the k (scalar) component is ever non-zero for this type.
	const HAS_K: bool = false;
	/// Whether the t component is ever non-zero for this type.
	const HAS_T: bool = false;
	/// Whether the x component is ever non-zero for this type.
	const HAS_X: bool = false;
	/// Whether the y component is ever non-zero for this type.
	const HAS_Y: bool = false;
	/// Whether the tx component is ever non-zero for this type.
	const HAS_TX: bool = false;
	/// Whether the ty component is ever non-zero for this type.
	const HAS_TY: bool = false;
	/// Whether the xy component is ever non-zero for this type.
	const HAS_XY: bool = false;
	/// Whether the txy component is ever non-zero for this type.
	const HAS_TXY: bool = false;
	/// Whether elements of this type are always unitized up to rounding.
	/// i.e. `x.antidot(x)` always approximately equals ±1.
	const UNITIZED: bool;

	/// The k (scalar) component.
	#[inline(always)]
	fn k(self) -> Self::K { Self::K::ZERO }
	/// The t component.
	#[inline(always)]
	fn t(self) -> Self::K { Self::K::ZERO }
	/// The x component.
	#[inline(always)]
	fn x(self) -> Self::K { Self::K::ZERO }
	/// The y component.
	#[inline(always)]
	fn y(self) -> Self::K { Self::K::ZERO }
	/// The tx component.
	#[inline(always)]
	fn tx(self) -> Self::K { Self::K::ZERO }
	/// The ty component.
	#[inline(always)]
	fn ty(self) -> Self::K { Self::K::ZERO }
	/// The xy component.
	#[inline(always)]
	fn xy(self) -> Self::K { Self::K::ZERO }
	/// The txy component.
	#[inline(always)]
	fn txy(self) -> Self::K { Self::K::ZERO }

	/// Collect a value of `Self` from the relevant components.
	/// 
	/// Components are given as `[k, t, x, y, tx, ty, xy, txy]`.
	fn from_all_components(components: [Self::K; 8]) -> Self;

	/// The dot product `self ∙ rhs`.
	/// 
	/// See [the Rigid Geometric Algebra wiki] for the related Euclidean case.
	/// The main difference is that the zero-valued components are given a scale of −1 instead.
	/// 
	/// [the Rigid Geometric Algebra wiki]: https://rigidgeometricalgebra.org/wiki/index.php?title=Dot_products
	fn dot<M: Multivector<K = Self::K>>(self, rhs: M) -> Self::K {
		let mut res = Self::K::ZERO;
		if Self::HAS_K && M::HAS_K { res += self.k() * rhs.k() }
		if Self::HAS_T && M::HAS_T { res -= self.t() * rhs.t() }
		if Self::HAS_X && M::HAS_X { res += self.x() * rhs.x() }
		if Self::HAS_Y && M::HAS_Y { res += self.y() * rhs.y() }
		if Self::HAS_TX && M::HAS_TX { res += self.tx() * rhs.tx() }
		if Self::HAS_TY && M::HAS_TY { res += self.ty() * rhs.ty() }
		if Self::HAS_XY && M::HAS_XY { res -= self.xy() * rhs.xy() }
		if Self::HAS_TXY && M::HAS_TXY { res += self.txy() * rhs.txy() }
		res
	}
	/// The antidot product `self ∘ rhs`.
	/// 
	/// See [the Rigid Geometric Algebra wiki] for the related Euclidean case.
	/// The main difference is that the zero-valued components are given a scale of −1 instead.
	/// 
	/// [the Rigid Geometric Algebra wiki]: https://rigidgeometricalgebra.org/wiki/index.php?title=Dot_products
	fn antidot<M: Multivector<K = Self::K>>(self, rhs: M) -> Self::K {
		let mut res = Self::K::ZERO;
		if Self::HAS_K && M::HAS_K { res -= self.k() * rhs.k() }
		if Self::HAS_T && M::HAS_T { res -= self.t() * rhs.t() }
		if Self::HAS_X && M::HAS_X { res += self.x() * rhs.x() }
		if Self::HAS_Y && M::HAS_Y { res += self.y() * rhs.y() }
		if Self::HAS_TX && M::HAS_TX { res -= self.tx() * rhs.tx() }
		if Self::HAS_TY && M::HAS_TY { res -= self.ty() * rhs.ty() }
		if Self::HAS_XY && M::HAS_XY { res += self.xy() * rhs.xy() }
		if Self::HAS_TXY && M::HAS_TXY { res += self.txy() * rhs.txy() }
		res
	}
}
fn all_components<M: Multivector>(m: M) -> [M::K; 8] {
	[
		m.k(), m.t(), m.x(), m.y(),
		m.tx(), m.ty(), m.xy(), m.txy(),
	]
}
fn project_all_components<M1, M2>(m: M1) -> M2
where M1: Multivector, M2: Multivector<K = M1::K> {
	M2::from_all_components(all_components(m))
}
macro_rules! multivector_impl {
	(comp => $T:ident,) => {};
	(comp => $T:ident, k $($p:ident)*) => {
		const HAS_K: bool = true;
		#[inline(always)]
		fn k(self) -> K { self.k }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, t $($p:ident)*) => {
		const HAS_T: bool = true;
		#[inline(always)]
		fn t(self) -> K { self.t }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, x $($p:ident)*) => {
		const HAS_X: bool = true;
		#[inline(always)]
		fn x(self) -> K { self.x }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, y $($p:ident)*) => {
		const HAS_Y: bool = true;
		#[inline(always)]
		fn y(self) -> K { self.y }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, tx $($p:ident)*) => {
		const HAS_TX: bool = true;
		#[inline(always)]
		fn tx(self) -> K { self.tx }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, ty $($p:ident)*) => {
		const HAS_TY: bool = true;
		#[inline(always)]
		fn ty(self) -> K { self.ty }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, xy $($p:ident)*) => {
		const HAS_XY: bool = true;
		#[inline(always)]
		fn xy(self) -> K { self.xy }
		multivector_impl!(comp => $T, $($p)*);
	};
	(comp => $T:ident, txy $($p:ident)*) => {
		const HAS_TXY: bool = true;
		#[inline(always)]
		fn txy(self) -> K { self.txy }
		multivector_impl!(comp => $T, $($p)*);
	};
	(unit => weighted) => {
		const UNITIZED: bool = false;
	};
	(unit => unitized) => {
		const UNITIZED: bool = true;
	};
	($T:ident, $($p:ident)*, $comp:expr, $u:ident $(,)?) => {
		impl<K: Scalar> Multivector for $T<K> {
			type K = K;
			multivector_impl!(comp => $T, $($p)*);
			multivector_impl!(unit => $u);
			fn from_all_components(components: [K; 8]) -> Self {
				#[inline(always)]
				fn typeassist<S, K, F: FnOnce([K; 8]) -> S>(comp: F) -> F {
					comp
				}
				let comp = typeassist::<Self, K, _>(#[inline] $comp);
				comp(components)
			}
		}
	}
}

/// A common implementation for unitizing (scaling to weight ±1) a multivector.
pub trait Unitize: Multivector {
	/// The result of unitizing `Self`. This type may require unitization or not.
	type Output: Multivector<K = Self::K>;

	/// Unitize (scale to weight ±1) this multivector.
	/// 
	///  # Panics
	/// 
	/// Panics if the division factor √|`self` ∙ `self`| is not normal (as defined by [`Scalar::is_normal`]).
	/// i.e. if dividing by zero, infinity, or other edge cases of the number system where dividing the components is likely to exacerbate issues.
	/// 
	/// To handle these cases yourself, see [`Unitize::try_unitize`].
	fn unitize(self) -> Self::Output {
		let Some(res) = self.try_unitize() else {
			let q = self.antidot(self).abs().sqrt();
			panic!("failed to unitize multivector with length {q:?}");
		};
		res
	}

	/// Try to unitize (scale to weight ±1) this multivector.
	/// 
	/// Returns `None` only in the case that the division factor √|`self` ∙ `self`| is not normal (as defined by [`Scalar::is_normal`]).
	fn try_unitize(self) -> Option<Self::Output> {
		let q = self.antidot(self).abs().sqrt();
		q.is_normal().then(||
			Self::Output::from_all_components([
				if Self::HAS_K && Self::Output::HAS_K { self.k() / q } else { Self::K::ZERO },
				if Self::HAS_T && Self::Output::HAS_T { self.t() / q } else { Self::K::ZERO },
				if Self::HAS_X && Self::Output::HAS_X { self.x() / q } else { Self::K::ZERO },
				if Self::HAS_Y && Self::Output::HAS_Y { self.y() / q } else { Self::K::ZERO },
				if Self::HAS_TX && Self::Output::HAS_TX { self.tx() / q } else { Self::K::ZERO },
				if Self::HAS_TY && Self::Output::HAS_TY { self.ty() / q } else { Self::K::ZERO },
				if Self::HAS_XY && Self::Output::HAS_XY { self.xy() / q } else { Self::K::ZERO },
				if Self::HAS_TXY && Self::Output::HAS_TXY { self.txy() / q } else { Self::K::ZERO },
			])
		)
	}
}

/*
pub trait Classify<K: Scalar>: Multivector {
	type Object: Multivector<K = Self::K>;
	type Ideal: Multivector<K = Self::K>;
	type Dual: Multivector<K = Self::K>;

	fn classify(self) -> Classification<Self::Object, Self::Ideal, Self::Dual> {
		let q = self.dot(self);
		// note: this doesn't make sense in the general case (because signs flip)
	}
}
*/

macro_rules! weight_scale_impl {
	($T:ident, $U:ident) => {
		impl<K: Scalar> ops::Mul<K> for $T<K> {
			type Output = $U<K>;
			fn mul(self, other: K) -> Self::Output {
				Self::Output::from_all_components([
					if Self::HAS_K && Self::Output::HAS_K { self.k() * other } else { K::ZERO },
					if Self::HAS_T && Self::Output::HAS_T { self.t() * other } else { K::ZERO },
					if Self::HAS_X && Self::Output::HAS_X { self.x() * other } else { K::ZERO },
					if Self::HAS_Y && Self::Output::HAS_Y { self.y() * other } else { K::ZERO },
					if Self::HAS_TX && Self::Output::HAS_TX { self.tx() * other } else { K::ZERO },
					if Self::HAS_TY && Self::Output::HAS_TY { self.ty() * other } else { K::ZERO },
					if Self::HAS_XY && Self::Output::HAS_XY { self.xy() * other } else { K::ZERO },
					if Self::HAS_TXY && Self::Output::HAS_TXY { self.txy() * other } else { K::ZERO },
				])
			}
		}
	};
	($T:ident) => {
		weight_scale_impl!($T, $T);
		impl<K: Scalar> ops::MulAssign<K> for $T<K> {
			fn mul_assign(&mut self, other: K) {
				*self = *self * other;
			}
		}
	};
}
macro_rules! weight_add_impl {
	($T:ident + $U:ident -> $V:ident) => {
		impl<K: Scalar> ops::Add<$U<K>> for $T<K> {
			type Output = $V<K>;
			fn add(self, other: $U<K>) -> Self::Output {
				Self::Output::from_all_components([
					self.k() + other.k(),
					self.t() + other.t(),
					self.x() + other.x(),
					self.y() + other.y(),
					self.tx() + other.tx(),
					self.ty() + other.ty(),
					self.xy() + other.xy(),
					self.txy() + other.txy(),
				])
			}
		}
		impl<K: Scalar> ops::Sub<$U<K>> for $T<K> {
			type Output = $V<K>;
			fn sub(self, other: $U<K>) -> Self::Output {
				Self::Output::from_all_components([
					self.k() - other.k(),
					self.t() - other.t(),
					self.x() - other.x(),
					self.y() - other.y(),
					self.tx() - other.tx(),
					self.ty() - other.ty(),
					self.xy() - other.xy(),
					self.txy() - other.txy(),
				])
			}
		}
	};
	($T:ident -> $V:ident) => {
		weight_add_impl!($T + $T -> $V);
	};
	($T:ident += $U:ident) => {
		weight_add_impl!($T + $U -> $T);
		impl<K: Scalar> ops::AddAssign<$U<K>> for $T<K> {
			fn add_assign(&mut self, other: $U<K>) {
				*self = *self + other
			}
		}
		impl<K: Scalar> ops::SubAssign<$U<K>> for $T<K> {
			fn sub_assign(&mut self, other: $U<K>) {
				*self = *self - other
			}
		}
	};
	($T:ident) => {
		weight_add_impl!($T += $T);
	};
}



/// The kind of a [`Center`] as determined by [`Center::classify`].
pub enum CenterKind<K: Scalar> {
	/// A point center, used for rotations.
	Point(WeightedPoint<K>),
	/// A center at infinity, used for horolations.
	Ideal(IdealPoint<K>),
	/// An axis line, used for translations.
	Line(WeightedLine<K>),
}



/// A center of some transformation in ℍ²: a point, ideal point, or line.
#[derive(Clone, Copy, Debug)]
pub struct Center<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> Center<K> {
	/// Determine the kind of center.
	/// 
	/// In the case an operation gives a generic center object where you expect otherwise,
	/// see that individual operation for specific interpretations of each result.
	/// e.g. the intersection of two ultraparallel lines gives their unique common perpendicular.
	pub fn classify(self) -> CenterKind<K> {
		let t2 = self.t * self.t;
		let q = self.x * self.x + self.y * self.y - t2;
		if q.abs() <= t2 * K::FEW_EPSILON {
			CenterKind::Ideal(project_all_components(self))
		} else if q > K::ZERO {
			CenterKind::Line(project_all_components(self))
		} else {
			CenterKind::Point(project_all_components(self))
		}
	}

	/*
	fn klein(self) -> (K, K) {
		(self.x / self.t, self.y / self.t)
	}
	fn poincaré(self) -> (K, K) {
		let (kx, ky) = self.klein();
		let h2 = K::ONE - kx * kx - ky * ky;
		// unsure if this is correct in general,
		// though we do have the nice property that the klein coordinate of a line
		// is the center of the circle which represents it in the poincaré disk.
		if h2 < K::ZERO { return (kx, ky) };
		let h = h2.sqrt();
		// now in hemisphere
		let ih = K::recip(K::ONE + h);
		(kx * ih, ky * ih)
	}
	*/
}
impl<K: Scalar> From<WeightedPoint<K>> for Center<K> {
	fn from(value: WeightedPoint<K>) -> Self {
		value.into_center()
	}
}
impl<K: Scalar> From<Point<K>> for Center<K> {
	fn from(value: Point<K>) -> Self {
		value.into_center()
	}
}
impl<K: Scalar> From<IdealPoint<K>> for Center<K> {
	fn from(value: IdealPoint<K>) -> Self {
		value.into_center()
	}
}
impl<K: Scalar> From<WeightedLine<K>> for Center<K> {
	fn from(value: WeightedLine<K>) -> Self {
		value.into_center()
	}
}
impl<K: Scalar> From<Line<K>> for Center<K> {
	fn from(value: Line<K>) -> Self {
		value.into_center()
	}
}

multivector_impl!(Center, t x y, |[_, t, x, y, ..]| Self { t, x, y }, weighted);
weight_scale_impl!(Center);
weight_add_impl!(Center);
weight_add_impl!(Center += WeightedPoint);
weight_add_impl!(Center += Point);
weight_add_impl!(Center += IdealPoint);
weight_add_impl!(Center += WeightedLine);
weight_add_impl!(Center += Line);



/// A point in 2D hyperbolic space.
#[derive(Clone, Copy, Debug)]
pub struct WeightedPoint<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> WeightedPoint<K> {
	/// Upcast to a generic center.
	pub const fn into_center(self) -> Center<K> {
		let Self { t, x, y } = self;
		Center { t, x, y }
	}

	/// Reflect this point across another point.
	/// 
	/// As in Euclidean, this is the same as a 180° rotation around the point `p`.
	pub fn invert(self, p: WeightedPoint<K>) -> WeightedPoint<K> {
		let k = p.t * self.t - p.x * self.x - p.y * self.y;
		let tx = p.x * self.t - p.t * self.x;
		let ty = p.y * self.t - p.t * self.y;
		let xy = p.y * self.x - p.x * self.y;
		WeightedPoint {
			t: self.t * k - self.x * tx - self.y * ty,
			x: self.x * k - self.t * tx - self.y * xy,
			y: self.y * k - self.t * ty + self.x * xy,
		}
	}

	/// The Klein projection of this point (within the unit circle).
	pub fn klein(self) -> (K, K) {
		(self.x / self.t, self.y / self.t)
	}
	/// The Poincaré disk projection of this point (within the unit circle).
	pub fn poincaré(self) -> (K, K) {
		let (kx, ky) = self.klein();
		let h2 = K::ONE - kx * kx - ky * ky;
		let h = h2.max(K::ZERO).sqrt();
		// now in hemisphere
		let ih = K::recip(K::ONE + h);
		(kx * ih, ky * ih)
	}
}
impl<K: Scalar> ops::Mul for WeightedPoint<K> {
	type Output = Motor<K>;
	fn mul(self, rhs: Self) -> Self::Output {
		Motor {
			k: self.dot(rhs),
			tx: self.t * rhs.x - self.x * rhs.t,
			ty: self.t * rhs.y - self.y * rhs.t,
			xy: self.x * rhs.y - self.y * rhs.x,
		}
	}
}
impl<K: Scalar> From<Point<K>> for WeightedPoint<K> {
	fn from(value: Point<K>) -> Self {
		value.into_weighted()
	}
}
multivector_impl!(WeightedPoint, t x y, |[_, t, x, y, ..]| Self { t, x, y }, weighted);
impl<K: Scalar> Unitize for WeightedPoint<K> { type Output = Point<K>; }
weight_scale_impl!(WeightedPoint);
weight_add_impl!(WeightedPoint);
weight_add_impl!(WeightedPoint + Center -> Center);
weight_add_impl!(WeightedPoint + Point -> Center);
weight_add_impl!(WeightedPoint + IdealPoint -> Center);
weight_add_impl!(WeightedPoint + WeightedLine -> Center);
weight_add_impl!(WeightedPoint + Line -> Center);

/// A unitized point in hyperbolic space. Enables some additional calculations compared to [`WeightedPoint`].
#[derive(Clone, Copy, Debug)]
pub struct Point<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> Point<K> {
	/// Origin for the space, intersection of the x and y axes.
	pub const ORIGIN: Self = Self { t: K::ONE, x: K::ZERO, y: K::ZERO };

	/// Upcast to a generic center.
	pub const fn into_center(self) -> Center<K> {
		let Self { t, x, y } = self;
		Center { t, x, y }
	}
	/// Upcast to a point with weight 1.
	pub const fn into_weighted(self) -> WeightedPoint<K> {
		let Self { t, x, y } = self;
		WeightedPoint { t, x, y }
	}

	/// Distance of this point from the origin.
	/// 
	/// This is the distance in hyperbolic space, not the distance in the model.
	pub fn length(self) -> K {
		K::acosh(self.t)
	}
	/// Distance between this point and another.
	/// 
	/// This is the distance in hyperbolic space, not the distance in the model.
	pub fn dist(self, rhs: Self) -> K {
		let d = K::max(K::ONE, -self.dot(rhs));
		K::acosh(d)
	}

	/// Perform unitization again to prevent additional error buildup.
	pub fn unitize(mut self) -> Self {
		let q = self.t * self.t - self.x * self.x - self.y * self.y;
		let q = q.recip().sqrt();
		self.t *= q;
		self.x *= q;
		self.y *= q;
		self
	}

	/// A point at the given offset along the X axis.
	pub fn on_x(x: K) -> Self {
		let (sh, ch) = x.sinh_cosh();
		Self {
			t: ch,
			x: sh,
			y: K::ZERO,
		}
	}
	/// A point at the given offset along the Y axis.
	pub fn on_y(y: K) -> Self {
		let (sh, ch) = y.sinh_cosh();
		Self {
			t: ch,
			x: K::ZERO,
			y: sh,
		}
	}
	/// A point at the given polar coordinates.
	pub fn polar(r: K, θ: K) -> Self {
		let (sh, ch) = r.sinh_cosh();
		let (s, c) = θ.sin_cos();
		Self {
			t: ch,
			x: c * sh,
			y: s * sh,
		}
	}

	/// The midpoint between two points, the point halfway along the line between them.
	pub fn midpoint(self, other: Self) -> WeightedPoint<K> {
		if self.t.is_sign_positive() == other.t.is_sign_positive() {
			project_all_components(self + other)
		} else {
			project_all_components(self - other)
		}
	}
	/// The perpendicular bisector of the line segment between the points.
	pub fn bisector(self, other: Self) -> WeightedLine<K> {
		if self.t.is_sign_positive() == other.t.is_sign_positive() {
			project_all_components(self - other)
		} else {
			project_all_components(self + other)
		}
	}

	/// Reflect this point across another point.
	/// 
	/// As in Euclidean, this is the same as a 180° rotation around the point `p`.
	pub fn invert(self, p: Point<K>) -> Point<K> {
		let k = p.t * self.t - p.x * self.x - p.y * self.y;
		let tx = p.x * self.t - p.t * self.x;
		let ty = p.y * self.t - p.t * self.y;
		let xy = p.y * self.x - p.x * self.y;
		Point {
			t: self.t * k - self.x * tx - self.y * ty,
			x: self.x * k - self.t * tx - self.y * xy,
			y: self.y * k - self.t * ty + self.x * xy,
		}
	}

	/// The Klein projection of this point (within the unit circle).
	pub fn klein(self) -> (K, K) {
		(self.x / self.t, self.y / self.t)
	}
	/// The Poincaré disk projection of this point (within the unit circle).
	pub fn poincaré(self) -> (K, K) {
		let t1 = self.t + K::ONE;
		(self.x / t1, self.y / t1)
	}
}
impl<K: Scalar> ops::Mul for Point<K> {
	type Output = Motor<K>;
	fn mul(self, rhs: Self) -> Self::Output {
		Motor {
			k: self.dot(rhs),
			tx: self.t * rhs.x - self.x * rhs.t,
			ty: self.t * rhs.y - self.y * rhs.t,
			xy: self.x * rhs.y - self.y * rhs.x,
		}
	}
}
multivector_impl!(Point, t x y, |[_, t, x, y, ..]| Self { t, x, y }, unitized);
impl<K: Scalar> Unitize for Point<K> { type Output = Self; }
weight_scale_impl!(Point, WeightedPoint);
weight_add_impl!(Point -> Center);
weight_add_impl!(Point + Center -> Center);
weight_add_impl!(Point + WeightedPoint -> Center);
weight_add_impl!(Point + IdealPoint -> Center);
weight_add_impl!(Point + WeightedLine -> Center);
weight_add_impl!(Point + Line -> Center);

/// An ideal point (‘point at infinity’) in hyperbolic space.
#[derive(Clone, Copy, Debug)]
pub struct IdealPoint<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> IdealPoint<K> {
	/// Ideal point along the positive y axis.
	pub const NORTH: Self = Self { t: K::ONE, x: K::ZERO, y: K::ONE };
	/// Ideal point along the negative y axis.
	pub const SOUTH: Self = Self { t: K::ONE, x: K::ZERO, y: K::NEG_ONE };
	/// Ideal point along the positive x axis.
	pub const EAST: Self = Self { t: K::ONE, x: K::ONE, y: K::ZERO };
	/// Ideal point along the negative x axis.
	pub const WEST: Self = Self { t: K::ONE, x: K::NEG_ONE, y: K::ZERO };

	/// Upcast to a generic center.
	pub const fn into_center(self) -> Center<K> {
		let Self { t, x, y } = self;
		Center { t, x, y }
	}

	/// An ideal point at infinity in the given direction.
	pub fn polar(θ: K) -> Self {
		let (s, c) = θ.sin_cos();
		Self {
			t: K::ONE,
			x: c,
			y: s,
		}
	}

	/// Angle of this ideal point relative to the X axis.
	/// 
	/// # Examples
	/// 
	/// ```
	/// # type IdealPoint = hyperbolic_drawing::IdealPoint::<f32>;
	/// assert_eq!(IdealPoint::EAST.angle(), 0.0);
	/// assert_eq!(IdealPoint::NORTH.angle(), std::f32::consts::FRAC_PI_2);
	/// assert_eq!(IdealPoint::WEST.angle(), std::f32::consts::PI);
	/// assert_eq!(IdealPoint::polar(1.0).angle(), 1.0);
	/// ```
	pub fn angle(self) -> K {
		if self.t >= K::ZERO {
			K::atan2(self.y, self.x)
		} else {
			K::atan2(-self.y, -self.x)
		}
	}
}
multivector_impl!(IdealPoint, t x y, |[_, t, x, y, ..]| Self { t, x, y }, weighted);
weight_scale_impl!(IdealPoint);
weight_add_impl!(IdealPoint -> Center);
weight_add_impl!(IdealPoint + Center -> Center);
weight_add_impl!(IdealPoint + WeightedPoint -> Center);
weight_add_impl!(IdealPoint + Point -> Center);
weight_add_impl!(IdealPoint + WeightedLine -> Center);
weight_add_impl!(IdealPoint + Line -> Center);

/// A line in hyperbolic geometry.
#[derive(Clone, Copy, Debug)]
pub struct WeightedLine<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> WeightedLine<K> {
	/// Upcast to a generic center.
	pub const fn into_center(self) -> Center<K> {
		let Self { t, x, y } = self;
		Center { t, x, y }
	}

	/// The line which connects the two points.
	pub fn between(a: WeightedPoint<K>, b: WeightedPoint<K>) -> Self {
		let xt = a.x * b.t - a.t * b.x;
		let ty = a.t * b.y - a.y * b.t;
		let xy = a.x * b.y - a.y * b.x;
		Self { t: xy, x: ty, y: xt }
	}

	/// Reflect a point across this line.
	pub fn reflect(self, p: WeightedPoint<K>) -> WeightedPoint<K> {
		let k = p.t * self.t - p.x * self.x - p.y * self.y;
		let tx = p.x * self.t - p.t * self.x;
		let ty = p.y * self.t - p.t * self.y;
		let xy = p.y * self.x - p.x * self.y;
		WeightedPoint {
			t: self.t * k - self.x * tx - self.y * ty,
			x: self.x * k - self.t * tx - self.y * xy,
			y: self.y * k - self.t * ty + self.x * xy,
		}
	}

	/// Radius of the arc representing this line in the Poincaré disk projection.
	fn poincaré_radius(self) -> K {
		let r2 = self.x * self.x + self.y * self.y;
		let r2 = r2 / (self.t * self.t);
		(r2 - K::ONE).max(K::ZERO).sqrt()
	}
	/*
	fn poincaré_center(self) -> (K, K) {
		(self.x / self.t, self.y / self.t)
	}
	*/
}
impl<K: Scalar> From<Line<K>> for WeightedLine<K> {
	fn from(value: Line<K>) -> Self {
		let Line { t, x, y } = value;
		Self { t, x, y }
	}
}

multivector_impl!(WeightedLine, t x y, |[_, t, x, y, ..]| Self { t, x, y }, weighted);
impl<K: Scalar> Unitize for WeightedLine<K> { type Output = Line<K>; }
weight_scale_impl!(WeightedLine);
weight_add_impl!(WeightedLine -> Center);
weight_add_impl!(WeightedLine + Center -> Center);
weight_add_impl!(WeightedLine + WeightedPoint -> Center);
weight_add_impl!(WeightedLine + Point -> Center);
weight_add_impl!(WeightedLine + IdealPoint -> Center);
weight_add_impl!(WeightedLine + Line -> Center);

/// A unitized line in hyperbolic geometry. Enables some additional calculations compared to [`WeightedLine`].
#[derive(Clone, Copy, Debug)]
pub struct Line<K: Scalar = f64> {
	pub t: K,
	pub x: K,
	pub y: K,
}
impl<K: Scalar> Line<K> {
	/// Ultraideal point denoting a horizontal line at the origin
	pub const X: Self = Self { t: K::ZERO, x: K::ZERO, y: K::NEG_ONE };
	/// Ultraideal point denoting a vertical line at the origin
	pub const Y: Self = Self { t: K::ZERO, x: K::ONE, y: K::ZERO };

	/// Upcast to a generic center.
	pub const fn into_center(self) -> Center<K> {
		let Self { t, x, y } = self;
		Center { t, x, y }
	}
	/// Upcast to a line with weight 1.
	pub const fn into_weighted(self) -> WeightedLine<K> {
		let Self { t, x, y } = self;
		WeightedLine { t, x, y }
	}
	/*
	pub fn into_line(self) -> Line<K> {
		Line::from(self)
	}
	*/

	/// A line perpendicular to the X axis at the given offset.
	pub fn on_x(x: K) -> Self {
		let (sh, ch) = x.sinh_cosh();
		Self {
			t: sh,
			x: ch,
			y: K::ZERO,
		}
	}
	/// A line perpendicular to the Y axis at the given offset.
	pub fn on_y(y: K) -> Self {
		let (sh, ch) = y.sinh_cosh();
		Self {
			t: sh,
			x: K::ZERO,
			y: ch,
		}
	}
	/// A line closest to the origin at the given polar coordinates.
	pub fn polar(r: K, θ: K) -> Self {
		let (sh, ch) = r.sinh_cosh();
		let (s, c) = θ.sin_cos();
		Self {
			t: sh,
			x: c * ch,
			y: s * ch,
		}
	}

	/// Reflect a point across this line.
	pub fn reflect(self, p: Point<K>) -> Point<K> {
		let k = p.t * self.t - p.x * self.x - p.y * self.y;
		let tx = p.x * self.t - p.t * self.x;
		let ty = p.y * self.t - p.t * self.y;
		let xy = p.y * self.x - p.x * self.y;
		Point {
			t: self.t * k - self.x * tx - self.y * ty,
			x: self.x * k - self.t * tx - self.y * xy,
			y: self.y * k - self.t * ty + self.x * xy,
		}
	}

	/// The line halfway between this and another line.
	pub fn midline(self, other: Self) -> WeightedLine<K> {
		if self.dot(other).is_sign_positive() {
			project_all_components(self + other)
		} else {
			project_all_components(self - other)
		}
	}
	/// A line perpendicular to the two given at their closest points to eachother.
	pub fn bisector(self, other: Self) -> WeightedLine<K> {
		if self.dot(other).is_sign_positive() {
			project_all_components(self - other)
		} else {
			project_all_components(self + other)
		}
	}

	/*
	fn poincaré_radius(self) -> K {
		let r2 = self.x * self.x + self.y * self.y;
		let r2 = r2 / (self.t * self.t);
		(r2 - K::ONE).max(K::ZERO).sqrt()
	}
	fn poincaré_center(self) -> (K, K) {
		(self.x / self.t, self.y / self.t)
	}
	*/
}

multivector_impl!(Line, t x y, |[_, t, x, y, ..]| Self { t, x, y }, unitized);
impl<K: Scalar> Unitize for Line<K> { type Output = Self; }
weight_scale_impl!(Line, WeightedLine);
weight_add_impl!(Line -> Center);
weight_add_impl!(Line + Center -> Center);
weight_add_impl!(Line + WeightedPoint -> Center);
weight_add_impl!(Line + Point -> Center);
weight_add_impl!(Line + IdealPoint -> Center);
weight_add_impl!(Line + WeightedLine -> Center);



#[derive(Clone, Copy, Debug)]
struct CenterBivector<K: Scalar = f64> {
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> CenterBivector<K> {
	pub fn at(c: Center<K>) -> Self {
		let Center { t, x, y } = c;
		Self { tx: -y, ty: x, xy: t }
	}
}

/*
	pub fn classify(self) -> Classification<Center<K>, Horocenter<K>, Line<K>> {
		let xy2 = self.xy * self.xy;
		let q = self.tx * self.tx + self.ty * self.ty - xy2;
		if q.abs() <= xy2 * K::FEW_EPSILON {
			Classification::Ideal(project_all_components(self))
		} else if q > K::ZERO {
			Classification::Dual(project_all_components(self))
		} else {
			Classification::Object(project_all_components(self))
		}
	}
multivector_impl!(AbstractCenter, tx ty xy, |[_, _, _, _, tx, ty, xy, _]| Self { tx, ty, xy }, weighted);
weight_scale_impl!(AbstractCenter);
weight_add_impl!(AbstractCenter);
weight_add_impl!(AbstractCenter + Center -> AbstractCenter);
weight_add_impl!(AbstractCenter + Horocenter -> AbstractCenter);
weight_add_impl!(AbstractCenter + Line -> AbstractCenter);

#[derive(Clone, Copy, Debug)]
pub struct WeightedCenter<K: Scalar> {
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> WeightedCenter<K> {
	pub fn at(c: WeightedPoint<K>) -> Self {
		let WeightedPoint { t, x, y } = c;
		Self { tx: -y, ty: x, xy: t }
	}
}
#[derive(Clone, Copy, Debug)]
pub struct Center<K: Scalar> {
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> Center<K> {
	pub fn at(c: Point<K>) -> Self {
		let Point { t, x, y } = c;
		Self { tx: -y, ty: x, xy: t }
	}
}
multivector_impl!(Center, tx ty xy, |[_, _, _, _, tx, ty, xy, _]| Self { tx, ty, xy }, weighted);

#[derive(Clone, Copy, Debug)]
pub struct Horocenter<K: Scalar = f64> {
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> Horocenter<K> {
	pub fn at_right(r: K) -> Self {
		let re = r.exp();
		Self {
			tx: K::ZERO,
			ty: re,
			xy: re,
		}
	}
	pub fn at_up(u: K) -> Self {
		let ue = u.exp();
		Self {
			tx: -ue,
			ty: K::ZERO,
			xy: ue,
		}
	}
	pub fn at_left(l: K) -> Self {
		let le = l.exp();
		Self {
			tx: K::ZERO,
			ty: -le,
			xy: le,
		}
	}
	pub fn at_down(d: K) -> Self {
		let de = d.exp();
		Self {
			tx: de,
			ty: K::ZERO,
			xy: de,
		}
	}
	pub fn polar(r: K, θ: K) -> Self {
		let re = r.exp();
		let (s, c) = θ.sin_cos();
		Self {
			tx: -s * re,
			ty: c * re,
			xy: re,
		}
	}

	pub fn ideal_center(self) -> K {
		K::atan2(-self.tx, self.ty)
	}
	pub fn surface_distance(self) -> K {
		K::ln(self.xy)
	}
}
multivector_impl!(Horocenter, tx ty xy, |[_, _, _, _, tx, ty, xy, _]| Self { tx, ty, xy }, weighted);

#[derive(Clone, Copy, Debug)]
pub struct Line<K: Scalar = f64> {
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> Line<K> {
	pub fn between(u: WeightedPoint<K>, v: WeightedPoint<K>) -> Self {
		Self {
			tx: u.t * v.x - u.x * v.t,
			ty: u.t * v.y - u.y * v.t,
			xy: u.x * v.y - u.y * v.x,
		}
	}

	pub fn normalize(mut self) -> Self {
		let q = self.tx * self.tx + self.ty * self.ty - self.xy * self.xy;
		let q = q.recip().sqrt();
		self.tx *= q;
		self.ty *= q;
		self.xy *= q;
		self
	}

	/*
	fn limit_center(self) -> K {
		K::atan2(-self.tx, self.ty)
	}
	fn limit_spread_radius(self) -> (K, K) {
		let r2 = self.tx * self.tx + self.ty * self.ty;
		let r2 = r2 / (self.xy * self.xy);
		//let radius = (r2 - 1.0).recip().sqrt();
		//let spread = radius.atan();
		let radius = (r2 - K::ONE).sqrt();
		let spread = radius.atan();
		(spread, radius)
	}
	*/
	fn poincaré_radius(self) -> K {
		let r2 = self.tx * self.tx + self.ty * self.ty;
		let r2 = r2 / (self.xy * self.xy);
		(r2 - K::ONE).sqrt()
	}
}
impl<K: Scalar> From<DualPoint<K>> for Line<K> {
	fn from(value: DualPoint<K>) -> Self {
		let DualPoint { t, x, y } = value;
		Self { tx: -y, ty: x, xy: t }
	}
}
multivector_impl!(Line, tx ty xy, |[_, _, _, _, tx, ty, xy, _]| Self { tx, ty, xy }, weighted);
*/

/// Precomputed values for performing a motion on other objects.
#[derive(Clone, Copy, Debug)]
pub struct Motor<K: Scalar = f64> {
	pub k: K,
	pub tx: K,
	pub ty: K,
	pub xy: K,
}
impl<K: Scalar> Motor<K> {
	/// The identity motor. Has no effect on any objects.
	pub const IDENTITY: Self = Self { k: K::ONE, tx: K::ZERO, ty: K::ZERO, xy: K::ZERO };

	/// Apply this motion to the given point.
	pub fn apply(self, v: WeightedPoint<K>) -> WeightedPoint<K> {
		/*
		let t = v.t * self.k - v.x * self.tx - v.y * self.ty;
		let x = v.x * self.k - v.t * self.tx - v.y * self.xy;
		let y = v.y * self.k - v.t * self.ty + v.x * self.xy;
		let txy = v.t * self.xy - v.x * self.ty + v.y * self.tx;
		Point {
			t: self.k * t - self.tx * x - self.ty * y + self.xy * txy,
			x: -self.tx * t + self.k * x - self.xy * y + self.ty * txy,
			y: -self.ty * t + self.xy * x + self.k * y - self.tx * txy,
			//txy: -self.xy * t + self.ty * x - self.tx * y + self.k * txy,
		}
		*/
		let t = v.t * self.k + v.x * self.tx + v.y * self.ty;
		let x = v.x * self.k + v.t * self.tx + v.y * self.xy;
		let y = v.y * self.k + v.t * self.ty - v.x * self.xy;
		let txy = -v.t * self.xy + v.x * self.ty - v.y * self.tx;
		WeightedPoint {
			t: self.k * t + self.tx * x + self.ty * y - self.xy * txy,
			x: self.tx * t + self.k * x + self.xy * y - self.ty * txy,
			y: self.ty * t - self.xy * x + self.k * y + self.tx * txy,
			//txy: self.xy * t - self.ty * x + self.tx * y + self.k * txy,
		}
	}

	/// A motor representing a horolation moving the origin by distance 2θ around the given ideal point.
	/// 
	/// Very slightly faster to produce than [`Motor::horolation`].
	pub fn horolation2(c: IdealPoint<K>, θ: K) -> Self {
		let c = CenterBivector::at(c.into_center());
		Self {
			k: K::NEG_ONE,
			tx: θ * c.tx,
			ty: θ * c.ty,
			xy: θ * c.xy,
		}
	}
	/// A motor representing a horolation moving the origin by distance θ around the given ideal point.
	pub fn horolation(c: IdealPoint<K>, θ: K) -> Self {
		Self::horolation2(c, K::HALF * θ)
	}
	/// A motor representating a movement along the line connecting `u` and `v` by distance *twice* that between `u` and `v`.
	/// 
	/// Somewhat faster to produce than [`Motor::along`].
	pub fn along2(u: Point<K>, v: Point<K>) -> Self {
		u * v
	}
	/// A motor representating a movement along the line connecting `u` and `v` by the distance between `u` and `v`.
	pub fn along(u: Point<K>, v: Point<K>) -> Self {
		u * u.midpoint(v).unitize()
	}

	pub fn normalize(mut self) -> Self {
		let q = self.k * self.k - self.tx * self.tx - self.ty * self.ty + self.xy * self.xy;
		let q = q.recip().sqrt();
		self.k *= q;
		self.tx *= q;
		self.ty *= q;
		self.xy *= q;
		self
	}
}
impl<K: Scalar> ops::Mul for Motor<K> {
	type Output = Self;
	fn mul(self, rhs: Self) -> Self::Output {
		Self {
			k: self.k * rhs.k + self.tx * rhs.tx + self.ty * rhs.ty - self.xy * rhs.xy,
			tx: self.k * rhs.tx + self.tx * rhs.k - self.ty * rhs.xy + self.xy * rhs.ty,
			ty: self.k * rhs.ty + self.tx * rhs.xy + self.ty * rhs.k - self.xy * rhs.tx,
			xy: self.k * rhs.xy + self.tx * rhs.ty - self.ty * rhs.tx + self.xy * rhs.k,
		}
	}
}



#[cfg(test)]
mod tests {
	use super::*;

	/// that a motor and its inverse compose to the identity
	#[test]
	fn compose_inverses() {
		fn dist2ish(m1: Motor, m2: Motor) -> f64 {
			let dk = m1.k - m2.k;
			let dtx = m1.tx - m2.tx;
			let dty = m1.ty - m2.ty;
			let dxy = m1.xy - m2.xy;
			dk * dk + dtx * dtx + dty * dty + dxy * dxy
		}
		fn assert_inverses(m1: Motor, m2: Motor) {
			let c = (m1 * m2).normalize();
			const NEG_IDENTITY: Motor = Motor { k: -1.0, tx: 0.0, ty: 0.0, xy: 0.0 };
			let d1 = dist2ish(Motor::IDENTITY, c);
			let d2 = dist2ish(NEG_IDENTITY, c);
			assert!(d1 <= 1e-5 || d2 <= 1e-5);
		}
		let p1 = Point::on_x(2.0);
		let p2 = Point::on_y(3.0);
		assert_inverses(Motor::along2(p1, p2), Motor::along2(p2, p1));
		let h = IdealPoint::polar(1.0);
		assert_inverses(Motor::horolation2(h, 5.0), Motor::horolation2(h, -5.0));
	}

	/// that applying two motors in sequence is equivalent to applying their composition
	#[test]
	fn apply_action_of_compose() {
		let m1 = Motor::horolation(IdealPoint::polar(5.0), 6.0);
		let m2 = Motor::along2(Point::on_x(2.0), Point::on_y(3.0));
		for i in 1..=5 {
			for j in 0..5 {
				let p = Point::polar(i as f64 / 5.0, j as f64).into_weighted();
				let composed = (m1 * m2).apply(p);
				let applied = m1.apply(m2.apply(p));
				let dist = dbg!(composed.unitize().dist(applied.unitize()));
				assert!(dist <= 1e-2);
			}
		}
	}
}