diff --git a/.travis.yml b/.travis.yml index cc196c3..f34d104 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,6 +11,9 @@ env: - NO_STD=1 matrix: + allow_failures: + - rust: nightly + env: exclude: - rust: beta env: NO_STD=1 diff --git a/CHANGELOG.md b/CHANGELOG.md index cd1ca2b..fe6b1f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,11 @@ documented here. This project adheres to [Semantic Versioning](http://semver.org/). +## [0.9] + * Add the `ComplexField` trait. + * Rename the `Real` trait to `RealField` for more clarity. + * Rename `Real::powf` to `ComplexField::pow`. + ## [0.5.2] * Add feature-gated support for the `decimal` crate. * Add attributes to automatically generate quickcheck checks. diff --git a/alga/Cargo.toml b/alga/Cargo.toml index 0e4d2b7..41efc51 100644 --- a/alga/Cargo.toml +++ b/alga/Cargo.toml @@ -9,6 +9,7 @@ readme = "README.md" repository = "https://github.com/rustsim/alga" documentation = "https://docs.rs/alga" license = "Apache-2.0" +edition = "2018" [lib] name = "alga" diff --git a/alga/src/general/complex.rs b/alga/src/general/complex.rs new file mode 100644 index 0000000..988411b --- /dev/null +++ b/alga/src/general/complex.rs @@ -0,0 +1,935 @@ +use num::{FromPrimitive, Num, NumAssign, Signed, Zero, One}; +use std::any::Any; +use std::fmt::{Debug, Display}; +use std::ops::Neg; +use std::{f32, f64}; + +use crate::general::{Field, SubsetOf, SupersetOf, RealField, MeetSemilattice, JoinSemilattice}; + +#[cfg(not(feature = "std"))] +use libm::F32Ext; +#[cfg(not(feature = "std"))] +use libm::F64Ext; +#[cfg(not(feature = "std"))] +use num; +//#[cfg(feature = "decimal")] +//use decimal::d128; + +#[allow(missing_docs)] + +/// Trait shared by all complex fields and its subfields (like real numbers). +/// +/// Complex numbers are equipped with functions that are commonly used on complex numbers and reals. +/// The results of those functions only have to be approximately equal to the actual theoretical values. +// FIXME: SubsetOf should be removed when specialization will be supported by rustc. This will +// allow a blanket impl: impl SubsetOf for T { ... } +// NOTE: make all types debuggable/'static/Any ? This seems essential for any kind of generic programming. +pub trait ComplexField: + SubsetOf + + SupersetOf + + Field + + Copy + + Num + + NumAssign + + FromPrimitive + + Neg + + MeetSemilattice + + JoinSemilattice +// + RelativeEq +// + UlpsEq + + Send + + Sync + + Any + + 'static + + Debug + + Display +{ + /// Type of the coefficients of a complex number. + type RealField: RealField; + + /// Builds a pure-real complex number from the given value. + fn from_real(re: Self::RealField) -> Self; + + /// The real part of this complex number. + fn real(self) -> Self::RealField; + + /// The imaginary part of this complex number. + fn imaginary(self) -> Self::RealField; + + /// The modulus of this complex number. + fn modulus(self) -> Self::RealField; + + /// The squared modulus of this complex number. + fn modulus_squared(self) -> Self::RealField; + + /// The argument of this complex number. + fn argument(self) -> Self::RealField; + + /// The sum of the absolute value of this complex number's real and imaginary part. + fn norm1(self) -> Self::RealField; + + /// Multiplies this complex number by `factor`. + fn scale(self, factor: Self::RealField) -> Self; + + /// Multiplies this complex number by `factor`. + fn unscale(self, factor: Self::RealField) -> Self; + + /// The polar form of this complex number: (modulus, arg) + fn to_polar(self) -> (Self::RealField, Self::RealField) { + (self.modulus(), self.argument()) + } + + /// The exponential form of this complex number: (modulus, e^{i arg}) + fn to_exp(self) -> (Self::RealField, Self) { + let m = self.modulus(); + + if !m.is_zero() { + (m, self.unscale(m)) + } else { + (Self::RealField::zero(), Self::one()) + } + } + + /// The exponential part of this complex number: `self / self.modulus()` + fn signum(self) -> Self { + self.to_exp().1 + } + + + fn floor(self) -> Self; + fn ceil(self) -> Self; + fn round(self) -> Self; + fn trunc(self) -> Self; + fn fract(self) -> Self; + fn mul_add(self, a: Self, b: Self) -> Self; + + /// The absolute value of this complex number: `self / self.signum()`. + /// + /// This is equivalent to `self.modulus()`. + fn abs(self) -> Self::RealField; + + /// Computes (self.conjugate() * self + other.conjugate() * other).sqrt() + fn hypot(self, other: Self) -> Self::RealField; + + fn recip(self) -> Self; + fn conjugate(self) -> Self; + fn sin(self) -> Self; + fn cos(self) -> Self; + fn sin_cos(self) -> (Self, Self); + #[inline] + fn sinh_cosh(self) -> (Self, Self) { + (self.sinh(), self.cosh()) + } + fn tan(self) -> Self; + fn asin(self) -> Self; + fn acos(self) -> Self; + fn atan(self) -> Self; + fn sinh(self) -> Self; + fn cosh(self) -> Self; + fn tanh(self) -> Self; + fn asinh(self) -> Self; + fn acosh(self) -> Self; + fn atanh(self) -> Self; + + fn is_finite(&self) -> bool; + + /// Cardinal sine + #[inline] + fn sinc(self) -> Self { + if self.is_zero() { + Self::one() + } else { + self.sin() / self + } + } + + #[inline] + fn sinhc(self) -> Self { + if self.is_zero() { + Self::one() + } else { + self.sinh() / self + } + } + + /// Cardinal cos + #[inline] + fn cosc(self) -> Self { + if self.is_zero() { + Self::one() + } else { + self.cos() / self + } + } + + #[inline] + fn coshc(self) -> Self { + if self.is_zero() { + Self::one() + } else { + self.cosh() / self + } + } + + fn log(self, base: Self::RealField) -> Self; + fn log2(self) -> Self; + fn log10(self) -> Self; + fn ln(self) -> Self; + fn ln_1p(self) -> Self; + fn sqrt(self) -> Self; + fn try_sqrt(self) -> Option; + fn exp(self) -> Self; + fn exp2(self) -> Self; + fn exp_m1(self) -> Self; + fn powi(self, n: i32) -> Self; + fn powf(self, n: Self::RealField) -> Self; + fn powc(self, n: Self) -> Self; + fn cbrt(self) -> Self; +} + +macro_rules! impl_complex( + ($($T:ty, $M:ident, $libm: ident);*) => ($( + impl ComplexField for $T { + type RealField = $T; + + #[inline] + fn from_real(re: Self::RealField) -> Self { + re + } + + #[inline] + fn real(self) -> Self::RealField { + self + } + + #[inline] + fn imaginary(self) -> Self::RealField { + Self::zero() + } + + #[inline] + fn norm1(self) -> Self::RealField { + $libm::abs(self) + } + + #[inline] + fn modulus(self) -> Self::RealField { + $libm::abs(self) + } + + #[inline] + fn modulus_squared(self) -> Self::RealField { + self * self + } + + #[inline] + fn argument(self) -> Self::RealField { + if self >= Self::zero() { + Self::zero() + } else { + Self::pi() + } + } + + #[inline] + fn to_exp(self) -> (Self, Self) { + if self >= Self::zero() { + (self, Self::one()) + } else { + (-self, -Self::one()) + } + } + + #[inline] + fn recip(self) -> Self { + $M::recip(self) + } + + #[inline] + fn conjugate(self) -> Self { + self + } + + #[inline] + fn scale(self, factor: Self::RealField) -> Self { + self * factor + } + + #[inline] + fn unscale(self, factor: Self::RealField) -> Self { + self / factor + } + + #[inline] + fn floor(self) -> Self { + $libm::floor(self) + } + + #[inline] + fn ceil(self) -> Self { + $libm::ceil(self) + } + + #[inline] + fn round(self) -> Self { + $libm::round(self) + } + + #[inline] + fn trunc(self) -> Self { + $libm::trunc(self) + } + + #[inline] + fn fract(self) -> Self { + $libm::fract(self) + } + + #[inline] + fn abs(self) -> Self { + $libm::abs(self) + } + + #[inline] + fn signum(self) -> Self { + Signed::signum(&self) + } + + #[inline] + fn mul_add(self, a: Self, b: Self) -> Self { + $libm::mul_add(self, a, b) + } + + #[cfg(feature = "std")] + #[inline] + fn powi(self, n: i32) -> Self { + self.powi(n) + } + + #[cfg(not(feature = "std"))] + #[inline] + fn powi(self, n: i32) -> Self { + // FIXME: is there a more accurate solution? + $libm::powf(self, n as $T) + } + + #[inline] + fn powf(self, n: Self) -> Self { + $libm::powf(self, n) + } + + #[inline] + fn powc(self, n: Self) -> Self { + // Same as powf. + $libm::powf(self, n) + } + + #[inline] + fn sqrt(self) -> Self { + $libm::sqrt(self) + } + + #[inline] + fn try_sqrt(self) -> Option { + if self >= Self::zero() { + Some($libm::sqrt(self)) + } else { + None + } + } + + #[inline] + fn exp(self) -> Self { + $libm::exp(self) + } + + #[inline] + fn exp2(self) -> Self { + $libm::exp2(self) + } + + + #[inline] + fn exp_m1(self) -> Self { + $libm::exp_m1(self) + } + + #[inline] + fn ln_1p(self) -> Self { + $libm::ln_1p(self) + } + + #[inline] + fn ln(self) -> Self { + $libm::ln(self) + } + + #[inline] + fn log(self, base: Self) -> Self { + $libm::log(self, base) + } + + #[inline] + fn log2(self) -> Self { + $libm::log2(self) + } + + #[inline] + fn log10(self) -> Self { + $libm::log10(self) + } + + #[inline] + fn cbrt(self) -> Self { + $libm::cbrt(self) + } + + #[inline] + fn hypot(self, other: Self) -> Self::RealField { + $libm::hypot(self, other) + } + + #[inline] + fn sin(self) -> Self { + $libm::sin(self) + } + + #[inline] + fn cos(self) -> Self { + $libm::cos(self) + } + + #[inline] + fn tan(self) -> Self { + $libm::tan(self) + } + + #[inline] + fn asin(self) -> Self { + $libm::asin(self) + } + + #[inline] + fn acos(self) -> Self { + $libm::acos(self) + } + + #[inline] + fn atan(self) -> Self { + $libm::atan(self) + } + + #[inline] + fn sin_cos(self) -> (Self, Self) { + $libm::sin_cos(self) + } + +// #[inline] +// fn exp_m1(self) -> Self { +// $libm::exp_m1(self) +// } +// +// #[inline] +// fn ln_1p(self) -> Self { +// $libm::ln_1p(self) +// } +// + #[inline] + fn sinh(self) -> Self { + $libm::sinh(self) + } + + #[inline] + fn cosh(self) -> Self { + $libm::cosh(self) + } + + #[inline] + fn tanh(self) -> Self { + $libm::tanh(self) + } + + #[inline] + fn asinh(self) -> Self { + $libm::asinh(self) + } + + #[inline] + fn acosh(self) -> Self { + $libm::acosh(self) + } + + #[inline] + fn atanh(self) -> Self { + $libm::atanh(self) + } + + #[inline] + fn is_finite(&self) -> bool { + $M::is_finite(*self) + } + } + )*) +); + +#[cfg(not(feature = "std"))] +impl_complex!( + f32, f32, F32Ext; + f64, f64, F64Ext +); + +#[cfg(feature = "std")] +impl_complex!( + f32,f32,f32; + f64,f64,f64 +); + +//#[cfg(feature = "decimal")] +//impl_real!(d128, d128, d128); + + +impl ComplexField for num_complex::Complex { + type RealField = N; + + #[inline] + fn from_real(re: Self::RealField) -> Self { + Self::new(re, Self::RealField::zero()) + } + + #[inline] + fn real(self) -> Self::RealField { + self.re + } + + #[inline] + fn imaginary(self) -> Self::RealField { + self.im + } + + #[inline] + fn argument(self) -> Self::RealField { + self.im.atan2(self.re) + } + + #[inline] + fn modulus(self) -> Self::RealField { + self.re.hypot(self.im) + } + + #[inline] + fn modulus_squared(self) -> Self::RealField { + self.re * self.re + self.im * self.im + } + + #[inline] + fn norm1(self) -> Self::RealField { + self.re.abs() + self.im.abs() + } + + #[inline] + fn recip(self) -> Self { + Self::one() / self + } + + #[inline] + fn conjugate(self) -> Self { + self.conj() + } + + #[inline] + fn scale(self, factor: Self::RealField) -> Self { + self * factor + } + + #[inline] + fn unscale(self, factor: Self::RealField) -> Self { + self / factor + } + + #[inline] + fn floor(self) -> Self { + Self::new(self.re.floor(), self.im.floor()) + } + + #[inline] + fn ceil(self) -> Self { + Self::new(self.re.ceil(), self.im.ceil()) + } + + #[inline] + fn round(self) -> Self { + Self::new(self.re.round(), self.im.round()) + } + + #[inline] + fn trunc(self) -> Self { + Self::new(self.re.trunc(), self.im.trunc()) + } + + #[inline] + fn fract(self) -> Self { + Self::new(self.re.fract(), self.im.fract()) + } + + #[inline] + fn mul_add(self, a: Self, b: Self) -> Self { + self * a + b + } + + + #[inline] + fn abs(self) -> Self::RealField { + self.modulus() + } + + #[inline] + fn exp2(self) -> Self { + let _2 = N::one() + N::one(); + num_complex::Complex::new(_2, N::zero()).powc(self) + } + + + #[inline] + fn exp_m1(self) -> Self { + self.exp() - Self::one() + } + + #[inline] + fn ln_1p(self) -> Self { + (Self::one() + self).ln() + } + + #[inline] + fn log2(self) -> Self { + let _2 = N::one() + N::one(); + self.log(_2) + } + + #[inline] + fn log10(self) -> Self { + let _10 = N::from_subset(&10.0f64); + self.log(_10) + } + + #[inline] + fn cbrt(self) -> Self { + let one_third = N::from_subset(&(1.0 / 3.0)); + self.powf(one_third) + } + + #[inline] + fn powi(self, n: i32) -> Self { + // FIXME: is there a more accurate solution? + let n = N::from_subset(&(n as f64)); + self.powf(n) + } + + #[inline] + fn is_finite(&self) -> bool { + self.re.is_finite() && self.im.is_finite() + } + + /* + * + * + * Unfortunately we are forced to copy-paste all + * those impls from https://github.com/rust-num/num-complex/blob/master/src/lib.rs + * to avoid requiring `std`. + * + * + */ + /// Computes `e^(self)`, where `e` is the base of the natural logarithm. + #[inline] + fn exp(self) -> Self { + // formula: e^(a + bi) = e^a (cos(b) + i*sin(b)) + // = from_polar(e^a, b) + complex_from_polar(self.re.exp(), self.im) + } + + /// Computes the principal value of natural logarithm of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 0]`, continuous from above. + /// + /// The branch satisfies `-π ≤ arg(ln(z)) ≤ π`. + #[inline] + fn ln(self) -> Self { + // formula: ln(z) = ln|z| + i*arg(z) + let (r, theta) = self.to_polar(); + Self::new(r.ln(), theta) + } + + /// Computes the principal value of the square root of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 0)`, continuous from above. + /// + /// The branch satisfies `-π/2 ≤ arg(sqrt(z)) ≤ π/2`. + #[inline] + fn sqrt(self) -> Self { + // formula: sqrt(r e^(it)) = sqrt(r) e^(it/2) + let two = N::one() + N::one(); + let (r, theta) = self.to_polar(); + complex_from_polar(r.sqrt(), theta / two) + } + + #[inline] + fn try_sqrt(self) -> Option { + Some(self.sqrt()) + } + + #[inline] + fn hypot(self, b: Self) -> Self::RealField { + (self.modulus_squared() + b.modulus_squared()).sqrt() + } + + /// Raises `self` to a floating point power. + #[inline] + fn powf(self, exp: Self::RealField) -> Self { + // formula: x^y = (ρ e^(i θ))^y = ρ^y e^(i θ y) + // = from_polar(ρ^y, θ y) + let (r, theta) = self.to_polar(); + complex_from_polar(r.powf(exp), theta * exp) + } + + /// Returns the logarithm of `self` with respect to an arbitrary base. + #[inline] + fn log(self, base: N) -> Self { + // formula: log_y(x) = log_y(ρ e^(i θ)) + // = log_y(ρ) + log_y(e^(i θ)) = log_y(ρ) + ln(e^(i θ)) / ln(y) + // = log_y(ρ) + i θ / ln(y) + let (r, theta) = self.to_polar(); + Self::new(r.log(base), theta / base.ln()) + } + + /// Raises `self` to a complex power. + #[inline] + fn powc(self, exp: Self) -> Self { + // formula: x^y = (a + i b)^(c + i d) + // = (ρ e^(i θ))^c (ρ e^(i θ))^(i d) + // where ρ=|x| and θ=arg(x) + // = ρ^c e^(−d θ) e^(i c θ) ρ^(i d) + // = p^c e^(−d θ) (cos(c θ) + // + i sin(c θ)) (cos(d ln(ρ)) + i sin(d ln(ρ))) + // = p^c e^(−d θ) ( + // cos(c θ) cos(d ln(ρ)) − sin(c θ) sin(d ln(ρ)) + // + i(cos(c θ) sin(d ln(ρ)) + sin(c θ) cos(d ln(ρ)))) + // = p^c e^(−d θ) (cos(c θ + d ln(ρ)) + i sin(c θ + d ln(ρ))) + // = from_polar(p^c e^(−d θ), c θ + d ln(ρ)) + let (r, theta) = self.to_polar(); + complex_from_polar( + r.powf(exp.re) * (-exp.im * theta).exp(), + exp.re * theta + exp.im * r.ln(), + ) + } + +/* + /// Raises a floating point number to the complex power `self`. + #[inline] + fn expf(&self, base: T) -> Self { + // formula: x^(a+bi) = x^a x^bi = x^a e^(b ln(x) i) + // = from_polar(x^a, b ln(x)) + Self::from_polar(&base.powf(self.re), &(self.im * base.ln())) + } + */ + + /// Computes the sine of `self`. + #[inline] + fn sin(self) -> Self { + // formula: sin(a + bi) = sin(a)cosh(b) + i*cos(a)sinh(b) + Self::new( + self.re.sin() * self.im.cosh(), + self.re.cos() * self.im.sinh(), + ) + } + + /// Computes the cosine of `self`. + #[inline] + fn cos(self) -> Self { + // formula: cos(a + bi) = cos(a)cosh(b) - i*sin(a)sinh(b) + Self::new( + self.re.cos() * self.im.cosh(), + -self.re.sin() * self.im.sinh(), + ) + } + + #[inline] + fn sin_cos(self) -> (Self, Self) { + let (rsin, rcos) = self.re.sin_cos(); + let (isinh, icosh) = self.im.sinh_cosh(); + let sin = Self::new( + rsin * icosh, + rcos * isinh, + ); + let cos = Self::new( + rcos * icosh, + -rsin * isinh, + ); + + (sin, cos) + } + + /// Computes the tangent of `self`. + #[inline] + fn tan(self) -> Self { + // formula: tan(a + bi) = (sin(2a) + i*sinh(2b))/(cos(2a) + cosh(2b)) + let (two_re, two_im) = (self.re + self.re, self.im + self.im); + Self::new(two_re.sin(), two_im.sinh()).unscale(two_re.cos() + two_im.cosh()) + } + + /// Computes the principal value of the inverse sine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1)`, continuous from above. + /// * `(1, ∞)`, continuous from below. + /// + /// The branch satisfies `-π/2 ≤ Re(asin(z)) ≤ π/2`. + #[inline] + fn asin(self) -> Self { + // formula: arcsin(z) = -i ln(sqrt(1-z^2) + iz) + let i = Self::i(); + -i * ((Self::one() - self * self).sqrt() + i * self).ln() + } + + /// Computes the principal value of the inverse cosine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1)`, continuous from above. + /// * `(1, ∞)`, continuous from below. + /// + /// The branch satisfies `0 ≤ Re(acos(z)) ≤ π`. + #[inline] + fn acos(self) -> Self { + // formula: arccos(z) = -i ln(i sqrt(1-z^2) + z) + let i = Self::i(); + -i * (i * (Self::one() - self * self).sqrt() + self).ln() + } + + /// Computes the principal value of the inverse tangent of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞i, -i]`, continuous from the left. + /// * `[i, ∞i)`, continuous from the right. + /// + /// The branch satisfies `-π/2 ≤ Re(atan(z)) ≤ π/2`. + #[inline] + fn atan(self) -> Self { + // formula: arctan(z) = (ln(1+iz) - ln(1-iz))/(2i) + let i = Self::i(); + let one = Self::one(); + let two = one + one; + + if self == i { + return Self::new(N::zero(), N::one() / N::zero()); + } else if self == -i { + return Self::new(N::zero(), -N::one() / N::zero()); + } + + ((one + i * self).ln() - (one - i * self).ln()) / (two * i) + } + + /// Computes the hyperbolic sine of `self`. + #[inline] + fn sinh(self) -> Self { + // formula: sinh(a + bi) = sinh(a)cos(b) + i*cosh(a)sin(b) + Self::new( + self.re.sinh() * self.im.cos(), + self.re.cosh() * self.im.sin(), + ) + } + + /// Computes the hyperbolic cosine of `self`. + #[inline] + fn cosh(self) -> Self { + // formula: cosh(a + bi) = cosh(a)cos(b) + i*sinh(a)sin(b) + Self::new( + self.re.cosh() * self.im.cos(), + self.re.sinh() * self.im.sin(), + ) + } + + #[inline] + fn sinh_cosh(self) -> (Self, Self) { + let (rsinh, rcosh) = self.re.sinh_cosh(); + let (isin, icos) = self.im.sin_cos(); + let sin = Self::new( + rsinh * icos, + rcosh * isin, + ); + let cos = Self::new( + rcosh * icos, + rsinh * isin, + ); + + (sin, cos) + } + + /// Computes the hyperbolic tangent of `self`. + #[inline] + fn tanh(self) -> Self { + // formula: tanh(a + bi) = (sinh(2a) + i*sin(2b))/(cosh(2a) + cos(2b)) + let (two_re, two_im) = (self.re + self.re, self.im + self.im); + Self::new(two_re.sinh(), two_im.sin()).unscale(two_re.cosh() + two_im.cos()) + } + + /// Computes the principal value of inverse hyperbolic sine of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞i, -i)`, continuous from the left. + /// * `(i, ∞i)`, continuous from the right. + /// + /// The branch satisfies `-π/2 ≤ Im(asinh(z)) ≤ π/2`. + #[inline] + fn asinh(self) -> Self { + // formula: arcsinh(z) = ln(z + sqrt(1+z^2)) + let one = Self::one(); + (self + (one + self * self).sqrt()).ln() + } + + /// Computes the principal value of inverse hyperbolic cosine of `self`. + /// + /// This function has one branch cut: + /// + /// * `(-∞, 1)`, continuous from above. + /// + /// The branch satisfies `-π ≤ Im(acosh(z)) ≤ π` and `0 ≤ Re(acosh(z)) < ∞`. + #[inline] + fn acosh(self) -> Self { + // formula: arccosh(z) = 2 ln(sqrt((z+1)/2) + sqrt((z-1)/2)) + let one = Self::one(); + let two = one + one; + two * (((self + one) / two).sqrt() + ((self - one) / two).sqrt()).ln() + } + + /// Computes the principal value of inverse hyperbolic tangent of `self`. + /// + /// This function has two branch cuts: + /// + /// * `(-∞, -1]`, continuous from above. + /// * `[1, ∞)`, continuous from below. + /// + /// The branch satisfies `-π/2 ≤ Im(atanh(z)) ≤ π/2`. + #[inline] + fn atanh(self) -> Self { + // formula: arctanh(z) = (ln(1+z) - ln(1-z))/2 + let one = Self::one(); + let two = one + one; + if self == one { + return Self::new(N::one() / N::zero(), N::zero()); + } else if self == -one { + return Self::new(-N::one() / N::zero(), N::zero()); + } + ((one + self).ln() - (one - self).ln()) / two + } +} + +#[inline] +fn complex_from_polar(r: N, theta: N) -> num_complex::Complex { + num_complex::Complex::new(r * theta.cos(), r * theta.sin()) +} \ No newline at end of file diff --git a/alga/src/general/identity.rs b/alga/src/general/identity.rs index b78810d..2fb630d 100644 --- a/alga/src/general/identity.rs +++ b/alga/src/general/identity.rs @@ -11,7 +11,7 @@ use num_complex::Complex; use approx::{AbsDiffEq, RelativeEq, UlpsEq}; -use general::{ +use crate::general::{ AbstractGroup, AbstractGroupAbelian, AbstractLoop, AbstractMagma, AbstractMonoid, AbstractQuasigroup, AbstractSemigroup, Additive, TwoSidedInverse, JoinSemilattice, Lattice, MeetSemilattice, Multiplicative, Operator, SubsetOf, diff --git a/alga/src/general/lattice.rs b/alga/src/general/lattice.rs index 40a78da..109bf48 100644 --- a/alga/src/general/lattice.rs +++ b/alga/src/general/lattice.rs @@ -122,3 +122,24 @@ macro_rules! impl_lattice( impl_lattice!(u8, u16, u32, u64, usize, i8, i16, i32, i64, isize, f32, f64); #[cfg(feature = "decimal")] impl_lattice!(d128); + + +impl MeetSemilattice for num_complex::Complex { + #[inline] + fn meet(&self, other: &Self) -> Self { + Self { + re: self.re.meet(&other.re), + im: self.im.meet(&other.im) + } + } +} + +impl JoinSemilattice for num_complex::Complex { + #[inline] + fn join(&self, other: &Self) -> Self { + Self { + re: self.re.join(&other.re), + im: self.im.join(&other.im) + } + } +} \ No newline at end of file diff --git a/alga/src/general/mod.rs b/alga/src/general/mod.rs index 975ada0..ba4683f 100644 --- a/alga/src/general/mod.rs +++ b/alga/src/general/mod.rs @@ -14,7 +14,7 @@ //! Fundamental algebraic structures. //! -//! For most applications requiring an abstraction over the reals, `Real` +//! For most applications requiring an abstraction over the reals, `RealField` //! should be sufficient. //! //! ## Algebraic properties @@ -171,7 +171,8 @@ pub use self::specialized::{AdditiveGroup, AdditiveGroupAbelian, AdditiveLoop, A MultiplicativeGroup, MultiplicativeGroupAbelian, MultiplicativeLoop, MultiplicativeMagma, MultiplicativeMonoid, MultiplicativeQuasigroup, MultiplicativeSemigroup, Ring, RingCommutative}; -pub use self::real::Real; +pub use self::real::RealField; +pub use self::complex::ComplexField; #[macro_use] mod one_operator; @@ -180,8 +181,15 @@ mod module; mod identity; mod operator; mod real; +mod complex; mod lattice; mod subset; mod specialized; #[doc(hidden)] pub mod wrapper; + +#[deprecated(note = "This has been renamed `RealField`.")] +/// The field of reals. This has been renamed to `RealField`. +pub trait Real: RealField {} + +impl Real for T {} \ No newline at end of file diff --git a/alga/src/general/module.rs b/alga/src/general/module.rs index 8c78599..0a89c4d 100644 --- a/alga/src/general/module.rs +++ b/alga/src/general/module.rs @@ -1,4 +1,4 @@ -use general::{AbstractGroupAbelian, AbstractRingCommutative, Additive, Multiplicative, Operator}; +use crate::general::{AbstractGroupAbelian, AbstractRingCommutative, Additive, Multiplicative, Operator}; /// A module combines two sets: one with an Abelian group structure and another with a /// commutative ring structure. diff --git a/alga/src/general/one_operator.rs b/alga/src/general/one_operator.rs index b93a253..3871f2b 100644 --- a/alga/src/general/one_operator.rs +++ b/alga/src/general/one_operator.rs @@ -6,7 +6,7 @@ use decimal::d128; use approx::RelativeEq; -use general::{Additive, ClosedNeg, Identity, TwoSidedInverse, Multiplicative, Operator}; +use crate::general::{Additive, ClosedNeg, Identity, TwoSidedInverse, Multiplicative, Operator}; /// A magma is an algebraic structure which consists of a set equipped with a binary operation, ∘, /// which must be closed. diff --git a/alga/src/general/real.rs b/alga/src/general/real.rs index 5e094ff..d183577 100644 --- a/alga/src/general/real.rs +++ b/alga/src/general/real.rs @@ -1,12 +1,9 @@ -use num::{Bounded, FromPrimitive, Num, NumAssign, Signed, Zero, One}; -use std::any::Any; -use std::fmt::{Debug, Display}; -use std::ops::Neg; +use num::{Bounded, Signed}; use std::{f32, f64}; use approx::{RelativeEq, UlpsEq}; -use general::{Field, Lattice, SubsetOf, SupersetOf}; +use crate::general::{ComplexField, Lattice}; #[cfg(not(feature = "std"))] use libm::F32Ext; @@ -26,68 +23,21 @@ use num; // FIXME: SubsetOf should be removed when specialization will be supported by rustc. This will // allow a blanket impl: impl SubsetOf for T { ... } // NOTE: make all types debuggable/'static/Any ? This seems essential for any kind of generic programming. -pub trait Real: - SubsetOf - + SupersetOf - + Field - + Copy - + Num - + NumAssign - + FromPrimitive - + Neg +pub trait RealField: + ComplexField + RelativeEq + UlpsEq + Lattice + Signed - + Send - + Sync - + Any - + 'static - + Debug - + Display + Bounded { // NOTE: a real must be bounded because, no matter the chosen representation, being `Copy` implies that it occupies a statically-known size, meaning that it must have min/max values. - fn floor(self) -> Self; - fn ceil(self) -> Self; - fn round(self) -> Self; - fn trunc(self) -> Self; - fn fract(self) -> Self; - fn abs(self) -> Self; - fn signum(self) -> Self; + fn is_sign_positive(self) -> bool; fn is_sign_negative(self) -> bool; - fn mul_add(self, a: Self, b: Self) -> Self; - fn recip(self) -> Self; - fn powi(self, n: i32) -> Self; - fn powf(self, n: Self) -> Self; - fn sqrt(self) -> Self; - fn exp(self) -> Self; - fn exp2(self) -> Self; - fn ln(self) -> Self; - fn log(self, base: Self) -> Self; - fn log2(self) -> Self; - fn log10(self) -> Self; fn max(self, other: Self) -> Self; fn min(self, other: Self) -> Self; - fn cbrt(self) -> Self; - fn hypot(self, other: Self) -> Self; - fn sin(self) -> Self; - fn cos(self) -> Self; - fn tan(self) -> Self; - fn asin(self) -> Self; - fn acos(self) -> Self; - fn atan(self) -> Self; fn atan2(self, other: Self) -> Self; - fn sin_cos(self) -> (Self, Self); - fn exp_m1(self) -> Self; - fn ln_1p(self) -> Self; - fn sinh(self) -> Self; - fn cosh(self) -> Self; - fn tanh(self) -> Self; - fn asinh(self) -> Self; - fn acosh(self) -> Self; - fn atanh(self) -> Self; fn pi() -> Self; fn two_pi() -> Self; @@ -105,69 +55,11 @@ pub trait Real: fn log10_e() -> Self; fn ln_2() -> Self; fn ln_10() -> Self; - - fn is_finite(&self) -> bool; - - /// Cardinal sine - #[inline] - fn sinc(self) -> Self { - if self == Self::zero() { - Self::one() - } - else { - self.sin() / self - } - } - - #[inline] - fn sinhc(self) -> Self { - if self == Self::zero() { - Self::one() - } - else { - self.sinh() / self - } - } } macro_rules! impl_real( ($($T:ty, $M:ident, $libm: ident);*) => ($( - impl Real for $T { - #[inline] - fn floor(self) -> Self { - $libm::floor(self) - } - - #[inline] - fn ceil(self) -> Self { - $libm::ceil(self) - } - - #[inline] - fn round(self) -> Self { - $libm::round(self) - } - - #[inline] - fn trunc(self) -> Self { - $libm::trunc(self) - } - - #[inline] - fn fract(self) -> Self { - $libm::fract(self) - } - - #[inline] - fn abs(self) -> Self { - $libm::abs(self) - } - - #[inline] - fn signum(self) -> Self { - Signed::signum(&self) - } - + impl RealField for $T { #[inline] fn is_sign_positive(self) -> bool { $M::is_sign_positive(self) @@ -178,69 +70,6 @@ macro_rules! impl_real( $M::is_sign_negative(self) } - #[inline] - fn mul_add(self, a: Self, b: Self) -> Self { - $libm::mul_add(self, a, b) - } - - #[inline] - fn recip(self) -> Self { - $M::recip(self) - } - - #[cfg(feature = "std")] - #[inline] - fn powi(self, n: i32) -> Self { - self.powi(n) - } - - #[cfg(not(feature = "std"))] - #[inline] - fn powi(self, n: i32) -> Self { - // FIXME: is there a more efficient solution? - num::pow(self, n as usize) - } - - #[inline] - fn powf(self, n: Self) -> Self { - $libm::powf(self, n) - } - - #[inline] - fn sqrt(self) -> Self { - $libm::sqrt(self) - } - - #[inline] - fn exp(self) -> Self { - $libm::exp(self) - } - - #[inline] - fn exp2(self) -> Self { - $libm::exp2(self) - } - - #[inline] - fn ln(self) -> Self { - $libm::ln(self) - } - - #[inline] - fn log(self, base: Self) -> Self { - $libm::log(self, base) - } - - #[inline] - fn log2(self) -> Self { - $libm::log2(self) - } - - #[inline] - fn log10(self) -> Self { - $libm::log10(self) - } - #[inline] fn max(self, other: Self) -> Self { $M::max(self, other) @@ -251,96 +80,11 @@ macro_rules! impl_real( $M::min(self, other) } - #[inline] - fn cbrt(self) -> Self { - $libm::cbrt(self) - } - - #[inline] - fn hypot(self, other: Self) -> Self { - $libm::hypot(self, other) - } - - #[inline] - fn sin(self) -> Self { - $libm::sin(self) - } - - #[inline] - fn cos(self) -> Self { - $libm::cos(self) - } - - #[inline] - fn tan(self) -> Self { - $libm::tan(self) - } - - #[inline] - fn asin(self) -> Self { - $libm::asin(self) - } - - #[inline] - fn acos(self) -> Self { - $libm::acos(self) - } - - #[inline] - fn atan(self) -> Self { - $libm::atan(self) - } - #[inline] fn atan2(self, other: Self) -> Self { $libm::atan2(self, other) } - #[inline] - fn sin_cos(self) -> (Self, Self) { - $libm::sin_cos(self) - } - - #[inline] - fn exp_m1(self) -> Self { - $libm::exp_m1(self) - } - - #[inline] - fn ln_1p(self) -> Self { - $libm::ln_1p(self) - } - - #[inline] - fn sinh(self) -> Self { - $libm::sinh(self) - } - - #[inline] - fn cosh(self) -> Self { - $libm::cosh(self) - } - - #[inline] - fn tanh(self) -> Self { - $libm::tanh(self) - } - - #[inline] - fn asinh(self) -> Self { - $libm::asinh(self) - } - - #[inline] - fn acosh(self) -> Self { - $libm::acosh(self) - } - - #[inline] - fn atanh(self) -> Self { - $libm::atanh(self) - } - /// Archimedes' constant. #[inline] fn pi() -> Self { @@ -431,10 +175,6 @@ macro_rules! impl_real( fn ln_10() -> Self { $M::consts::LN_10 } - - fn is_finite(&self) -> bool { - $M::is_finite(*self) - } } )*) ); diff --git a/alga/src/general/specialized.rs b/alga/src/general/specialized.rs index 756e581..ad242db 100644 --- a/alga/src/general/specialized.rs +++ b/alga/src/general/specialized.rs @@ -1,5 +1,5 @@ use num::{One, Zero}; -use general::{AbstractField, AbstractGroup, AbstractGroupAbelian, AbstractLoop, AbstractMagma, +use crate::general::{AbstractField, AbstractGroup, AbstractGroupAbelian, AbstractLoop, AbstractMagma, AbstractModule, AbstractMonoid, AbstractQuasigroup, AbstractRing, AbstractRingCommutative, AbstractSemigroup, Additive, ClosedAdd, ClosedDiv, ClosedMul, ClosedNeg, ClosedSub, Multiplicative}; diff --git a/alga/src/general/two_operators.rs b/alga/src/general/two_operators.rs index 3e792e9..6738204 100644 --- a/alga/src/general/two_operators.rs +++ b/alga/src/general/two_operators.rs @@ -4,8 +4,8 @@ use num_complex::Complex; #[cfg(feature = "decimal")] use decimal::d128; -use general::wrapper::Wrapper as W; -use general::{ +use crate::general::wrapper::Wrapper as W; +use crate::general::{ AbstractGroupAbelian, AbstractMonoid, Additive, ClosedNeg, Multiplicative, Operator, }; diff --git a/alga/src/general/wrapper.rs b/alga/src/general/wrapper.rs index c4bae35..f945e47 100644 --- a/alga/src/general/wrapper.rs +++ b/alga/src/general/wrapper.rs @@ -7,9 +7,9 @@ use std::ops::{Add, Div, Mul, Neg, Sub}; use approx::{AbsDiffEq, RelativeEq, UlpsEq}; -use general::AbstractMagma; -use general::AbstractQuasigroup; -use general::{TwoSidedInverse, Operator}; +use crate::general::AbstractMagma; +use crate::general::AbstractQuasigroup; +use crate::general::{TwoSidedInverse, Operator}; /// Wrapper that allows to use operators on algebraic types. #[derive(Debug)] diff --git a/alga/src/lib.rs b/alga/src/lib.rs index 820432d..ffb657b 100644 --- a/alga/src/lib.rs +++ b/alga/src/lib.rs @@ -26,13 +26,10 @@ extern crate approx; #[cfg(feature = "decimal")] #[macro_use] extern crate decimal; -#[cfg(not(feature = "std"))] -extern crate libm; -extern crate num_complex; extern crate num_traits as num; #[cfg(not(feature = "std"))] -use core as std; +extern crate core as std; #[macro_use] mod macros; diff --git a/alga/src/linear/id.rs b/alga/src/linear/id.rs index 565e21c..041348c 100644 --- a/alga/src/linear/id.rs +++ b/alga/src/linear/id.rs @@ -1,7 +1,7 @@ use num; -use general::{Id, Identity}; -use linear::{AffineTransformation, DirectIsometry, EuclideanSpace, InnerSpace, Isometry, +use crate::general::{Id, Identity}; +use crate::linear::{AffineTransformation, DirectIsometry, EuclideanSpace, InnerSpace, Isometry, OrthogonalTransformation, ProjectiveTransformation, Rotation, Scaling, Similarity, Transformation, Translation}; @@ -99,7 +99,7 @@ impl OrthogonalTransformation for Id {} impl Rotation for Id { #[inline] - fn powf(&self, _: E::Real) -> Option { + fn powf(&self, _: E::RealField) -> Option { Some(Id::new()) } @@ -113,7 +113,7 @@ impl Rotation for Id { } #[inline] - fn scaled_rotation_between(a: &E::Coordinates, b: &E::Coordinates, _: E::Real) -> Option { + fn scaled_rotation_between(a: &E::Coordinates, b: &E::Coordinates, _: E::RealField) -> Option { Rotation::::rotation_between(a, b) } } diff --git a/alga/src/linear/matrix.rs b/alga/src/linear/matrix.rs index 3c8bf65..cd50c3a 100644 --- a/alga/src/linear/matrix.rs +++ b/alga/src/linear/matrix.rs @@ -1,7 +1,7 @@ use std::ops::Mul; -use general::{Field, MultiplicativeGroup, MultiplicativeMonoid}; -use linear::FiniteDimVectorSpace; +use crate::general::{Field, MultiplicativeGroup, MultiplicativeMonoid}; +use crate::linear::FiniteDimVectorSpace; /// The space of all matrices. pub trait Matrix: diff --git a/alga/src/linear/transformation.rs b/alga/src/linear/transformation.rs index 1d016ee..01701b9 100644 --- a/alga/src/linear/transformation.rs +++ b/alga/src/linear/transformation.rs @@ -1,6 +1,6 @@ -use general::{ClosedDiv, ClosedMul, ClosedNeg, Id, TwoSidedInverse, MultiplicativeGroup, - MultiplicativeMonoid, Real, SubsetOf}; -use linear::{EuclideanSpace, NormedSpace}; +use crate::general::{ClosedDiv, ClosedMul, ClosedNeg, Id, TwoSidedInverse, MultiplicativeGroup, + MultiplicativeMonoid, RealField, SubsetOf, ComplexField}; +use crate::linear::{EuclideanSpace, NormedSpace}; // NOTE: A subgroup trait inherit from its parent groups. @@ -16,7 +16,7 @@ pub trait Transformation: MultiplicativeMonoid { fn transform_vector(&self, v: &E::Coordinates) -> E::Coordinates; } -/// The most general form of inversible transformations on an euclidean space. +/// The most general form of invertible transformations on an euclidean space. pub trait ProjectiveTransformation : MultiplicativeGroup + Transformation { /// Applies this group's two_sided_inverse action on a point from the euclidean space. @@ -194,24 +194,24 @@ pub trait OrthogonalTransformation /// Subgroups of the (signed) uniform scaling group. pub trait Scaling : AffineTransformation - + SubsetOf { + + SubsetOf { /// Converts this scaling factor to a real. Same as `self.to_superset()`. #[inline] - fn to_real(&self) -> E::Real { + fn to_real(&self) -> E::RealField { self.to_superset() } /// Attempts to convert a real to an element of this scaling subgroup. Same as /// `Self::from_superset()`. Returns `None` if no such scaling is possible for this subgroup. #[inline] - fn from_real(r: E::Real) -> Option { + fn from_real(r: E::RealField) -> Option { Self::from_superset(&r) } /// Raises the scaling to a power. The result must be equivalent to /// `self.to_superset().powf(n)`. Returns `None` if the result is not representable by `Self`. #[inline] - fn powf(&self, n: E::Real) -> Option { + fn powf(&self, n: E::RealField) -> Option { Self::from_superset(&self.to_superset().powf(n)) } @@ -240,7 +240,7 @@ pub trait Translation /// Raises the translation to a power. The result must be equivalent to /// `self.to_superset() * n`. Returns `None` if the result is not representable by `Self`. #[inline] - fn powf(&self, n: E::Real) -> Option { + fn powf(&self, n: E::RealField) -> Option { Self::from_vector(self.to_vector() * n) } @@ -257,7 +257,7 @@ pub trait Rotation { /// Raises this rotation to a power. If this is a simple rotation, the result must be /// equivalent to multiplying the rotation angle by `n`. - fn powf(&self, n: E::Real) -> Option; + fn powf(&self, n: E::RealField) -> Option; /// Computes a simple rotation that makes the angle between `a` and `b` equal to zero, i.e., /// `b.angle(a * delta_rotation(a, b)) = 0`. If `a` and `b` are collinear, the computed @@ -270,7 +270,7 @@ pub trait Rotation /// This is equivalent to calling `self.rotation_between(a, b)` followed by `.powf(n)` but will /// usually be much more efficient. #[inline] - fn scaled_rotation_between(a: &E::Coordinates, b: &E::Coordinates, s: E::Real) -> Option; + fn scaled_rotation_between(a: &E::Coordinates, b: &E::Coordinates, s: E::RealField) -> Option; // FIXME: add a function that computes the rotation with the axis orthogonal to Span(a, b) and // with angle equal to `n`? @@ -284,8 +284,8 @@ pub trait Rotation impl Transformation for R where - R: Real, - E: EuclideanSpace, + R: RealField, + E: EuclideanSpace, E::Coordinates: ClosedMul + ClosedDiv + ClosedNeg, { #[inline] @@ -301,8 +301,8 @@ where impl ProjectiveTransformation for R where - R: Real, - E: EuclideanSpace, + R: RealField, + E: EuclideanSpace, E::Coordinates: ClosedMul + ClosedDiv + ClosedNeg, { #[inline] @@ -320,8 +320,8 @@ where impl AffineTransformation for R where - R: Real, - E: EuclideanSpace, + R: RealField, + E: EuclideanSpace, E::Coordinates: ClosedMul + ClosedDiv + ClosedNeg, { type Rotation = Id; @@ -366,22 +366,22 @@ where impl Scaling for R where - R: Real + SubsetOf, - E: EuclideanSpace, + R: RealField + SubsetOf, + E: EuclideanSpace, E::Coordinates: ClosedMul + ClosedDiv + ClosedNeg, { #[inline] - fn to_real(&self) -> E::Real { + fn to_real(&self) -> E::RealField { *self } #[inline] - fn from_real(r: E::Real) -> Option { + fn from_real(r: E::RealField) -> Option { Some(r) } #[inline] - fn powf(&self, n: E::Real) -> Option { + fn powf(&self, n: E::RealField) -> Option { Some(n.powf(n)) } @@ -393,8 +393,8 @@ where impl Similarity for R where - R: Real + SubsetOf, - E: EuclideanSpace, + R: RealField + SubsetOf, + E: EuclideanSpace, E::Coordinates: ClosedMul + ClosedDiv + ClosedNeg, { type Scaling = R; diff --git a/alga/src/linear/vector.rs b/alga/src/linear/vector.rs index efd5411..a3639ed 100644 --- a/alga/src/linear/vector.rs +++ b/alga/src/linear/vector.rs @@ -5,7 +5,7 @@ use std::ops::{ Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign, }; -use general::{ClosedAdd, ClosedDiv, ClosedMul, Field, Module, Real}; +use crate::general::{ClosedAdd, ClosedDiv, ClosedMul, Field, Module, RealField, ComplexField}; /// A vector space has a module structure over a field instead of a ring. pub trait VectorSpace: Module::Field> @@ -17,42 +17,44 @@ pub trait VectorSpace: Module::Field> } /// A normed vector space. -pub trait NormedSpace: VectorSpace { +pub trait NormedSpace: VectorSpace::ComplexField> { + /// The result of the norm (not necessarily the same same as the field used by this vector space). + type RealField: RealField; + /// The field of this space must be this complex number. + type ComplexField: ComplexField; + /// The squared norm of this vector. - fn norm_squared(&self) -> Self::Field; + fn norm_squared(&self) -> Self::RealField; /// The norm of this vector. - fn norm(&self) -> Self::Field; + fn norm(&self) -> Self::RealField; /// Returns a normalized version of this vector. fn normalize(&self) -> Self; /// Normalizes this vector in-place and returns its norm. - fn normalize_mut(&mut self) -> Self::Field; + fn normalize_mut(&mut self) -> Self::RealField; /// Returns a normalized version of this vector unless its norm as smaller or equal to `eps`. - fn try_normalize(&self, eps: Self::Field) -> Option; + fn try_normalize(&self, eps: Self::RealField) -> Option; /// Normalizes this vector in-place or does nothing if its norm is smaller or equal to `eps`. /// /// If the normalization succeeded, returns the old normal of this vector. - fn try_normalize_mut(&mut self, eps: Self::Field) -> Option; + fn try_normalize_mut(&mut self, eps: Self::RealField) -> Option; } /// A vector space equipped with an inner product. /// /// It must be a normed space as well and the norm must agree with the inner product. /// The inner product must be symmetric, linear in its first argument, and positive definite. -pub trait InnerSpace: NormedSpace::Real> { - /// The result of inner product (same as the field used by this vector space). - type Real: Real; - +pub trait InnerSpace: NormedSpace { /// Computes the inner product of `self` with `other`. - fn inner_product(&self, other: &Self) -> Self::Real; + fn inner_product(&self, other: &Self) -> Self::ComplexField; /// Measures the angle between two vectors. #[inline] - fn angle(&self, other: &Self) -> Self::Real { + fn angle(&self, other: &Self) -> Self::RealField { let prod = self.inner_product(other); let n1 = self.norm(); let n2 = other.norm(); @@ -60,12 +62,12 @@ pub trait InnerSpace: NormedSpace::Real> { if n1 == num::zero() || n2 == num::zero() { num::zero() } else { - let cang = prod / (n1 * n2); + let cang = prod.real() * n1 * n2; if cang > num::one() { num::zero() - } else if cang < -num::one::() { - Self::Real::pi() + } else if cang < -num::one::() { + Self::RealField::pi() } else { cang.acos() } @@ -109,7 +111,7 @@ pub trait FiniteDimVectorSpace: /// A finite-dimensional vector space equipped with an inner product that must coincide /// with the dot product. pub trait FiniteDimInnerSpace: - InnerSpace + FiniteDimVectorSpace::Real> + InnerSpace + FiniteDimVectorSpace::ComplexField> { /// Orthonormalizes the given family of vectors. The largest free family of vectors is moved at /// the beginning of the array and its size is returned. Vectors at an indices larger or equal to @@ -153,13 +155,13 @@ pub trait AffineSpace: /// The finite-dimensional affine space based on the field of reals. pub trait EuclideanSpace: AffineSpace::Coordinates> + // Equivalent to `.scale_by`. - ClosedMul<::Real> + + ClosedMul<::RealField> + // Equivalent to `.scale_by`. - ClosedDiv<::Real> + + ClosedDiv<::RealField> + // Equivalent to `.scale_by(-1.0)`. Neg { /// The underlying finite vector space. - type Coordinates: FiniteDimInnerSpace + + type Coordinates: FiniteDimInnerSpace + // XXX: the following bounds should not be necessary but the compiler does not // seem to be able to find them (from supertraits of VectorSpace)… Also, it won't // find them even if we add ClosedMul instead of Mul and MulAssign separately… @@ -167,19 +169,19 @@ pub trait EuclideanSpace: AffineSpace::Co AddAssign + Sub + SubAssign + - Mul + - MulAssign + - Div + - DivAssign + + Mul + + MulAssign + + Div + + DivAssign + Neg; // XXX: we can't write the following =( : - // type Vector: FiniteDimInnerSpace + InnerSpace; - // The compiler won't recognize that VectorSpace::Field = Self::Real. + // type Vector: FiniteDimInnerSpace + InnerSpace; + // The compiler won't recognize that VectorSpace::Field = Self::RealField. // Though it will work if only one bound is used… looks like a compiler bug. /// The underlying reals. - type Real: Real; + type RealField: RealField; /// The preferred origin of this euclidean space. /// @@ -191,7 +193,7 @@ pub trait EuclideanSpace: AffineSpace::Co /// /// Same as self * s. #[inline] - fn scale_by(&self, s: Self::Real) -> Self { + fn scale_by(&self, s: Self::RealField) -> Self { Self::from_coordinates(self.coordinates() * s) } @@ -210,13 +212,13 @@ pub trait EuclideanSpace: AffineSpace::Co /// The distance between two points. #[inline] - fn distance_squared(&self, b: &Self) -> Self::Real { + fn distance_squared(&self, b: &Self) -> Self::RealField { self.subtract(b).norm_squared() } /// The distance between two points. #[inline] - fn distance(&self, b: &Self) -> Self::Real { + fn distance(&self, b: &Self) -> Self::RealField { self.subtract(b).norm() } } @@ -227,14 +229,17 @@ impl VectorSpace for Complex { } -impl NormedSpace for Complex { +impl NormedSpace for Complex { + type RealField = N; + type ComplexField = N; + #[inline] - fn norm_squared(&self) -> Self::Field { + fn norm_squared(&self) -> Self::RealField { self.norm_sqr() } #[inline] - fn norm(&self) -> Self::Field { + fn norm(&self) -> Self::RealField { self.norm_sqr().sqrt() } @@ -244,14 +249,14 @@ impl NormedSpace for Complex { } #[inline] - fn normalize_mut(&mut self) -> Self::Field { + fn normalize_mut(&mut self) -> Self::RealField { let norm = self.norm(); *self /= norm; norm } #[inline] - fn try_normalize(&self, eps: Self::Field) -> Option { + fn try_normalize(&self, eps: Self::RealField) -> Option { let norm = self.norm_sqr(); if norm > eps * eps { Some(*self / norm.sqrt()) @@ -261,7 +266,7 @@ impl NormedSpace for Complex { } #[inline] - fn try_normalize_mut(&mut self, eps: Self::Field) -> Option { + fn try_normalize_mut(&mut self, eps: Self::RealField) -> Option { let sq_norm = self.norm_sqr(); if sq_norm > eps * eps { let norm = sq_norm.sqrt(); @@ -273,15 +278,5 @@ impl NormedSpace for Complex { } } - -impl InnerSpace for Complex { - type Real = N; - - #[inline] - fn inner_product(&self, other: &Self) -> Self::Real { - self.re * other.re + self.im * other.im - } -} - // Note: we can't implement FiniteDimVectorSpace for Complex because // the `Complex` type does not implement Index. \ No newline at end of file