diff --git a/src/lib.rs b/src/lib.rs
index 8226c85..dd73a44 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,6 +1,5 @@
 extern crate core;
 
-pub mod linalg;
-
-
+use crate::linalg::rarray::{Dim, Rarray};
 
+pub mod linalg;
\ No newline at end of file
diff --git a/src/linalg.rs b/src/linalg.rs
index e1eaef7..524eceb 100644
--- a/src/linalg.rs
+++ b/src/linalg.rs
@@ -8,3 +8,4 @@ pub mod rarray1d_create;
 mod rarray2d_impl;
 pub mod rarray2d_ops;
 mod dimension;
+mod numeric_trait;
diff --git a/src/linalg/dimension.rs b/src/linalg/dimension.rs
index 96dba59..63a01b4 100644
--- a/src/linalg/dimension.rs
+++ b/src/linalg/dimension.rs
@@ -29,7 +29,7 @@ impl PartialEq<Self> for D1 {
     }
 
     fn ne(&self, other: &Self) -> bool {
-        self.width == other.width || self.height == other.height
+        self.width != other.width || self.height != other.height
     }
 }
 
@@ -60,7 +60,7 @@ impl PartialEq<Self> for D2 {
     }
 
     fn ne(&self, other: &Self) -> bool {
-        self.width == other.width || self.height == other.height
+        self.width != other.width || self.height != other.height
     }
 }
 
@@ -93,7 +93,7 @@ impl PartialEq<Self> for D3 {
     }
 
     fn ne(&self, other: &Self) -> bool {
-        self.width == other.width || self.height == other.height || self.depth == other.depth
+        self.width != other.width || self.height != other.height || self.depth != other.depth
     }
 }
 
diff --git a/src/linalg/numeric_trait.rs b/src/linalg/numeric_trait.rs
new file mode 100644
index 0000000..d4fb05e
--- /dev/null
+++ b/src/linalg/numeric_trait.rs
@@ -0,0 +1,25 @@
+use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign};
+
+/// Specifies that a generic type is a numerical type
+pub trait Numeric: 
+    Add<Output = Self> 
+    + AddAssign 
+    + Sub<Output = Self>
+    + SubAssign 
+    + Div<Output = Self>
+    + DivAssign 
+    + Mul<Output = Self>
+    + MulAssign 
+    + Copy 
+    + Default 
+    + ToString
+    + From<i32>
+    {}
+
+impl<T> Numeric for T where
+    T:
+    Add<Output = T> + AddAssign +
+    Sub<Output = T> + SubAssign +
+    Div<Output = T> + DivAssign +
+    Mul<Output = T> + MulAssign + 
+    ToString + Copy + Default + From<i32> {}
\ No newline at end of file
diff --git a/src/linalg/rarray.rs b/src/linalg/rarray.rs
index bf2b09d..d47ebc3 100644
--- a/src/linalg/rarray.rs
+++ b/src/linalg/rarray.rs
@@ -1,26 +1,25 @@
 use std::fmt::Debug;
-use std::ops::{Add, Mul, MulAssign, Sub};
+use std::ops::Neg;
 use std::usize;
-use rand::seq::IndexedRandom;
-
+use crate::linalg::numeric_trait::Numeric;
 pub(crate) use super::dimension::{Dim, D1, D2, D3};
 
-// Base array struct 
+/// Base array struct
+/// Consists of data field which is a 1 dimensional `Vec<T>` and of a shape field which is of type `Dim`
 #[derive(Debug)]
