diff --git a/src/lib.rs b/src/lib.rs index 215c222..caa0518 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -38,13 +38,15 @@ use core::str::FromStr; #[cfg(feature = "std")] use std::error::Error; -use traits::{Inv, MulAdd, Num, One, Signed, Zero}; +use traits::{Inv, MulAdd, Num, One, Pow, Signed, Zero}; #[cfg(feature = "std")] use traits::float::Float; use traits::float::FloatCore; mod cast; +mod pow; + #[cfg(feature = "rand")] mod crand; #[cfg(feature = "rand")] @@ -121,6 +123,12 @@ impl<T: Clone + Num> Complex<T> { pub fn unscale(&self, t: T) -> Self { Self::new(self.re.clone() / t.clone(), self.im.clone() / t) } + + /// Raises `self` to an unsigned integer power. + #[inline] + pub fn powu(&self, exp: u32) -> Self { + Pow::pow(self, exp) + } } impl<T: Clone + Num + Neg<Output = T>> Complex<T> { @@ -139,6 +147,12 @@ impl<T: Clone + Num + Neg<Output = T>> Complex<T> { -self.im.clone() / norm_sqr, ) } + + /// Raises `self` to a signed integer power. + #[inline] + pub fn powi(&self, exp: i32) -> Self { + Pow::pow(self, exp) + } } impl<T: Clone + Signed> Complex<T> { @@ -1487,10 +1501,29 @@ mod test { assert_eq!(_4_2i.l1_norm(), 6.0); } + #[test] + fn test_pow() { + for c in all_consts.iter() { + assert_eq!(c.powi(0), _1_0i); + let mut pos = _1_0i; + let mut neg = _1_0i; + for i in 1i32..20 { + pos *= c; + assert_eq!(pos, c.powi(i)); + if c.is_zero() { + assert!(c.powi(-i).is_nan()); + } else { + neg /= c; + assert_eq!(neg, c.powi(-i)); + } + } + } + } + #[cfg(feature = "std")] mod float { use super::*; - use traits::Float; + use traits::{Float, Pow}; #[test] #[cfg_attr(target_arch = "x86", ignore)] @@ -1535,7 +1568,11 @@ mod test { fn close_to_tol(a: Complex64, b: Complex64, tol: f64) -> bool { // returns true if a and b are reasonably close - (a == b) || (a - b).norm() < tol + let close = (a == b) || (a - b).norm() < tol; + if !close { + println!("{:?} != {:?}", a, b); + } + close } #[test] @@ -1592,9 +1629,11 @@ mod test { #[test] fn test_powf() { - let c = Complex::new(2.0, -1.0); - let r = c.powf(3.5); - assert!(close_to_tol(r, Complex::new(-0.8684746, -16.695934), 1e-5)); + let c = Complex64::new(2.0, -1.0); + let expected = Complex64::new(-0.8684746, -16.695934); + assert!(close_to_tol(c.powf(3.5), expected, 1e-5)); + assert!(close_to_tol(Pow::pow(c, 3.5_f64), expected, 1e-5)); + assert!(close_to_tol(Pow::pow(c, 3.5_f32), expected, 1e-5)); } #[test] diff --git a/src/pow.rs b/src/pow.rs new file mode 100644 index 0000000..2f6b5ba --- /dev/null +++ b/src/pow.rs @@ -0,0 +1,187 @@ +use super::Complex; + +use core::ops::Neg; +#[cfg(feature = "std")] +use traits::Float; +use traits::{Num, One, Pow}; + +macro_rules! pow_impl { + ($U:ty, $S:ty) => { + impl<'a, T: Clone + Num> Pow<$U> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, mut exp: $U) -> Self::Output { + if exp == 0 { + return Complex::one(); + } + let mut base = self.clone(); + + while exp & 1 == 0 { + base = base.clone() * base; + exp >>= 1; + } + + if exp == 1 { + return base; + } + + let mut acc = base.clone(); + while exp > 1 { + exp >>= 1; + base = base.clone() * base; + if exp & 1 == 1 { + acc = acc * base.clone(); + } + } + acc + } + } + + impl<'a, 'b, T: Clone + Num> Pow<&'b $U> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: &$U) -> Self::Output { + self.pow(*exp) + } + } + + impl<'a, T: Clone + Num + Neg<Output = T>> Pow<$S> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: $S) -> Self::Output { + if exp < 0 { + Pow::pow(&self.inv(), exp.wrapping_neg() as $U) + } else { + Pow::pow(self, exp as $U) + } + } + } + + impl<'a, 'b, T: Clone + Num + Neg<Output = T>> Pow<&'b $S> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: &$S) -> Self::Output { + self.pow(*exp) + } + } + }; +} + +pow_impl!(u8, i8); +pow_impl!(u16, i16); +pow_impl!(u32, i32); +pow_impl!(u64, i64); +pow_impl!(usize, isize); +#[cfg(has_i128)] +pow_impl!(u128, i128); + +// Note: we can't add `impl<T: Float> Pow<T> for Complex<T>` because new blanket impls are a +// breaking change. Someone could already have their own `F` and `impl Pow<F> for Complex<F>` +// which would conflict. We can't even do this in a new semantic version, because we have to +// gate it on the "std" feature, and features can't add breaking changes either. + +macro_rules! powf_impl { + ($F:ty) => { + #[cfg(feature = "std")] + impl<'a, T: Float> Pow<$F> for &'a Complex<T> + where + $F: Into<T>, + { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: $F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(feature = "std")] + impl<'a, 'b, T: Float> Pow<&'b $F> for &'a Complex<T> + where + $F: Into<T>, + { + type Output = Complex<T>; + + #[inline] + fn pow(self, &exp: &$F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(feature = "std")] + impl<T: Float> Pow<$F> for Complex<T> + where + $F: Into<T>, + { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: $F) -> Self::Output { + self.powf(exp.into()) + } + } + + #[cfg(feature = "std")] + impl<'b, T: Float> Pow<&'b $F> for Complex<T> + where + $F: Into<T>, + { + type Output = Complex<T>; + + #[inline] + fn pow(self, &exp: &$F) -> Self::Output { + self.powf(exp.into()) + } + } + }; +} + +powf_impl!(f32); +powf_impl!(f64); + +// These blanket impls are OK, because both the target type and the trait parameter would be +// foreign to anyone else trying to implement something that would overlap, raising E0117. + +#[cfg(feature = "std")] +impl<'a, T: Float> Pow<Complex<T>> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: Complex<T>) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(feature = "std")] +impl<'a, 'b, T: Float> Pow<&'b Complex<T>> for &'a Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, &exp: &'b Complex<T>) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(feature = "std")] +impl<T: Float> Pow<Complex<T>> for Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, exp: Complex<T>) -> Self::Output { + self.powc(exp) + } +} + +#[cfg(feature = "std")] +impl<'b, T: Float> Pow<&'b Complex<T>> for Complex<T> { + type Output = Complex<T>; + + #[inline] + fn pow(self, &exp: &'b Complex<T>) -> Self::Output { + self.powc(exp) + } +}