diff --git a/src/systems/polynom.rs b/src/systems/polynom.rs
index 7e0e01c..47c43df 100644
--- a/src/systems/polynom.rs
+++ b/src/systems/polynom.rs
@@ -112,13 +112,13 @@ impl<T> Polynomial<T> {
     }
 }
 
-impl<T> Mul for Polynomial<T>
+impl<'a, T> Mul<&'a Polynomial<T>> for &Polynomial<T>
 where
     T: Default + Mul<Output = T> + AddAssign + Clone,
 {
     type Output = Polynomial<T>;
 
-    fn mul(self, rhs: Self) -> Self::Output {
+    fn mul(self, rhs: &'a Polynomial<T>) -> Self::Output {
         let mut result =
             vec![T::default(); self.coeffs.len() + rhs.coeffs.len() - 1];
         for (idx_l, val_l) in self.coeffs.iter().enumerate() {
@@ -131,13 +131,33 @@ where
     }
 }
 
-impl<T> Add for Polynomial<T>
+impl<T> Mul for Polynomial<T>
+where
+    T: Default + Mul<Output = T> + AddAssign + Clone,
+{
+    type Output = Polynomial<T>;
+
+    fn mul(self, rhs: Self) -> Self::Output {
+        &self * &rhs
+        // let mut result =
+        //     vec![T::default(); self.coeffs.len() + rhs.coeffs.len() - 1];
+        // for (idx_l, val_l) in self.coeffs.iter().enumerate() {
+        //     for (idx_r, val_r) in rhs.coeffs.iter().enumerate() {
+        //         result[idx_l + idx_r] += val_l.clone() * val_r.clone();
+        //     }
+        // }
+
+        // Polynomial { coeffs: result }
+    }
+}
+
+impl<'a, T> Add<&'a Polynomial<T>> for &Polynomial<T>
 where
     T: AddAssign + Default + Clone,
 {
     type Output = Polynomial<T>;
 
-    fn add(self, rhs: Self) -> Self::Output {
+    fn add(self, rhs: &'a Polynomial<T>) -> Self::Output {
         let mut result =
             vec![T::default(); self.coeffs.len().max(rhs.coeffs.len())];
         for (idx, val) in self.coeffs.iter().enumerate() {
@@ -151,27 +171,49 @@ where
     }
 }
 
-impl<T> Neg for Polynomial<T>
+impl<T> Add for Polynomial<T>
 where
-    T: Neg<Output = T> + Clone,
+    T: AddAssign + Default + Clone,
+{
+    type Output = Polynomial<T>;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        &self + &rhs
+    }
+}
+
+impl<T> Neg for &Polynomial<T>
+where
+    T: Neg<Output = T> + Clone + Copy,
 {
     type Output = Polynomial<T>;
     fn neg(self) -> Self::Output {
         Polynomial {
-            coeffs: self.coeffs.into_iter().map(|x| -x).collect(),
+            coeffs: self.coeffs.iter().map(|&x| -x).collect(),
         }
     }
 }
 
+impl<'a, T> Sub<&'a Polynomial<T>> for &Polynomial<T>
+where
+    T: Neg<Output = T> + Clone + AddAssign + Default + Copy,
+{
+    type Output = Polynomial<T>;
+
+    fn sub(self, rhs: &'a Polynomial<T>) -> Self::Output {
+        let neg_rhs = -rhs;
+        self + &neg_rhs
+    }
+}
+
 impl<T> Sub for Polynomial<T>
 where
-    T: Neg<Output = T> + Clone + AddAssign + Default,
+    T: Neg<Output = T> + Clone + AddAssign + Default + Copy,
 {
     type Output = Polynomial<T>;
 
     fn sub(self, rhs: Self) -> Self::Output {
-        let neg_rhs = -rhs;
-        self + neg_rhs
+        &self - &rhs
     }
 }
 
@@ -180,12 +222,21 @@ macro_rules! impl_compound_assign {
     $(
         impl<T> $assign_trait for $struct_type<T>
         where
-            T: $trait<Output = T> + $assign_trait + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add,
+            T: $trait<Output = T> + $assign_trait + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add + Copy,
         {
             fn $assign_method(&mut self, rhs: Self) {
                 *self = self.clone().$method(rhs)
             }
         }
+
+        impl<'a, T> $assign_trait<&'a $struct_type<T>> for $struct_type<T>
+        where
+            T: $trait<Output = T> + $assign_trait + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add + Copy,
+        {
+            fn $assign_method(&mut self, rhs: &'a $struct_type<T>) {
+                *self = <&Self as $trait<&Self>>::$method(self, rhs);
+            }
+        }
     )*
     };
 }
@@ -426,13 +477,13 @@ where
 
 impl<T> RationalFunction<T>
 where
-    T: One + Clone + Zero + Add + Mul<Output = T> + AddAssign + Default,
+    T: One + Clone + Zero + Add + Mul<Output = T> + AddAssign + Default + Copy,
 {
-    pub fn powi(self, exp: i32) -> Self {
+    pub fn powi(&self, exp: i32) -> Self {
         let base = self;
         let mut result = RationalFunction::new_from_scalar(T::one());
         for _ in 0..exp.abs() {
-            result = result * base.clone();
+            result = &result * base;
         }
         if exp < 0 {
             result = RationalFunction::new_from_scalar(T::one()) / result
@@ -459,32 +510,67 @@ where
     }
 }
 
+impl<'a, T> Mul<&'a RationalFunction<T>> for &RationalFunction<T>
+where
+    T: Add + Default + Mul<Output = T> + AddAssign + Clone + Copy,
+{
+    type Output = RationalFunction<T>;
+    fn mul(self, rhs: &'a RationalFunction<T>) -> Self::Output {
+        let new_num = &self.num * &rhs.num;
+        let new_den = &self.den * &rhs.den;
+
+        RationalFunction {
+            num: new_num,
+            den: new_den,
+        }
+    }
+}
+
 impl<T> Mul for RationalFunction<T>
 where
-    T: Add + Default + Mul<Output = T> + AddAssign + Clone,
+    T: Add + Default + Mul<Output = T> + AddAssign + Clone + Copy,
 {
     type Output = RationalFunction<T>;
     fn mul(self, rhs: Self) -> Self::Output {
-        let new_num = self.num * rhs.num;
-        let new_den = self.den * rhs.den;
+        &self * &rhs
+    }
+}
 
-        Self {
+impl<'a, T> Div<&'a RationalFunction<T>> for &RationalFunction<T>
+where
+    T: Clone + Default + Mul<Output = T> + AddAssign + Copy,
+{
+    type Output = RationalFunction<T>;
+    fn div(self, rhs: &'a RationalFunction<T>) -> Self::Output {
+        let new_num = &self.num * &rhs.den;
+        let new_den = &self.den * &rhs.num;
+
+        RationalFunction {
             num: new_num,
             den: new_den,
         }
     }
 }
-
 impl<T> Div for RationalFunction<T>
 where
-    T: Clone + Default + Mul<Output = T> + AddAssign,
+    T: Clone + Default + Mul<Output = T> + AddAssign + Copy,
 {
     type Output = RationalFunction<T>;
     fn div(self, rhs: Self) -> Self::Output {
-        let new_num = self.num * rhs.den;
-        let new_den = self.den * rhs.num;
+        &self / &rhs
+    }
+}
 
-        Self {
+impl<'a, T> Add<&'a RationalFunction<T>> for &RationalFunction<T>
+where
+    T: AddAssign + Default + Clone + Mul<Output = T>,
+{
+    type Output = RationalFunction<T>;
+    fn add(self, rhs: &'a RationalFunction<T>) -> Self::Output {
+        let new_den = &self.den * &rhs.den;
+        let new_num = &self.num * &rhs.den + &rhs.num * &self.den;
+
+        RationalFunction {
             num: new_num,
             den: new_den,
         }
@@ -497,38 +583,42 @@ where
 {
     type Output = RationalFunction<T>;
     fn add(self, rhs: Self) -> Self::Output {
-        let new_den = self.den.clone() * rhs.den.clone();
-        let new_num = self.num * rhs.den + rhs.num * self.den;
-
-        Self {
-            num: new_num,
-            den: new_den,
-        }
+        &self + &rhs
     }
 }
 
-impl<T> Neg for RationalFunction<T>
+impl<T> Neg for &RationalFunction<T>
 where
-    T: Neg<Output = T> + Clone,
+    T: Neg<Output = T> + Clone + Copy,
 {
     type Output = RationalFunction<T>;
     fn neg(self) -> Self::Output {
-        let new_num = -self.num;
+        let new_num = -&self.num;
 
-        Self {
+        // TODO: should not need to use clone here
+        RationalFunction {
             num: new_num,
-            den: self.den,
+            den: self.den.clone(),
         }
     }
 }
 
+impl<'a, T> Sub<&'a RationalFunction<T>> for &RationalFunction<T>
+where
+    T: Neg<Output = T> + Clone + AddAssign + Default + Mul<Output = T> + Copy,
+{
+    type Output = RationalFunction<T>;
+    fn sub(self, rhs: &'a RationalFunction<T>) -> Self::Output {
+        self + &(-rhs)
+    }
+}
 impl<T> Sub for RationalFunction<T>
 where
-    T: Neg<Output = T> + Clone + AddAssign + Default + Mul<Output = T>,
+    T: Neg<Output = T> + Clone + AddAssign + Default + Mul<Output = T> + Copy,
 {
     type Output = RationalFunction<T>;
     fn sub(self, rhs: Self) -> Self::Output {
-        self + (-rhs)
+        &self - &rhs
     }
 }
 
@@ -867,4 +957,17 @@ mod tests {
             RationalFunction::new_from_coeffs(&[1., 2., 4.], &[1., 4., 2., 8.]);
         assert_abs_diff_eq!(rf.clone(), rf);
     }
+
+    #[test]
+    fn reference_arithemtic() {
+        let mut rf1 =
+            RationalFunction::new_from_coeffs(&[1.0, 2.0], &[2.0, 1.0]);
+        let rf2 = rf1.clone();
+
+        rf1 += &rf2;
+        rf1 += &rf2;
+
+        let ans: RationalFunction<f64> = 3.0 * rf2;
+        assert_abs_diff_eq!(rf1.eval(&1.2), ans.eval(&1.2));
+    }
 }
diff --git a/src/systems/state_space.rs b/src/systems/state_space.rs
index b4407c0..81374a9 100644
--- a/src/systems/state_space.rs
+++ b/src/systems/state_space.rs
@@ -204,7 +204,7 @@ impl<U: Time + 'static> Ss<U> {
         let is_valid = self.is_valid();
         assert!(
             is_valid.is_ok(),
-            "Ss is invalid: Message: {}, Ss: {:?}",
+            "Ss is invalid: Message: {}, Ss: {}",
             is_valid.unwrap_err(),
             self
         );
@@ -278,7 +278,7 @@ impl<U: Time + 'static> Ss<U> {
     ///
     /// # Returns
     /// A new system representing the series connection of `self` and `sys2`.
-    pub fn series(self, sys2: Self) -> Self {
+    pub fn series(&self, sys2: &Self) -> Self {
         assert_eq!(
             self.noutputs(),
             sys2.ninputs(),
@@ -286,9 +286,37 @@ impl<U: Time + 'static> Ss<U> {
         );
         self.assert_valid();
         sys2.assert_valid();
-        let ret = sys2 * self;
-        ret.assert_valid();
-        ret
+
+        let nx2 = sys2.order();
+        let nu2 = sys2.ninputs();
+        let ny2 = sys2.noutputs();
+        let nx1 = self.order();
+        let nu1 = self.ninputs();
+        let ny1 = self.noutputs();
+        assert_eq!(nu2, ny1);
+
+        let mut a_new = DMatrix::zeros(nx1 + nx2, nx1 + nx2);
+        a_new.view_mut((0, 0), (nx1, nx1)).copy_from(self.a());
+        a_new
+            .view_mut((nx1, 0), (nx2, nx1))
+            .copy_from(&(sys2.b() * self.c()));
+        a_new.view_mut((nx1, nx1), (nx2, nx2)).copy_from(sys2.a());
+
+        let mut b_new = DMatrix::zeros(nx1 + nx2, nu1);
+        b_new.view_mut((0, 0), (nx1, nu1)).copy_from(self.b());
+        b_new
+            .view_mut((nx1, 0), (nx2, nu1))
+            .copy_from(&(sys2.b() * self.d()));
+
+        let mut c_new = DMatrix::zeros(ny2, nx1 + nx2);
+        c_new
+            .view_mut((0, 0), (ny2, nx1))
+            .copy_from(&(sys2.d() * self.c()));
+        c_new.view_mut((0, nx1), (ny2, nx2)).copy_from(sys2.c());
+
+        let d_new = sys2.d() * self.d();
+
+        Ss::<U>::new(a_new, b_new, c_new, d_new).unwrap()
     }
 
     /// Connects two systems in **parallel**.
@@ -305,7 +333,7 @@ impl<U: Time + 'static> Ss<U> {
     ///
     /// # Returns
     /// A new system representing the parallel connection of `self` and `sys2`.
-    pub fn parallel(self, sys2: Self) -> Self {
+    pub fn parallel(&self, sys2: &Self) -> Self {
         assert_eq!(
             self.shape(),
             sys2.shape(),
@@ -313,9 +341,28 @@ impl<U: Time + 'static> Ss<U> {
         );
         self.assert_valid();
         sys2.assert_valid();
-        let ret = self + sys2;
-        ret.assert_valid();
-        ret
+
+        let nu = self.ninputs();
+        let ny = self.noutputs();
+
+        let nx1 = self.order();
+        let nx2 = sys2.order();
+
+        let mut a_new = DMatrix::zeros(nx1 + nx2, nx1 + nx2);
+        a_new.view_mut((0, 0), (nx1, nx1)).copy_from(self.a());
+        a_new.view_mut((nx1, nx1), (nx2, nx2)).copy_from(sys2.a());
+
+        let mut b_new = DMatrix::zeros(nx1 + nx2, nu);
+        b_new.view_mut((0, 0), (nx1, nu)).copy_from(self.b());
+        b_new.view_mut((nx1, 0), (nx2, nu)).copy_from(sys2.b());
+
+        let mut c_new = DMatrix::zeros(ny, nx1 + nx2);
+        c_new.view_mut((0, 0), (ny, nx1)).copy_from(self.c());
+
+        c_new.view_mut((0, nx1), (ny, nx2)).copy_from(sys2.c());
+        let d_new = self.d() + sys2.d();
+
+        Ss::<U>::new(a_new, b_new, c_new, d_new).unwrap()
     }
 
     /// **Negative** Feedback connection of `self` with `sys2`.
@@ -414,7 +461,7 @@ impl<U: Time + 'static> Ss<U> {
     }
 }
 
-impl<U: Time + 'static> Add for Ss<U> {
+impl<'b, U: Time + 'static> Add<&'b Ss<U>> for &Ss<U> {
     type Output = Ss<U>;
 
     /// Adds two state-space systems in parallel.
@@ -423,32 +470,26 @@ impl<U: Time + 'static> Add for Ss<U> {
     ///
     /// # Panics
     /// Panics if the input-output dimensions of both systems do not match.
-    fn add(self, rhs: Self) -> Self::Output {
+    fn add(self, rhs: &'b Ss<U>) -> Self::Output {
         assert_eq!(self.shape(), rhs.shape());
         self.assert_valid();
         rhs.assert_valid();
 
-        let nu = self.ninputs();
-        let ny = self.noutputs();
-
-        let nx1 = self.order();
-        let nx2 = rhs.order();
-
-        let mut a_new = DMatrix::zeros(nx1 + nx2, nx1 + nx2);
-        a_new.view_mut((0, 0), (nx1, nx1)).copy_from(self.a());
-        a_new.view_mut((nx1, nx1), (nx2, nx2)).copy_from(rhs.a());
-
-        let mut b_new = DMatrix::zeros(nx1 + nx2, nu);
-        b_new.view_mut((0, 0), (nx1, nu)).copy_from(self.b());
-        b_new.view_mut((nx1, 0), (nx2, nu)).copy_from(rhs.b());
-
-        let mut c_new = DMatrix::zeros(ny, nx1 + nx2);
-        c_new.view_mut((0, 0), (ny, nx1)).copy_from(self.c());
+        self.parallel(rhs)
+    }
+}
 
-        c_new.view_mut((0, nx1), (ny, nx2)).copy_from(rhs.c());
-        let d_new = self.d() + rhs.d();
+impl<U: Time + 'static> Add for Ss<U> {
+    type Output = Ss<U>;
 
-        Ss::<U>::new(a_new, b_new, c_new, d_new).unwrap()
+    /// Adds two state-space systems in parallel.
+    ///
+    /// y = y1 + y2
+    ///
+    /// # Panics
+    /// Panics if the input-output dimensions of both systems do not match.
+    fn add(self, rhs: Self) -> Self::Output {
+        &self + &rhs
     }
 }
 
@@ -459,13 +500,26 @@ impl<U: Time + 'static> Neg for Ss<U> {
     ///
     /// I.e. y_new = -y    
     fn neg(mut self) -> Self::Output {
-        self.c = -self.c();
-        self.d = -self.d();
+        self.c *= -1.0;
+        self.d *= -1.0;
         self.assert_valid();
         self
     }
 }
 
+impl<'b, U: Time + 'static> Sub<&'b Ss<U>> for &Ss<U> {
+    type Output = Ss<U>;
+
+    /// Subtracts two state-space systems by adding one and negating the other.
+    fn sub(self, rhs: &'b Ss<U>) -> Self::Output {
+        self.assert_valid();
+        rhs.assert_valid();
+
+        // TODO: should not need to use clone
+        self + &(-rhs.clone())
+    }
+}
+
 impl<U: Time + 'static> Sub for Ss<U> {
     type Output = Ss<U>;
 
@@ -477,7 +531,7 @@ impl<U: Time + 'static> Sub for Ss<U> {
     }
 }
 
-impl<U: Time + 'static> Mul for Ss<U> {
+impl<'b, U: Time + 'static> Mul<&'b Ss<U>> for &Ss<U> {
     type Output = Ss<U>;
 
     /// Multiplies two state-space systems in series (cascading connection).
@@ -488,40 +542,49 @@ impl<U: Time + 'static> Mul for Ss<U> {
     /// # Panics
     /// Panics if the output size of `rhs` does not match the input size of
     /// `self`.
-    fn mul(self, rhs: Self) -> Self::Output {
+    fn mul(self, rhs: &'b Ss<U>) -> Self::Output {
         self.assert_valid();
         rhs.assert_valid();
 
-        let nx2 = self.order();
-        let nu2 = self.ninputs();
-        let ny2 = self.noutputs();
-        let nx1 = rhs.order();
-        let nu1 = rhs.ninputs();
-        let ny1 = rhs.noutputs();
-        assert_eq!(nu2, ny1);
+        rhs.series(self)
+    }
+}
 
-        let mut a_new = DMatrix::zeros(nx1 + nx2, nx1 + nx2);
-        a_new.view_mut((0, 0), (nx1, nx1)).copy_from(rhs.a());
-        a_new
-            .view_mut((nx1, 0), (nx2, nx1))
-            .copy_from(&(self.b() * rhs.c()));
-        a_new.view_mut((nx1, nx1), (nx2, nx2)).copy_from(self.a());
+impl<U: Time + 'static> Mul for Ss<U> {
+    type Output = Ss<U>;
 
-        let mut b_new = DMatrix::zeros(nx1 + nx2, nu1);
-        b_new.view_mut((0, 0), (nx1, nu1)).copy_from(rhs.b());
-        b_new
-            .view_mut((nx1, 0), (nx2, nu1))
-            .copy_from(&(self.b() * rhs.d()));
+    /// Multiplies two state-space systems in series (cascading connection).
+    ///
+    /// The resulting system represents `y <- self <- rhs <- u`, where `self`
+    /// follows `rhs`.
+    ///
+    /// # Panics
+    /// Panics if the output size of `rhs` does not match the input size of
+    /// `self`.
+    fn mul(self, rhs: Self) -> Self::Output {
+        &self * &rhs
+    }
+}
 
-        let mut c_new = DMatrix::zeros(ny2, nx1 + nx2);
-        c_new
-            .view_mut((0, 0), (ny2, nx1))
-            .copy_from(&(self.d() * rhs.c()));
-        c_new.view_mut((0, nx1), (ny2, nx2)).copy_from(self.c());
+impl<'b, U: Time + 'static> Div<&'b Ss<U>> for &Ss<U> {
+    type Output = Ss<U>;
 
-        let d_new = self.d() * rhs.d();
+    /// Divides one state-space system by another (right multiplication by
+    /// inverse).
+    ///
+    /// # Warning
+    /// The need for division often arrised due to feedback connections. Such as
+    /// T = L/(1 + L). For such cases use the `feedback` function, which is more
+    /// numerically stable.
+    ///
+    /// # Panics
+    /// Panics if the inverse of `rhs` cannot be computed.
+    #[allow(clippy::suspicious_arithmetic_impl)]
+    fn div(self, rhs: &'b Ss<U>) -> Self::Output {
+        self.assert_valid();
+        rhs.assert_valid();
 
-        Ss::<U>::new(a_new, b_new, c_new, d_new).unwrap()
+        self * &rhs.inv().unwrap()
     }
 }
 
@@ -540,9 +603,7 @@ impl<U: Time + 'static> Div for Ss<U> {
     /// Panics if the inverse of `rhs` cannot be computed.
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn div(self, rhs: Self) -> Self::Output {
-        self.assert_valid();
-        rhs.assert_valid();
-        self * rhs.inv().unwrap()
+        &self / &rhs
     }
 }
 