-pub struct Rarray<T, D> {
+pub struct Rarray<T, D: Dim> {
     pub(crate) data: Vec<T>,
     pub(crate) shape: D
 }
 
 // Specific implementations 
-pub type Rarray1D = Rarray<f64, D1>;
-pub type Rarray2D = Rarray<f64, D2>;
-pub type Rarray3D = Rarray<f64, D3>;
+pub type Rarray1D<T> = Rarray<T, D1>;
+pub type Rarray2D<T> = Rarray<T, D2>;
+pub type Rarray3D<T> = Rarray<T, D3>;
 
 pub trait RarrayCreate<T, V, S> {
     fn new(data: &V) -> Self;
     fn zeros(shape: T) -> Self;
-    fn random(shape: T) -> Self;
     fn fill(value: S, shape: T) -> Self;
 }
 
@@ -40,9 +39,11 @@ pub trait RarraySub<T, V, S> {
     fn sub(one: &T, other: &V) -> S;
 }
 
-impl RarrayMul<Rarray1D, Rarray2D, Rarray1D> for Rarray2D {
+impl<T> RarrayMul<Rarray1D<T>, Rarray2D<T>, Rarray1D<T>> for Rarray2D<T> where
+    T: Numeric
+{
     /// Performs (1 x n) x (n x m) matrix multiplication
-    fn mul(one: &Rarray1D, other: &Rarray2D) -> Rarray1D {
+    fn mul(one: &Rarray1D<T>, other: &Rarray2D<T>) -> Rarray1D<T> {
         let mut major: usize = 1;
         if one.shape.width > one.shape.height {
             assert_eq!(one.shape.width, other.shape.height, "Rarray shape mismatch");
@@ -54,24 +55,24 @@ impl RarrayMul<Rarray1D, Rarray2D, Rarray1D> for Rarray2D {
 
         let mut result = Rarray1D {
             shape: D1 { width: one.shape.width, height: 1 },
-            data: vec![0.; major]
+            data: vec![T::default(); major]
         };
 
         for i in 0..one.shape.width {
-            let mut sum: f64 = 0.;
             for j in 0..major {
-                sum += one[j] * other[[i, j]];
+                result[i] += one[j] * other[[i, j]];
             }
-            result[i] = sum;
         }
 
         result
     }
 }
 
-impl RarrayMul<Rarray2D, Rarray1D, Rarray1D> for Rarray2D {
+impl<T> RarrayMul<Rarray2D<T>, Rarray1D<T>, Rarray1D<T>> for Rarray2D<T> where
+    T: Numeric
+{
     /// Performs (n x m) x (m x 1) matrix multiplication
-    fn mul(one: &Rarray2D, other: &Rarray1D) -> Rarray1D {
+    fn mul(one: &Rarray2D<T>, other: &Rarray1D<T>) -> Rarray1D<T> {
         let mut major: usize = 1;
         if one.shape.width > one.shape.height {
             assert_eq!(one.shape.width, other.shape.height, "Rarray shape mismatch");
@@ -83,38 +84,36 @@ impl RarrayMul<Rarray2D, Rarray1D, Rarray1D> for Rarray2D {
 
         let mut result = Rarray1D {
             shape: D1 { width: one.shape.height, height: 1 },
-            data: vec![0.; major]
+            data: vec![T::default(); major]
         };
 
         for i in 0..one.shape.height {
-            let mut sum: f64 = 0.;
             for j in 0..major {
-                sum += one[[i, j]] * other[j];
+                result[i] += one[[i, j]] * other[j];
             }
-            result[i] = sum;
         }
 
         result
     }
 }
 
-impl RarrayMul<Rarray2D, Rarray2D, Rarray2D> for Rarray2D {
+impl<T> RarrayMul<Rarray2D<T>, Rarray2D<T>, Rarray2D<T>> for Rarray2D<T> where
+    T: Numeric
+{
     /// Performs (n x m) x (m x l) matrix multiplication
-    fn mul(one: &Rarray2D, other: &Rarray2D) -> Rarray2D {
+    fn mul(one: &Rarray2D<T>, other: &Rarray2D<T>) -> Rarray2D<T> {
         assert_eq!(one.shape.height, other.shape.width, "Rarray shape mismatch");
 
         let mut result = Rarray2D {
             shape: D2 { height: one.shape.height, width: other.shape.width },
-            data: vec![0.; one.shape.height * other.shape.width]
+            data: vec![T::default(); one.shape.height * other.shape.width]
         };
 
         for i in 0..one.shape.height {
             for j in 0..other.shape.width {
-                let mut sum: f64 = 0.;
                 for k in 0..one.shape.width {
-                    sum += one.data[i * one.shape.width + k] * other.data[j * other.shape.height + k];
+                    result.data[i * result.shape.width + j * result.shape.height] += one.data[i * one.shape.width + k] * other.data[j * other.shape.height + k];
                 }
-                result.data[i * result.shape.width + j * result.shape.height] = sum;
             }
         }
 
@@ -123,7 +122,7 @@ impl RarrayMul<Rarray2D, Rarray2D, Rarray2D> for Rarray2D {
 }
 
 impl<T, D> RarrayScalMul<T, Rarray<T, D>> for Rarray<T, D> where
-    T: Copy + MulAssign,
+    T: Numeric,
     D : Copy + Dim + Debug,
 {
     fn scal_mul(scal: T, rarray: &Rarray<T, D>) -> Rarray<T, D> {
@@ -141,7 +140,7 @@ impl<T, D> RarrayScalMul<T, Rarray<T, D>> for Rarray<T, D> where
 }
 
 impl<T, D> RarrayAdd<Rarray<T, D>, Rarray<T, D>, Rarray<T, D>> for Rarray<T, D> where
-    T : Copy + Add<Output = T> + Default,
+    T : Numeric,
     D : Copy + Dim + Debug + Eq
 {
     fn add(one: &Rarray<T, D>, other: &Rarray<T, D>) -> Rarray<T, D> {
@@ -161,7 +160,7 @@ impl<T, D> RarrayAdd<Rarray<T, D>, Rarray<T, D>, Rarray<T, D>> for Rarray<T, D>
 }
 
 impl<T, D> RarraySub<Rarray<T, D>, Rarray<T, D>, Rarray<T, D>> for Rarray<T, D> where
-    T : Copy + Sub<Output = T> + Default,
+    T : Numeric,
     D : Copy + Dim + Debug + Eq
 {
     fn sub(one: &Rarray<T, D>, other: &Rarray<T, D>) -> Rarray<T, D> {
@@ -178,4 +177,4 @@ impl<T, D> RarraySub<Rarray<T, D>, Rarray<T, D>, Rarray<T, D>> for Rarray<T, D>
 
         result
     }
-}
+}
\ No newline at end of file
diff --git a/src/linalg/rarray1d_create.rs b/src/linalg/rarray1d_create.rs
index 0b5a41c..2757435 100644
--- a/src/linalg/rarray1d_create.rs
+++ b/src/linalg/rarray1d_create.rs
@@ -1,9 +1,12 @@
 use std::cmp::max;
 use crate::linalg::dimension::{D1, D2};
+use crate::linalg::numeric_trait::Numeric;
 use crate::linalg::rarray::{Rarray1D, Rarray2D, RarrayCreate};
 
 // Creation
-impl RarrayCreate<usize, Vec<f64>, f64> for Rarray1D {
+impl<T> RarrayCreate<usize, Vec<T>, T> for Rarray1D<T> where
+    T: Numeric
+{
     /// Create rarray1D vector using Vec
     ///
     /// # Examples
@@ -12,16 +15,17 @@ impl RarrayCreate<usize, Vec<f64>, f64> for Rarray1D {
     /// use rumpy::linalg::rarray::Rarray1D;
     /// use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let v = Rarray1D::new(&vec![1., 1., 1.]);
+    /// let v = Rarray1D::<f64>::new(&vec![1., 1., 1.]);
     /// println!("{}", v);
     /// ```
-    fn new(data: &Vec<f64>) -> Self {
+    fn new(data: &Vec<T>) -> Self {
         Rarray1D {
-            shape : D1 { width: data.len() as usize, height: 1 },
+            shape : D1 { width: data.len(), height: 1 },
             data: data.clone()
         }
     }
 
+    //TODO: Constraint to numeric values
     /// Create rarra1d vector of length `shape` filled with zeros
     ///
     /// # Examples
@@ -30,37 +34,13 @@ impl RarrayCreate<usize, Vec<f64>, f64> for Rarray1D {
     ///use rumpy::linalg::rarray::Rarray1D;
     ///use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let v = Rarray1D::zeros(3);
+    /// let v = Rarray1D::<f64>::zeros(3);
     /// println!("{}", v);
     /// ```
     fn zeros(shape: usize) -> Self {
         Rarray1D {
             shape: D1 { width: shape, height: 1},
-            data: vec![0.; shape]
-        }
-    }
-
-    /// Create rarray1D vector of length `shape` filled with random numbers
-    ///
-    /// # Examples
-    ///
-    /// ```
-    /// use rumpy::linalg::rarray::Rarray1D;
-    /// use crate::rumpy::linalg::rarray::RarrayCreate;
-    ///
-    /// let v = Rarray1D::random(3);
-    /// println!("{}", v);
-    /// ```
-    fn random(shape: usize) -> Self {
-        let mut data = Vec::<f64>::with_capacity(shape);
-        for _ in 0..shape {
-            data.push(rand::random::<f64>());
-        }
-
-
-        Rarray1D {
-            shape: D1 { width: shape, height: 1 },
-            data
+            data: vec![T::default(); shape]
         }
     }
 
@@ -75,7 +55,7 @@ impl RarrayCreate<usize, Vec<f64>, f64> for Rarray1D {
     /// let v = Rarray1D::fill(4., 3);
     /// println!("{}", v);
     /// ```
-    fn fill(value: f64, shape: usize) -> Self {
+    fn fill(value: T, shape: usize) -> Self {
         Rarray1D {
             shape: D1 { width: shape, height: 1},
             data: vec![value; shape]
@@ -83,7 +63,9 @@ impl RarrayCreate<usize, Vec<f64>, f64> for Rarray1D {
     }
 }
 
-impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
+impl<T> RarrayCreate<(usize, usize), Vec<Vec<T>>, T> for Rarray1D<T> where
+    T: Numeric
+{
     /// Create vector using Vec<Vec>
     ///
     /// # Examples
@@ -92,14 +74,14 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
     /// use rumpy::linalg::rarray::Rarray1D;
     /// use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let v = Rarray1D::new(&vec![vec![0.], vec![0.], vec![0.]]);
+    /// let v = Rarray1D::<f64>::new(&vec![vec![0.], vec![0.], vec![0.]]);
     /// println!("{}", v);
     /// ```
     ///
     /// # Panics
     ///
     /// If rows aren't of same length
-    fn new(data: &Vec<Vec<f64>>) -> Self {
+    fn new(data: &Vec<Vec<T>>) -> Self {
         let row = data.len();
         let col = data[0].len();
         for i in 0..row {
@@ -120,7 +102,7 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
     /// use rumpy::linalg::rarray::Rarray1D;
     /// use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    ///let v = Rarray1D::zeros((3, 1));
+    ///let v = Rarray1D::<f64>::zeros((3, 1));
     /// println!("{}", v);
     /// ```
     ///
@@ -131,35 +113,7 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
         assert!(shape.0 == 1 || shape.1 == 1, "Cannot create 2D array using 1D array type");
         Rarray1D {
             shape: D1 { width: shape.1, height: shape.0},
-            data: vec![0.; shape.0 * shape.1]
-        }
-    }
-
-    /// Create vector of shape` filled with random numbers
-    ///
-    /// # Examples
-    ///
-    /// ```
-    /// use rumpy::linalg::rarray::Rarray1D;
-    /// use crate::rumpy::linalg::rarray::RarrayCreate;
-    ///
-    /// let v = Rarray1D::random((3, 1));
-    /// println!("{}", v);
-    /// ```
-    /// # Panics
-    ///
-    /// If (x, y) with x > 1 and y > 1
-    fn random(shape: (usize, usize)) -> Self {
-        assert!(shape.0 == 1 || shape.1 == 1, "Cannot create 2D array using 1D array type");
-        let mut data = Vec::<f64>::with_capacity(shape.0 * shape.1);
-        for _ in 0..(shape.0 * shape.1) {
-            data.push(rand::random::<f64>());
-        }
-
-
-        Rarray1D {
-            shape: D1 { width: shape.1, height: shape.0 },
-            data
+            data: vec![T::default(); shape.0 * shape.1]
         }
     }
 
@@ -171,13 +125,13 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
     /// use rumpy::linalg::rarray::Rarray1D;
     /// use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let v = Rarray1D::fill(4., (3, 1));
+    /// let v = Rarray1D::<f64>::fill(4., (3, 1));
     /// println!("{}", v);
     /// ```
     /// # Panics
     ///
     /// If (x, y) with x > 1 and y > 1
-    fn fill(value: f64, shape: (usize, usize)) -> Self {
+    fn fill(value: T, shape: (usize, usize)) -> Self {
         assert!(shape.0 == 1 || shape.1 == 1, "Cannot create 2D array using 1D array type");
         Rarray1D {
             shape: D1 { width: shape.1, height: shape.0},
@@ -186,7 +140,9 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray1D {
     }
 }
 
-impl Rarray1D {
+impl<T> Rarray1D<T> where
+    T: Numeric
+{
     /// Create matrix where given vector are the diagonal elements
     ///
     /// # Examples
@@ -200,12 +156,12 @@ impl Rarray1D {
     /// println!("{}", res);
     /// // >> Rarray2D([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
     /// ```
-    pub fn diag(&self) -> Rarray2D {
+    pub fn diag(&self) -> Rarray2D<T> {
         let major = max(self.shape.height, self.shape.width);
 
         let mut result = Rarray2D {
             shape: D2 { height: major, width: major},
-            data: vec![0.; major * major]
+            data: vec![T::default(); major * major]
         };
 
         for i in 0..major {
@@ -215,6 +171,7 @@ impl Rarray1D {
         result
     }
 
+
     /// Create vector filled with values in given range with step size
     ///
     /// # Examples
@@ -222,19 +179,19 @@ impl Rarray1D {
     /// ```
     /// use rumpy::linalg::rarray::Rarray1D;
     ///
-    /// let a = Rarray1D::range(0, 10, 2);
+    /// let a = Rarray1D::<f64>::range(0, 10, 2);
     /// println!("{}", a);
     /// // >> Rarray1D([0., 2., 4., 6., 8.])
     /// ```
-    pub fn range(start: usize, stop: usize, step: usize) -> Rarray1D {
+    pub fn range(start: usize, stop: usize, step: usize) -> Rarray1D<T> {
         let shape = (stop - start) / step;
         let mut result = Rarray1D {
             shape: D1 { height: 1, width: shape },
-            data: Vec::<f64>::with_capacity(shape)
+            data: Vec::<T>::with_capacity(shape)
         };
 
         for i in (start..stop).step_by(step) {
-            result.data.push(i as f64);
+            result.data.push(T::from(i as i32));
         }
 
         result
diff --git a/src/linalg/rarray1d_impl.rs b/src/linalg/rarray1d_impl.rs
index 75f7b8c..5821539 100644
--- a/src/linalg/rarray1d_impl.rs
+++ b/src/linalg/rarray1d_impl.rs
@@ -5,12 +5,11 @@
 
 use core::fmt;
 use std::ops::{Index, IndexMut, Mul, MulAssign};
-
-use crate::linalg::dimension::D1;
+use crate::linalg::numeric_trait::Numeric;
 use super::rarray::{Rarray, Rarray1D, RarrayCreate, RarrayScalMul};
 
-impl Index<usize> for Rarray1D {
-    type Output = f64;
+impl<T> Index<usize> for Rarray1D<T> {
+    type Output = T;
 
     fn index(&self, index: usize) -> &Self::Output {
         if self.shape[0] < self.shape[1] {
@@ -24,8 +23,8 @@ impl Index<usize> for Rarray1D {
     }
 }
 
-impl IndexMut<usize> for Rarray1D {
-    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+impl<T> IndexMut<usize> for Rarray1D<T> {
+    fn index_mut(&mut self, index: usize) -> &mut T {
         if self.shape[0] < self.shape[1] {
             assert!(index < self.shape[1], "Index out of bounds");
         } 
@@ -37,31 +36,39 @@ impl IndexMut<usize> for Rarray1D {
     }
 }
 
-impl Mul<f64> for &Rarray1D {
-    type Output = Rarray1D;
+impl<T> Mul<T> for &Rarray1D<T> where
+    T: Numeric
+{
+    type Output = Rarray1D<T>;
 
-    fn mul(self, rhs: f64) -> Self::Output {
+    fn mul(self, rhs: T) -> Self::Output {
         Rarray::scal_mul(rhs, self)
     }
 }
 
-impl MulAssign<f64> for Rarray1D {
-    fn mul_assign(&mut self, rhs: f64) {
+impl<T> MulAssign<T> for Rarray1D<T> where
+    T: Numeric
+{
+    fn mul_assign(&mut self, rhs: T) {
         self.data = Rarray::scal_mul(rhs, self).data;
     }
 }
 
-impl Mul<&Rarray1D> for &Rarray1D {
-    type Output = f64;
+impl<T> Mul<&Rarray1D<T>> for &Rarray1D<T>
+    where T: Numeric
+{
+    type Output = T;
 
-    fn mul(self, rhs: &Rarray1D) -> Self::Output {
+    fn mul(self, rhs: &Rarray1D<T>) -> Self::Output {
         assert_eq!(self.shape.width, rhs.shape.height, "Rarray shape mismatch");
         Rarray1D::dot(self, rhs)
     }
 }
 
-impl fmt::Display for Rarray1D {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+impl<T> fmt::Display for Rarray1D<T> where
+    T: Numeric
+{
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         let mut formatted_string: String = String::new();
         for i in 0..self.shape.width {
             formatted_string.push_str(&self.data[i].to_string());
diff --git a/src/linalg/rarray1d_ops.rs b/src/linalg/rarray1d_ops.rs
index 90f7a29..5e370b6 100644
--- a/src/linalg/rarray1d_ops.rs
+++ b/src/linalg/rarray1d_ops.rs
@@ -1,9 +1,9 @@
-use std::cmp::max;
-use std::iter::zip;
+use crate::linalg::numeric_trait::Numeric;
+use super::rarray::{Rarray1D, Rarray2D, D1, D2};
 
-use super::rarray::{Rarray1D, Rarray2D, RarrayCreate, D1, D2};
-
-impl Rarray1D {
+impl<T> Rarray1D<T> where 
+    T: Numeric
+{
     /// Tranpose 1D matrix
     pub fn transpose(&self) -> Self {
         Rarray1D {
@@ -18,7 +18,9 @@ impl Rarray1D {
     } 
 }
 
-impl Rarray1D {
+impl<T> Rarray1D<T> where
+    T: Numeric
+{
     /// Calculate dot product of two vectors
     ///
     /// # Examples
@@ -37,11 +39,11 @@ impl Rarray1D {
     /// # Panics
     ///
     /// If vectors aren't of same length
-    pub fn dot(a: &Self, b: &Self) -> f64 {
+    pub fn dot(a: &Self, b: &Self) -> T {
         assert_eq!(a.data.len(), b.data.len(), "Vectors not of same size");
-        let mut result: f64 = 0.;
-        for (x, y) in zip(&a.data, &b.data) {
-            result += x * y;
+        let mut result: T = T::default();
+        for i in 0..a.data.len() {
+            result += a.data[i] * b.data[i];
         }
         result
     }
@@ -54,13 +56,13 @@ impl Rarray1D {
     /// use rumpy::linalg::rarray::Rarray1D;
     /// use rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let a = Rarray1D::new(&vec![1., 1., 1.]);
-    /// let b = Rarray1D::new(&vec![1., 1., 1.]);
+    /// let a = Rarray1D::<f64>::new(&vec![1., 1., 1.]);
+    /// let b = Rarray1D::<f64>::new(&vec![1., 1., 1.]);
     /// let res = Rarray1D::outer(&a, &b);
     /// println!("{}", res);
     /// // >> Rarray2D([[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]])
     /// ```
-    pub fn outer(a: &Self, b: &Self) -> Rarray2D {
+    pub fn outer(a: &Self, b: &Self) -> Rarray2D<T> {
         let mut row = a.shape.height;
         if a.shape.width > row {
             row = a.shape.width;
@@ -71,9 +73,9 @@ impl Rarray1D {
             col = b.shape.width;
         }
 
-        let mut result: Rarray2D = Rarray2D {
+        let mut result: Rarray2D<T> = Rarray2D {
             shape: D2 { height: row, width: col },
-            data: vec![0.; row * col]
+            data: vec![T::default(); row * col]
         }; 
 
         for i in 0..row {
@@ -85,13 +87,13 @@ impl Rarray1D {
         result
     }
 
+    /*
     /// Sum values of array
-    pub fn sum(&self) -> f64 {
+    pub fn sum(&self) -> T {
         self.data.iter().sum()
     }
    
     // TODO: yet to implement functionality
-    /*
     pub fn unique(&self) -> Rarray1D {
         
     }
diff --git a/src/linalg/rarray2d_impl.rs b/src/linalg/rarray2d_impl.rs
index 7d5ed15..03a4bab 100644
--- a/src/linalg/rarray2d_impl.rs
+++ b/src/linalg/rarray2d_impl.rs
@@ -5,11 +5,11 @@
 
 use core::fmt;
 use std::ops::{Mul, MulAssign, Index, IndexMut};
-
+use crate::linalg::numeric_trait::Numeric;
 use super::rarray::{Rarray2D, RarrayCreate, RarrayMul, D2};
 
-impl Index<[usize; 2]> for Rarray2D {
-    type Output = f64;
+impl<T> Index<[usize; 2]> for Rarray2D<T> {
+    type Output = T;
 
     fn index(&self, index: [usize; 2]) -> &Self::Output {
         assert!((index[0] < self.shape[0]) && index[1] < self.shape[1], "Index out of bounds");
@@ -17,28 +17,34 @@ impl Index<[usize; 2]> for Rarray2D {
     }
 }
 
-impl IndexMut<[usize; 2]> for Rarray2D {
+impl<T> IndexMut<[usize; 2]> for Rarray2D<T> {
     fn index_mut(&mut self, index: [usize; 2]) -> &mut Self::Output {
         assert!((index[0] < self.shape[0]) && index[1] < self.shape[1], "Index out of bounds");
         &mut self.data[index[0] * self.shape[0] + index[1]]
     }
 }
 
-impl Mul<&Rarray2D> for &Rarray2D {
-    type Output = Rarray2D;
+impl<T> Mul<&Rarray2D<T>> for &Rarray2D<T> where 
+    T: Numeric
+{
+    type Output = Rarray2D<T>;
 
-    fn mul(self, rhs: &Rarray2D) -> Self::Output {
+    fn mul(self, rhs: &Rarray2D<T>) -> Self::Output {
         Rarray2D::mul(self, rhs)
     }
 }
 
-impl MulAssign<&Rarray2D> for Rarray2D {
-    fn mul_assign(&mut self, rhs: &Rarray2D) {
+impl<T> MulAssign<&Rarray2D<T>> for Rarray2D<T> where
+    T: Numeric
+{
+    fn mul_assign(&mut self, rhs: &Rarray2D<T>) {
         self.data = Rarray2D::mul(self, rhs).data;
     }
 }
 
-impl fmt::Display for Rarray2D {
+impl<T> fmt::Display for Rarray2D<T> where 
+    T: Numeric
+{
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         let mut formatted_string: String = String::new();
         for i in 0..self.shape.height {
@@ -55,7 +61,9 @@ impl fmt::Display for Rarray2D {
     }
 }
 
-impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
+impl<T> RarrayCreate<(usize, usize), Vec<Vec<T>>, T> for Rarray2D<T> where
+    T: Numeric
+{
     /// Create 2D matrix using Vec<Vec>
     ///
     /// # Examples
@@ -64,7 +72,7 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
     ///use rumpy::linalg::rarray::Rarray2D;
     ///use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    ///let m = Rarray2D::new(&vec![
+    ///let m = Rarray2D::<f64>::new(&vec![
     ///     vec![1., 1., 1.],
     ///     vec![1., 1., 1.],
     ///     vec![1., 1., 1.]
@@ -76,7 +84,7 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
     /// # Panics
     ///
     /// If not all rows are of same length
-    fn new(data: &Vec<Vec<f64>>) -> Self {
+    fn new(data: &Vec<Vec<T>>) -> Self {
         let row = data.len();
         let col = data[0].len();
         for i in 0..row {
@@ -97,40 +105,16 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
     /// use rumpy::linalg::rarray::Rarray2D;
     /// use crate::rumpy::linalg::rarray::RarrayCreate;
     ///
-    /// let m = Rarray2D::zeros((3, 3));
+    /// let m = Rarray2D::<f64>::zeros((3, 3));
     /// println!("{}", m);
     /// // >> Rarray2D([0., 0., 0.], [0., 0., 0.], [0., 0., 0.])
     /// ```
     fn zeros(shape: (usize, usize)) -> Self {
         Rarray2D {
             shape: D2 { height: shape.0, width: shape.1 },
-            data: vec![0.; shape.0 * shape.1]
+            data: vec![T::default(); shape.0 * shape.1]
         } 
     }
-
-    /// Create 2D matrix of shape `m x n` filled with random numbers
-    ///
-    /// # Examples
-    ///
-    /// ```
-    /// use rumpy::linalg::rarray::Rarray2D;
-    /// use crate::rumpy::linalg::rarray::RarrayCreate;
-    ///
-    /// let m = Rarray2D::random((3, 3));
-    /// println!("{}", m);
-    /// // >> Rarray2D([])
-    /// ```
-    fn random(shape: (usize, usize)) -> Self {
-        let mut data = Vec::<f64>::with_capacity(shape.0 * shape.1);
-        for _ in 0..(shape.0 * shape.1) {
-            data.push(rand::random::<f64>());
-        }
-
-        Rarray2D {
-            shape: D2 { height: shape.0, width: shape.1 },
-            data
-        }
-    }
     
     /// Create 2D matrix of shape `m x n` filled with `value`
     ///
@@ -144,7 +128,7 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
     /// println!("{}", m);
     /// // >> Rarray2D([[4., 4., 4.], [4., 4., 4.], [4., 4., 4.]])
     /// ```
-    fn fill(value: f64, shape: (usize, usize)) -> Self {
+    fn fill(value: T, shape: (usize, usize)) -> Self {
         Rarray2D {
             shape: D2 { height: shape.0, width: shape.1 },
             data: vec![value; shape.0 * shape.1]
@@ -152,7 +136,9 @@ impl RarrayCreate<(usize, usize), Vec<Vec<f64>>, f64> for Rarray2D {
     }
 }
 
-impl Rarray2D {
+impl<T> Rarray2D<T> where
+    T: Numeric
+{
     /// Create 2D matrix of shape `m x n` filled with `value`
     ///
     /// # Examples
@@ -160,18 +146,18 @@ impl Rarray2D {
     /// ```
     /// use rumpy::linalg::rarray::Rarray2D;
     ///
-    /// let m = Rarray2D::ones(3);
+    /// let m = Rarray2D::<f64>::ones(3);
     /// println!("{}", m);
     /// // >> Rarray2D([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
     /// ```
     pub fn ones(shape: usize) -> Self {
         let mut result = Rarray2D {
             shape: D2 { height: shape, width: shape },
-            data: vec![0.; shape * shape]
+            data: vec![T::default(); shape * shape]
         };
 
         for i in 0..shape {
-            result.data[shape * i + i] = 1.
+            result.data[shape * i + i] = T::from(1);
         }
 
         result
diff --git a/src/linalg/rarray_impl.rs b/src/linalg/rarray_impl.rs
index 0a5a3a2..db39955 100644
--- a/src/linalg/rarray_impl.rs
+++ b/src/linalg/rarray_impl.rs
@@ -1,23 +1,22 @@
 use super::rarray::{Rarray, Rarray1D, RarrayAdd, RarraySub};
 use std::{ops::{Add, AddAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign}, usize};
 use std::fmt::Debug;
-use std::process::Output;
+use std::ops::Neg;
 use crate::linalg::dimension::Dim;
+use crate::linalg::numeric_trait::Numeric;
 
-impl<T, D> Rarray<T, D> {
+impl<T, D> Rarray<T, D> where
+    D: Dim
+{
     /// Return shape of array
     pub fn get_shape(&self) -> &D {
         &self.shape
     }
-    // TODO: Implement following functions
-
-    /// Move ownership of array
-    pub fn to_owned(self) {}
 }
 
 // Base operations for the Rarray abstract struct
 impl<T, D> Add<&Rarray<T, D>> for &Rarray<T, D> where
-    T : Add<Output = T> + Copy + Default,
+    T : Numeric,
     D : Copy + Dim + Debug + Eq
 {
     type Output = Rarray<T, D>;
@@ -28,7 +27,7 @@ impl<T, D> Add<&Rarray<T, D>> for &Rarray<T, D> where
 }
 
 impl<T, D> AddAssign<&Rarray<T, D>> for Rarray<T, D> where 
-    T: Add<Output = T> + Copy + Default,
+    T: Numeric,
     D: Copy + Dim + Debug + Eq
 {
     fn add_assign(&mut self, rhs: &Rarray<T, D>) {
@@ -38,7 +37,7 @@ impl<T, D> AddAssign<&Rarray<T, D>> for Rarray<T, D> where
 
 
 impl<T, D> Sub<&Rarray<T, D>> for &Rarray<T, D> where 
-    T : Sub<Output = T> + Copy + Default,
+    T : Numeric,
     D : Copy + Dim + Debug + Eq
 {
     type Output = Rarray<T, D>;
@@ -50,7 +49,7 @@ impl<T, D> Sub<&Rarray<T, D>> for &Rarray<T, D> where
 
 
 impl<T, D> SubAssign<&Rarray<T, D>> for Rarray<T, D> where
-    T : Sub<Output = T> + Copy + Default,
+    T : Numeric,
     D : Copy + Dim + Debug + Eq
 {
     fn sub_assign(&mut self, rhs: &Rarray<T, D>) {
diff --git a/tests/rarray1d_create_test.rs b/tests/rarray1d_create_test.rs
index 2cf6f3b..aa2d18f 100644
--- a/tests/rarray1d_create_test.rs
+++ b/tests/rarray1d_create_test.rs
@@ -1,16 +1,13 @@
 mod test {
-    use std::alloc::handle_alloc_error;
-    use std::env::var;
     use rstest::rstest;
-    use rumpy::linalg::rarray;
-    use rumpy::linalg::rarray::{Rarray1D, Rarray2D, RarrayCreate};
+    use rumpy::linalg::rarray::{Rarray1D, RarrayCreate};
 
     #[rstest]
     #[case(vec![1., 1., 1.])]
     #[case(vec![0.])]
     #[case(vec![2.4, 2., 1., 0.4, 6.,])]
     fn rarray1d_new_hor(#[case] a: Vec<f64>){
-        let rarray_a = Rarray1D::new(&a);
+        let rarray_a = Rarray1D::<f64>::new(&a);
         for i in 0..a.len() {
             assert_eq!(a[i], rarray_a[i]);
         }
@@ -21,7 +18,7 @@ mod test {
     #[case(vec![vec![0.]])]
     #[case(vec![vec![2.4], vec![2.], vec![1.], vec![0.4], vec![6.]])]
     fn rarray1d_new_ver(#[case] a: Vec<Vec<f64>>){
-        let rarray_a = Rarray1D::new(&a);
+        let rarray_a = Rarray1D::<f64>::new(&a);
         for i in 0..a[0].len() {
             assert_eq!(a[0][i], rarray_a[i]);
         }
@@ -32,7 +29,7 @@ mod test {
     #[case(1)]
     #[case(10)]
     fn rarray1d_zeros_hor(#[case] shape: usize){
-        let rarray_zeros = Rarray1D::zeros(shape);
+        let rarray_zeros = Rarray1D::<f64>::zeros(shape);
         for i in 0..shape {
             assert_eq!(rarray_zeros[i], 0.);
         }
@@ -43,7 +40,7 @@ mod test {
     #[case(1, 1)]
     #[case(10, 1)]
     fn rarray1d_zeros_ver(#[case] height: usize, #[case] width: usize){
-        let rarray_zeros = Rarray1D::zeros((height, width));
+        let rarray_zeros = Rarray1D::<f64>::zeros((height, width));
         for i in 0..height {
             assert_eq!(rarray_zeros[i], 0.);
         }