Browse Source

Add general Rem and Num implementations for Complex<T>

Isaac Carruthers 7 years ago
parent
commit
1b671ca43e
1 changed files with 227 additions and 32 deletions
  1. 227 32
      complex/src/lib.rs

+ 227 - 32
complex/src/lib.rs

@@ -26,7 +26,7 @@ use std::error::Error;
 use std::fmt;
 #[cfg(test)]
 use std::hash;
-use std::ops::{Add, Div, Mul, Neg, Sub};
+use std::ops::{Add, Div, Mul, Neg, Sub, Rem};
 use std::str::FromStr;
 
 use traits::{Zero, One, Num, Float};
@@ -261,8 +261,8 @@ impl<T: Clone + Float> Complex<T> {
     #[inline]
     pub fn asin(&self) -> Complex<T> {
         // formula: arcsin(z) = -i ln(sqrt(1-z^2) + iz)
-        let i = Complex::i();
-        -i*((Complex::one() - self*self).sqrt() + i*self).ln()
+        let i = Complex::<T>::i();
+        -i*((Complex::<T>::one() - self*self).sqrt() + i*self).ln()
     }
 
     /// Computes the principal value of the inverse cosine of `self`.
@@ -276,8 +276,8 @@ impl<T: Clone + Float> Complex<T> {
     #[inline]
     pub fn acos(&self) -> Complex<T> {
         // formula: arccos(z) = -i ln(i sqrt(1-z^2) + z)
-        let i = Complex::i();
-        -i*(i*(Complex::one() - self*self).sqrt() + self).ln()
+        let i = Complex::<T>::i();
+        -i*(i*(Complex::<T>::one() - self*self).sqrt() + self).ln()
     }
 
     /// Computes the principal value of the inverse tangent of `self`.
@@ -291,8 +291,8 @@ impl<T: Clone + Float> Complex<T> {
     #[inline]
     pub fn atan(&self) -> Complex<T> {
         // formula: arctan(z) = (ln(1+iz) - ln(1-iz))/(2i)
-        let i = Complex::i();
-        let one = Complex::one();
+        let i = Complex::<T>::i();
+        let one = Complex::<T>::one();
         let two = one + one;
         if *self == i {
             return Complex::new(T::zero(), T::infinity());
@@ -336,7 +336,7 @@ impl<T: Clone + Float> Complex<T> {
     #[inline]
     pub fn asinh(&self) -> Complex<T> {
         // formula: arcsinh(z) = ln(z + sqrt(1+z^2))
-        let one = Complex::one();
+        let one = Complex::<T>::one();
         (self + (one + self * self).sqrt()).ln()
     }
 
@@ -417,8 +417,8 @@ impl<'a, T: Clone + Num> From<&'a T> for Complex<T> {
 }
 
 macro_rules! forward_ref_ref_binop {
-    (impl $imp:ident, $method:ident) => {
-        impl<'a, 'b, T: Clone + Num> $imp<&'b Complex<T>> for &'a Complex<T> {
+    (impl $imp:ident, $method:ident, $($dep:ident),*) => {
+        impl<'a, 'b, T: Clone + Num $(+ $dep)*> $imp<&'b Complex<T>> for &'a Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -430,8 +430,8 @@ macro_rules! forward_ref_ref_binop {
 }
 
 macro_rules! forward_ref_val_binop {
-    (impl $imp:ident, $method:ident) => {
-        impl<'a, T: Clone + Num> $imp<Complex<T>> for &'a Complex<T> {
+    (impl $imp:ident, $method:ident, $($dep:ident),*) => {
+        impl<'a, T: Clone + Num $(+ $dep)*> $imp<Complex<T>> for &'a Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -443,8 +443,8 @@ macro_rules! forward_ref_val_binop {
 }
 
 macro_rules! forward_val_ref_binop {
-    (impl $imp:ident, $method:ident) => {
-        impl<'a, T: Clone + Num> $imp<&'a Complex<T>> for Complex<T> {
+    (impl $imp:ident, $method:ident, $($dep:ident),*) => {
+        impl<'a, T: Clone + Num $(+ $dep)*> $imp<&'a Complex<T>> for Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -456,15 +456,15 @@ macro_rules! forward_val_ref_binop {
 }
 
 macro_rules! forward_all_binop {
-    (impl $imp:ident, $method:ident) => {
-        forward_ref_ref_binop!(impl $imp, $method);
-        forward_ref_val_binop!(impl $imp, $method);
-        forward_val_ref_binop!(impl $imp, $method);
+    (impl $imp:ident, $method:ident, $($dep:ident),*) => {
+        forward_ref_ref_binop!(impl $imp, $method, $($dep),*);
+        forward_ref_val_binop!(impl $imp, $method, $($dep),*);
+        forward_val_ref_binop!(impl $imp, $method, $($dep),*);
     };
 }
 
 /* arithmetic */
-forward_all_binop!(impl Add, add);
+forward_all_binop!(impl Add, add, );
 
 // (a + i b) + (c + i d) == (a + c) + i (b + d)
 impl<T: Clone + Num> Add<Complex<T>> for Complex<T> {
@@ -476,7 +476,7 @@ impl<T: Clone + Num> Add<Complex<T>> for Complex<T> {
     }
 }
 
-forward_all_binop!(impl Sub, sub);
+forward_all_binop!(impl Sub, sub, );
 
 // (a + i b) - (c + i d) == (a - c) + i (b - d)
 impl<T: Clone + Num> Sub<Complex<T>> for Complex<T> {
@@ -488,7 +488,7 @@ impl<T: Clone + Num> Sub<Complex<T>> for Complex<T> {
     }
 }
 
-forward_all_binop!(impl Mul, mul);
+forward_all_binop!(impl Mul, mul, );
 
 // (a + i b) * (c + i d) == (a*c - b*d) + i (a*d + b*c)
 impl<T: Clone + Num> Mul<Complex<T>> for Complex<T> {
@@ -502,7 +502,7 @@ impl<T: Clone + Num> Mul<Complex<T>> for Complex<T> {
     }
 }
 
-forward_all_binop!(impl Div, div);
+forward_all_binop!(impl Div, div, );
 
 // (a + i b) / (c + i d) == [(a + i b) * (c - i d)] / (c*c + d*d)
 //   == [(a*c + b*d) / (c*c + d*d)] + i [(b*c - a*d) / (c*c + d*d)]
@@ -518,10 +518,50 @@ impl<T: Clone + Num> Div<Complex<T>> for Complex<T> {
     }
 }
 
+forward_all_binop!(impl Rem, rem, PartialOrd);
+
+// Attempts to identify the gaussian integer whose product with `modulus`
+// is closest to `self`.
+impl<T: Clone + Num + PartialOrd> Rem<Complex<T>> for Complex<T> where {
+    type Output = Complex<T>;
+
+    #[inline]
+    fn rem(self, modulus: Complex<T>) -> Self {
+        let Complex { re, im } = self.clone() / modulus.clone();
+        // This is the gaussian integer corresponding to the true ratio
+        // rounded towards zero.
+        let (re0, im0) = (re.clone() - re % T::one(), im.clone() - im % T::one());
+
+        let zero = T::zero();
+        let one  = T::one();
+        let neg  = zero.clone() - one.clone();
+        // Traverse the 3x3 square of gaussian integers surrounding our
+        // current approximation, and select the one whose product with
+        // `modulus` is closest to `self`.
+        let mut bestrem = self.clone() - modulus.clone() *
+                          Complex::new(re0.clone(), im0.clone());
+        let mut bestnorm = bestrem.norm_sqr();
+        for &(dr, di) in
+            vec![(&one, &zero), (&one, &one), (&zero, &one), (&neg, &one),
+                 (&neg, &zero), (&neg, &neg), (&zero, &neg), (&one, &neg)]
+                            .iter() {
+            let newrem = self.clone() - modulus.clone() *
+                         Complex::new(re0.clone() + dr.clone(),
+                                      im0.clone() + di.clone());
+            let newnorm = newrem.norm_sqr();
+            if newnorm < bestnorm {
+                bestrem = newrem;
+                bestnorm = newnorm;
+            }
+        }
+        bestrem
+    }
+}
+
 // Op Assign
 
 mod opassign {
-    use std::ops::{AddAssign, SubAssign, MulAssign, DivAssign};
+    use std::ops::{AddAssign, SubAssign, MulAssign, DivAssign, RemAssign};
 
     use traits::NumAssign;
 
@@ -553,6 +593,12 @@ mod opassign {
         }
     }
 
+    impl<T: Clone + NumAssign + PartialOrd> RemAssign for Complex<T> {
+        fn rem_assign(&mut self, other: Complex<T>) {
+            *self = self.clone() % other;
+        }
+    }
+
     impl<T: Clone + NumAssign> AddAssign<T> for Complex<T> {
         fn add_assign(&mut self, other: T) {
             self.re += other;
@@ -579,6 +625,12 @@ mod opassign {
         }
     }
 
+    impl<T: Clone + NumAssign + PartialOrd> RemAssign<T> for Complex<T> {
+        fn rem_assign(&mut self, other: T) {
+            *self = self.clone() % other;
+        }
+    }
+
     macro_rules! forward_op_assign {
         (impl $imp:ident, $method:ident) => {
             impl<'a, T: Clone + NumAssign> $imp<&'a Complex<T>> for Complex<T> {
@@ -600,6 +652,19 @@ mod opassign {
     forward_op_assign!(impl SubAssign, sub_assign);
     forward_op_assign!(impl MulAssign, mul_assign);
     forward_op_assign!(impl DivAssign, div_assign);
+
+    impl<'a, T: Clone + NumAssign + PartialOrd> RemAssign<&'a Complex<T>> for Complex<T> {
+        #[inline]
+        fn rem_assign(&mut self, other: &Complex<T>) {
+            self.rem_assign(other.clone())
+        }
+    }
+    impl<'a, T: Clone + NumAssign + PartialOrd> RemAssign<&'a T> for Complex<T> {
+        #[inline]
+        fn rem_assign(&mut self, other: &T) {
+            self.rem_assign(other.clone())
+        }
+    }
 }
 
 impl<T: Clone + Num + Neg<Output = T>> Neg for Complex<T> {
@@ -621,8 +686,8 @@ impl<'a, T: Clone + Num + Neg<Output = T>> Neg for &'a Complex<T> {
 }
 
 macro_rules! real_arithmetic {
-    (@forward $imp:ident::$method:ident for $($real:ident),*) => (
-        impl<'a, T: Clone + Num> $imp<&'a T> for Complex<T> {
+    (@forward $imp:ident::$method:ident for $($real:ident),*:  $($dep:ident),*) => (
+        impl<'a, T: Clone + Num $(+ $dep)*> $imp<&'a T> for Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -630,7 +695,7 @@ macro_rules! real_arithmetic {
                 self.$method(other.clone())
             }
         }
-        impl<'a, T: Clone + Num> $imp<T> for &'a Complex<T> {
+        impl<'a, T: Clone + Num $(+ $dep)*> $imp<T> for &'a Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -638,7 +703,7 @@ macro_rules! real_arithmetic {
                 self.clone().$method(other)
             }
         }
-        impl<'a, 'b, T: Clone + Num> $imp<&'a T> for &'b Complex<T> {
+        impl<'a, 'b, T: Clone + Num $(+ $dep)*> $imp<&'a T> for &'b Complex<T> {
             type Output = Complex<T>;
 
             #[inline]
@@ -674,10 +739,11 @@ macro_rules! real_arithmetic {
         )*
     );
     ($($real:ident),*) => (
-        real_arithmetic!(@forward Add::add for $($real),*);
-        real_arithmetic!(@forward Sub::sub for $($real),*);
-        real_arithmetic!(@forward Mul::mul for $($real),*);
-        real_arithmetic!(@forward Div::div for $($real),*);
+        real_arithmetic!(@forward Add::add for $($real),*: );
+        real_arithmetic!(@forward Sub::sub for $($real),*: );
+        real_arithmetic!(@forward Mul::mul for $($real),*: );
+        real_arithmetic!(@forward Div::div for $($real),*: );
+        real_arithmetic!(@forward Rem::rem for $($real),*: PartialOrd);
 
         $(
             impl Add<Complex<$real>> for $real {
@@ -718,6 +784,15 @@ macro_rules! real_arithmetic {
                                  $real::zero() - self * other.im / norm_sqr)
                 }
             }
+
+            impl Rem<Complex<$real>> for $real {
+                type Output = Complex<$real>;
+
+                #[inline]
+                fn rem(self, other: Complex<$real>) -> Complex<$real> {
+                    Complex::new(self, Self::zero()) % other
+                }
+            }
         )*
     );
 }
@@ -758,6 +833,15 @@ impl<T: Clone + Num> Div<T> for Complex<T> {
     }
 }
 
+impl<T: Clone + Num + PartialOrd> Rem<T> for Complex<T> {
+    type Output = Complex<T>;
+
+    #[inline]
+    fn rem(self, other: T) -> Complex<T> {
+        self % Complex::new(other, T::zero())
+    }
+}
+
 real_arithmetic!(usize, u8, u16, u32, u64, isize, i8, i16, i32, i64, f32, f64);
 
 /* constants */
@@ -970,6 +1054,96 @@ impl<T> FromStr for Complex<T> where
     }
 }
 
+impl<T: Num + Clone + PartialOrd> Num for Complex<T> where
+{
+    type FromStrRadixErr = ParseComplexError<T::FromStrRadixErr>;
+
+    /// Parses `a +/- bi`; `ai +/- b`; `a`; or `bi` where `a` and `b` are of type `T`
+    fn from_str_radix(s: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr>
+    {
+        let imag = match s.rfind('j') {
+            None => 'i',
+            _ => 'j'
+        };
+
+        let mut b = String::with_capacity(s.len());
+        let mut first = true;
+
+        let char_indices = s.char_indices();
+        let mut pc = ' ';
+        let mut split_index = s.len();
+
+        for (i, cc) in char_indices {
+            if cc == '+' && pc != 'e' && pc != 'E' && i > 0 {
+                // ignore '+' if part of an exponent
+                if first {
+                    split_index = i;
+                    first = false;
+                }
+                // don't carry '+' over into b
+                pc = ' ';
+                continue;
+            } else if cc == '-' && pc != 'e' && pc != 'E' && i > 0 {
+                // ignore '-' if part of an exponent or begins the string
+                if first {
+                    split_index = i;
+                    first = false;
+                }
+                // DO carry '-' over into b
+            }
+
+            if pc == '-' && cc == ' ' && !first {
+                // ignore whitespace between minus sign and next number
+                continue;
+            }
+
+            if !first {
+                b.push(cc);
+            }
+            pc = cc;
+        }
+
+        // split off real and imaginary parts, trim whitespace
+        let (a, _) = s.split_at(split_index);
+        let a = a.trim_right();
+        let mut b = b.trim_left();
+        // input was either pure real or pure imaginary
+        if b.is_empty() {
+            b = match a.ends_with(imag) {
+                false => "0i",
+                true => "0"
+            };
+        }
+
+        let re;
+        let im;
+        if a.ends_with(imag) {
+            im = a; re = b;
+        } else if b.ends_with(imag) {
+            re = a; im = b;
+        } else {
+            return Err(ParseComplexError::new());
+        }
+
+        // parse re
+        let re = try!(T::from_str_radix(re, radix).map_err(ParseComplexError::from_error));
+
+        // pop imaginary unit off
+        let mut im = &im[..im.len()-1];
+        // handle im == "i" or im == "-i"
+        if im.is_empty() || im == "+" {
+            im = "1";
+        } else if im == "-" {
+            im = "-1";
+        }
+
+        // parse im
+        let im = try!(T::from_str_radix(im, radix).map_err(ParseComplexError::from_error));
+
+        Ok(Complex::new(re, im))
+    }
+}
+
 #[cfg(feature = "serde")]
 impl<T> serde::Serialize for Complex<T>
     where T: serde::Serialize
@@ -1512,6 +1686,10 @@ mod test {
             assert_eq!($a / $b, $answer);
             assert_eq!({ let mut x = $a; x /= $b; x}, $answer);
         };
+        ($a:ident % $b:expr, $answer:expr) => {
+            assert_eq!($a % $b, $answer);
+            assert_eq!({ let mut x = $a; x %= $b; x}, $answer);
+        }
     }
 
     // Test both a + b and a + &b
@@ -1523,7 +1701,7 @@ mod test {
     }
 
     mod complex_arithmetic {
-        use super::{_0_0i, _1_0i, _1_1i, _0_1i, _neg1_1i, _05_05i, all_consts};
+        use super::{_0_0i, _1_0i, _1_1i, _0_1i, _neg1_1i, _05_05i, _4_2i, all_consts};
         use traits::Zero;
 
         #[test]
@@ -1575,6 +1753,16 @@ mod test {
             }
         }
 
+        #[test]
+        fn test_rem() {
+            test_op!(_neg1_1i % _0_1i, _0_0i);
+            test_op!(_4_2i % _0_1i, _0_0i);
+            test_op!(_05_05i % _0_1i, _05_05i);
+            test_op!(_05_05i % _1_1i, _05_05i);
+            assert_eq!((_4_2i + _05_05i) % _0_1i, _05_05i);
+            assert_eq!((_4_2i + _05_05i) % _1_1i, _05_05i);
+        }
+
         #[test]
         fn test_neg() {
             assert_eq!(-_1_0i + _0_1i, _neg1_1i);
@@ -1612,6 +1800,13 @@ mod test {
             assert_eq!(_4_2i / 0.5, Complex::new(8.0, 4.0));
             assert_eq!(0.5 / _4_2i, Complex::new(0.1, -0.05));
         }
+
+        #[test]
+        fn test_rem() {
+            assert_eq!(_4_2i % 2.0, Complex::new(0.0, 0.0));
+            assert_eq!(_4_2i % 3.0, Complex::new(1.0, -1.0));
+            assert_eq!(3.0 % _4_2i, Complex::new(-1.0, -2.0));
+        }
     }
 
     #[test]