@@ -557,6 +618,14 @@ macro_rules! impl_compound_assign {
                 *self = self.clone().$method(rhs)
             }
         }
+
+        impl<'a, U: Time + 'static> $assign_trait<&'a $struct_type<U>> for $struct_type<U>
+        {
+            fn $assign_method(&mut self, rhs: &'a $struct_type<U>) {
+                *self = <&Self as $trait<&Self>>::$method(self, rhs);
+            }
+        }
+
     )*
     };
 }
@@ -621,10 +690,12 @@ mod tests {
     use core::f64;
 
     use approx::assert_abs_diff_eq;
+    use num_complex::c64;
     use rand::Rng;
 
     use super::*;
     use crate::{
+        FrequencyResponse,
         analysis::frequency_response::lin_space,
         systems::Tf,
         transformations::SsRealization::{ControllableCF, ObservableCF},
@@ -854,11 +925,11 @@ mod tests {
         assert_abs_diff_eq!(tf_fb, 0.5 * Tf::s() / (Tf::s() + 0.5));
 
         let ss_add = ss1.clone() + ss2.clone();
-        let ss_parallel = ss1.clone().parallel(ss2.clone());
+        let ss_parallel = ss1.parallel(&ss2);
         assert_eq!(ss_add, ss_parallel);
 
         let ss_mul = ss2.clone() * ss1.clone();
-        let ss_series = ss1.series(ss2);
+        let ss_series = ss1.series(&ss2);
         assert_eq!(ss_mul, ss_series);
     }
 
@@ -875,7 +946,7 @@ mod tests {
 
             let ss_fb = ss1.feedback(ss2);
             let ss_fb = ss_fb.to_tf().unwrap();
-            let tf_fb = tf1.feedback(tf2);
+            let tf_fb = tf1.feedback(&tf2);
 
             assert_abs_diff_eq!(tf_fb, ss_fb, epsilon = 1e-1);
         }
@@ -894,4 +965,22 @@ mod tests {
             .unwrap();
         println!("3: {sys}after");
     }
+
+    #[test]
+    fn ss_reference_arithmetic() {
+        let mut ss = (1.0 / Tf::s().powi(2)).to_ss().unwrap();
+        let ss_org = ss.clone();
+
+        ss += &ss_org;
+        ss += &ss_org;
+
+        let ans = (ss_org.clone() + ss_org.clone()) + ss_org.clone();
+
+        let freq = c64(0.0, 1.2);
+        let resp1 = ss.freq_response(&freq);
+        let resp2 = ans.freq_response(&freq);
+
+        assert_abs_diff_eq!(resp1.re, resp2.re, epsilon = 1e-9);
+        assert_abs_diff_eq!(resp1.im, resp2.im, epsilon = 1e-9);
+    }
 }
diff --git a/src/systems/transfer_function.rs b/src/systems/transfer_function.rs
index 3719338..c0cc6cc 100644
--- a/src/systems/transfer_function.rs
+++ b/src/systems/transfer_function.rs
@@ -223,9 +223,10 @@ where
 macro_rules! impl_operator_tf {
     ([$(($trait:ident, $method:ident)), *]) => {
     $(
+        // Value implementation
         impl<T, U> $trait for Tf<T, U>
         where
-            T: $trait<Output = T> + Clone + Zero + One+ Default + Add + AddAssign + Mul<Output = T> + Neg<Output = T>,
+            T: $trait<Output = T> + Clone + Zero + One+ Default + Add + AddAssign + Mul<Output = T> + Neg<Output = T> + Copy,
             U: Time,
         {
             type Output = Tf<T, U>;
@@ -234,19 +235,44 @@ macro_rules! impl_operator_tf {
                 Tf::new_from_rf(new_rf)
             }
         }
+
+        // Reference Implementation
+        impl<'a, T, U> $trait<&'a Tf<T, U>> for &Tf<T, U>
+        where
+            T: $trait<Output = T> + Clone + Zero + One+ Default + Add + AddAssign + Mul<Output = T> + Neg<Output = T> + Copy,
+            U: Time,
+        {
+            type Output = Tf<T, U>;
+            fn $method(self, rhs: &'a Tf<T, U>) -> Self::Output {
+                let new_rf = <&RationalFunction<T> as $trait<&RationalFunction<T>>>::$method(&self.rf, &rhs.rf);
+                Tf::new_from_rf(new_rf)
+            }
+        }
     )*
     };
 }
 impl_operator_tf!([(Add, add), (Sub, sub), (Mul, mul), (Div, div)]);
 
+impl<T, U> Neg for &Tf<T, U>
+where
+    T: Neg<Output = T> + Clone + Zero + One + Copy,
+    U: Time,
+{
+    type Output = Tf<T, U>;
+    fn neg(self) -> Self::Output {
+        let new_rf = -&self.rf;
+        Tf::new_from_rf(new_rf)
+    }
+}
+
 impl<T, U> Neg for Tf<T, U>
 where
-    T: Neg<Output = T> + Clone + Zero + One,
+    T: Neg<Output = T> + Clone + Zero + One + Copy,
     U: Time,
 {
     type Output = Tf<T, U>;
     fn neg(self) -> Self::Output {
-        let new_rf = -self.rf;
+        let new_rf = -&self.rf;
         Tf::new_from_rf(new_rf)
     }
 }
@@ -256,12 +282,21 @@ macro_rules! impl_compound_assign {
     $(
         impl<T, U: Time> $assign_trait for $struct_type<T, U>
         where
-            T: $trait<Output = T> + $assign_trait + One + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add + Zero,
+            T: $trait<Output = T> + $assign_trait + One + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add + Zero + Copy,
         {
             fn $assign_method(&mut self, rhs: Self) {
                 *self = self.clone().$method(rhs)
             }
         }
+
+        impl<'a, T, U: Time> $assign_trait<&'a $struct_type<T, U>> for $struct_type<T, U>
+        where
+            T: $trait<Output = T> + $assign_trait + One + Clone + Default + AddAssign + Neg<Output = T> + Mul<Output = T> + Add + Zero + Copy,
+        {
+            fn $assign_method(&mut self, rhs: &$struct_type<T, U>) {
+                *self = <&Self as $trait<&Self>>::$method(self, rhs);
+            }
+        }
     )*
     };
 }
