Add Complex::powi and assorted Pow impls

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]

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)
+    }
+}