diff --git a/src/primary/covariance.rs b/src/primary/covariance.rs
new file mode 100644
index 0000000..29c9041
--- /dev/null
+++ b/src/primary/covariance.rs
@@ -0,0 +1,37 @@
+/// Returns the arithmetic mean of a sequence of `f64` values.
+///
+/// Empty input returns `None`.
+///
+/// # Example
+/// ```
+/// use rustats::primary::mean;
+///
+/// let values_x = [1.0, 2.0, 3.0];
+/// let values_y = [4.0, 5.0, 6.0];
+/// let cv = covariance(&values_x, &values_y);
+/// assert_eq!(cv, Some(1.0));
+/// ```
+use super::mean;
+
+pub fn covariance(x: &[f64], y: &[f64]) -> Option<f64> {
+    let mut denominator: f64 = 0.0;
+
+    if x.len() != y.len(){
+        return None;
+    }
+
+    if x.len() < 2 {
+            return None;
+    }
+
+    let mean_x = mean(x)?;
+    let mean_y = mean(y)?;
+
+        for i in 0..x.len() {
+            denominator += (x[i]-mean_x)*(y[i]-mean_y);
+        }
+        Some(denominator / (x.len() - 1) as f64)
+
+}
+    
+
diff --git a/src/primary/mod.rs b/src/primary/mod.rs
index 9ee4243..2ea7492 100644
--- a/src/primary/mod.rs
+++ b/src/primary/mod.rs
@@ -6,6 +6,7 @@ pub mod variance;
 pub mod std_deviation;
 pub mod interquartile_range;
 pub mod range;
+pub mod covariance;
 
 
 pub use mean::mean;
@@ -14,5 +15,6 @@ pub use mode::mode;
 pub use quartiles::{Quartiles, quartiles};
 pub use variance::variance;
 pub use std_deviation::std_deviation;
+pub use covariance::covariance;
+pub use range::range;
 pub use interquartile_range::interquartile_range;
-pub use range::range;
\ No newline at end of file
diff --git a/tests/covariance.rs b/tests/covariance.rs
new file mode 100644
index 0000000..fc14fff
--- /dev/null
+++ b/tests/covariance.rs
@@ -0,0 +1,96 @@
+use rustats::primary::covariance;
+
+fn approx_eq(left: f64, right: f64, epsilon: f64) -> bool {
+    (left - right).abs() <= epsilon
+}
+
+#[test]
+fn covariance_empty_returns_none() {
+    // Covariance requires at least 2 values in each input (denominator is n-1)
+    assert_eq!(covariance(&[], &[]), None);
+}
+
+#[test]
+fn covariance_different_lengths(){
+    //Covariance requires both input arrays to have the same length
+    let values_x = [1.0, 2.0, 3.0];
+    let values_y = [2.0, 3.0];
+    assert_eq!(covariance(&values_x, &values_y),None);
+}
+
+#[test]
+fn covariance_single_element_arrays_returns_none() {
+    // Covariance requires at least 2 values in each input (denominator is n-1)
+    assert_eq!(covariance(&[5.0], &[9.0]), None);
+}
+
+#[test]
+fn covariance_two_identical_sample_arrays() {
+    // If both input arrays are identical, the covariance should be equal to the variance of either array, which is 0 since all values are the same.
+    let values_x = [4.0, 4.0];
+    let values_y = [4.0, 4.0];
+    let cv = covariance(&values_x, &values_y).expect("two or more values should return Some");
+    assert!(approx_eq(cv, 0.0, 1e-12));
+}
+
+#[test]
+fn covariance_smallest_valid_input() {
+    // mean_x = 1.5, mean_y = 3.5, denominator= (-0.5*-0.5)+(0.5*0.5) = 0.5, covariance = denominator / (2-1) = 0.5 / 1 = 0.5
+    let values_x = [1.0, 2.0];
+    let values_y = [3.0, 4.0];
+    let cv = covariance(&values_x, &values_y).expect("two or more values should return Some");
+    assert!(approx_eq(cv, 0.5, 1e-12));
+}
+
+#[test]
+fn covariance_multiple_values() {
+    // mean_x = 3.0, mean_y = 7.0, deviations_x = [-2.0, -1.0, 0, 1.0, 2.0], deviations_y = [-2.0, -1.0, 0, 1.0, 2.0, denominator = (-2*-2) + (-1*-1) + (0*0) + (1*1) + (2*2) = 10.0 
+    //covariance = denominator / (5-1) = 2.5
+    let values_x = [1.0, 2.0, 3.0, 4.0, 5.0];
+    let values_y = [5.0, 6.0, 7.0, 8.0, 9.0];
+    let cv = covariance(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(cv, 2.5, 1e-12));
+}
+
+#[test]
+fn covariance_with_negative_values() {
+    // mean_x = 0.0, mean_y = 0.0, deviations_x = [-2, -1, 1, 2], deviations_y = [-5.0, -4.0, 4.0, 5.0]  denominator = (-2*-5)+(-1*-4)+(1*4)+(2*5) = 28.0 
+    //covariance = 28 / 3 ≈ 9.333...
+    let values_x = [-2.0, -1.0, 1.0, 2.0];
+    let values_y = [-5.0, -4.0, 4.0, 5.0];
+    let cv = covariance(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(cv, 28.0 / 3.0, 1e-12));
+}
+
+#[test]
+fn covariance_with_decimals() {
+    // mean_x = 2.0, mean_y = 2.0, deviations_x = [-1.5, -0.5, 0.5, 1.5], deviations_y = [-1.5, -0.5, 0.5, 1.5], denominator = (-1.5*-1.5)+(-0.5*-0.5)+(0.5*0.5)+(1.5*1.5) = 5.0
+    //covariance = 5 / 3 ≈ 1.6667
+    let values_x = [0.5, 1.5, 2.5, 3.5];
+    let values_y = [0.5, 1.5, 2.5, 3.5];
+    let v = covariance(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(v, 5.0 / 3.0, 1e-12));
+}
+
+#[test]
+fn covariance_negative() {
+    // covariance is negative when the deviations of x and y have opposite signs, which means that when x is above its mean, y tends to be below its mean, and vice versa. In this case, the covariance should be negative.
+    // mean_x = 0.0, mean_y = 0.0, deviations_x = [-10.0, 0.0, 10.0, 5.0, -5.0], deviations_y = [10.0, 0.0, -10.0, -5.0, 5.0], denominator = (-10*10)+(0*0)+(10*-10)+(5*-5)+(-5*5) = -250.0
+    //covariance = -250 / 4 = -62.5
+    let values_x = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let values_y = [10.0, 0.0, -10.0, -5.0, 5.0];
+    let cv = covariance(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(cv < 0.0);
+}
+
+#[test]
+fn covariance_positive() {
+    // covariance is positive when the deviations of x and y have the same signs, which means that when x is above its mean, y tends to be above its mean, and vice versa. In this case, the covariance should be positive.
+    // mean_x = 0.0, mean_y = 0.0, deviations_x = [-10.0, 0.0, 10.0, 5.0, -5.0], deviations_y = [-10.0, 0.0, 10.0, 5.0, -5.0], denominator = (-10*-10)+(0*0)+(10*10)+(5*5)+(-5*-5) = 250.0
+    //covariance = 250 / 4 = 62.5
+    let values_x = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let values_y = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let cv = covariance(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(cv > 0.0);
+
+}
\ No newline at end of file