@@ -368,7 +403,7 @@ where
 
 impl<T, U: Time> Tf<T, U>
 where
-    T: One + Clone + Zero + Add + Mul<Output = T> + AddAssign + Default,
+    T: One + Clone + Zero + Add + Mul<Output = T> + AddAssign + Default + Copy,
 {
     /// Raises the transfer function to a specified integer power.
     ///
@@ -377,7 +412,7 @@ where
     ///
     /// # Returns
     /// A new transfer function raised to the specified power.
-    pub fn powi(self, exp: i32) -> Self {
+    pub fn powi(&self, exp: i32) -> Self {
         let new_rf = self.rf.powi(exp);
         Tf::new_from_rf(new_rf)
     }
@@ -386,6 +421,7 @@ where
 impl<T, U: Time> Tf<T, U>
 where
     T: Clone
+        + Copy
         + Add<T, Output = T>
         + Zero
         + One
@@ -407,7 +443,7 @@ where
     ///
     /// # Returns
     /// A new system representing the series connection of `self` and `sys2`.
-    pub fn series(self, sys2: Self) -> Self {
+    pub fn series(&self, sys2: &Self) -> Self {
         sys2 * self
     }
 
@@ -425,7 +461,7 @@ where
     ///
     /// # Returns
     /// A new system representing the parallel connection of `self` and `sys2`.
-    pub fn parallel(self, sys2: Self) -> Self {
+    pub fn parallel(&self, sys2: &Self) -> Self {
         self + sys2
     }
 
@@ -448,7 +484,7 @@ where
     /// # Returns
     /// A new system representing the negative feedback connection of `self`
     /// with `sys2`.
-    pub fn feedback(self, sys2: Self) -> Self {
+    pub fn feedback(&self, sys2: &Self) -> Self {
         let n1 = self.numerator();
         let d1 = self.denominator();
         let n2 = sys2.numerator();
@@ -532,7 +568,7 @@ mod tests {
     #[test]
     fn tf_math() {
         let tf = Tf::s() + 1. / (Tf::s() + 1.0);
-        let tf_neg = -tf.clone();
+        let tf_neg = -&tf;
         let tf_neg_neg = -tf_neg;
         assert_eq!(tf, tf_neg_neg);
     }
@@ -571,11 +607,11 @@ mod tests {
         let tf2 = 1.0 / Tf::s();
 
         let tf_add = tf1.clone() + tf2.clone();
-        let tf_parallel = tf1.clone().parallel(tf2.clone());
+        let tf_parallel = tf1.clone().parallel(&tf2);
         assert_eq!(tf_add, tf_parallel);
 
         let tf_mul = tf2.clone() * tf1.clone();
-        let tf_series = tf1.clone().series(tf2.clone());
+        let tf_series = tf1.clone().series(&tf2);
         assert_eq!(tf_series, tf_mul);
     }
 
@@ -589,4 +625,17 @@ mod tests {
             2.1 / Tf::s() * (Tf::s() - 1.2).powi(5) / (Tf::s() + 1.0).powi(3);
         println!("3: {tf}");
     }
+
+    #[test]
+    fn reference_arithemtic() {
+        let mut tf = 1.0 / Tf::s();
+        let tf_org = tf.clone();
+
+        tf += &tf_org;
+        tf += &tf_org;
+
+        let ans = 3.0 * tf_org;
+
+        assert_abs_diff_eq!(tf.eval(&1.2), ans.eval(&1.2), epsilon = 1e-9);
+    }
 }