diff --git a/src/bindings/PyDP/algorithms/util.cpp b/src/bindings/PyDP/algorithms/util.cpp
index ae81dc27..796d625d 100644
--- a/src/bindings/PyDP/algorithms/util.cpp
+++ b/src/bindings/PyDP/algorithms/util.cpp
@@ -1,6 +1,7 @@
 // Provides bindings for Util
 
 #include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
 
 #include "differential_privacy/algorithms/util.h"
 
@@ -14,4 +15,31 @@ void init_algorithms_util(py::module& m) {
   util.def("default_epsilon", &dp::DefaultEpsilon);
   util.def("get_next_power_of_two", &dp::GetNextPowerOfTwo);
   util.def("qnorm", &dp::Qnorm);
-}
+  util.def("mean", &dp::Mean<double>);
+  util.def("mean", &dp::Mean<int>);
+  util.def("variance", &dp::Variance<double>);
+  util.def("standard_deviation", &dp::StandardDev<double>);
+  util.def("order_statistics", &dp::OrderStatistic<double>);
+  util.def("correlation", &dp::Correlation<double>);
+  util.def("vector_filter", &dp::VectorFilter<double>);
+  util.def("vector_to_string", &dp::VectorToString<double>);
+  util.def("round_to_nearest_multiple", &dp::RoundToNearestMultiple);
+  util.def("safe_add", [](int64_t i, int64_t j) {
+    int64_t k;
+    bool result = dp::SafeAdd(i, j, &k);
+    if (result) return k;
+    throw std::runtime_error("Result of addition will overflow.");
+  });
+  util.def("safe_subtract", [](int64_t i, int64_t j) {
+    int64_t k;
+    bool result = dp::SafeSubtract(i, j, &k);
+    if (result) return k;
+    throw std::runtime_error("Result of subtraction will overflow.");
+  });
+  util.def("safe_square", [](int64_t i) {
+    int64_t k;
+    bool result = dp::SafeSquare(i, &k);
+    if (result) return k;
+    throw std::runtime_error("Result of squaring will overflow.");
+  });
+}
\ No newline at end of file
diff --git a/tests/algorithms/test_util.py b/tests/algorithms/test_util.py
new file mode 100644
index 00000000..a6514797
--- /dev/null
+++ b/tests/algorithms/test_util.py
@@ -0,0 +1,98 @@
+import pytest
+import pydp as dp
+import math
+
+
+def test_default_epsilon():
+    assert dp.util.default_epsilon() == math.log(3)
+
+
+def test_next_power_positive():
+    kTolerance = 1e-5
+    npp1 = dp.util.get_next_power_of_two(3.0)
+    npp2 = dp.util.get_next_power_of_two(5.0)
+    npp3 = dp.util.get_next_power_of_two(7.9)
+    assert abs(npp1 - 4) < kTolerance
+    assert abs(npp2 - 8) < kTolerance
+    assert abs(npp3 - 8) < kTolerance
+
+
+def test_next_power_exact_positive():
+    kTolerance = 1e-5
+    npep1 = dp.util.get_next_power_of_two(2.0)
+    npep2 = dp.util.get_next_power_of_two(8.0)
+    assert abs(npep1 - 2) < kTolerance
+    assert abs(npep2 - 8) < kTolerance
+
+
+def test_next_power_one():
+    kTolerance = 1e-5
+    npo = dp.util.get_next_power_of_two(1.0)
+    assert abs(npo - 1) < kTolerance
+
+
+def test_next_power_negative():
+    kTolerance = 1e-5
+    npn1 = dp.util.get_next_power_of_two(0.4)
+    npn2 = dp.util.get_next_power_of_two(0.2)
+    assert abs(npn1 - 0.5) < kTolerance
+    assert abs(npn2 - 0.25) < kTolerance
+
+
+def test_next_power_exact_negative():
+    kTolerance = 1e-5
+    npn1 = dp.util.get_next_power_of_two(0.5)
+    npn2 = dp.util.get_next_power_of_two(0.125)
+    assert abs(npn1 - 0.5) < kTolerance
+    assert abs(npn2 - 0.125) < kTolerance
+
+
+def test_round_positive():
+    kTolerance = 1e-5
+    rp1 = dp.util.round_to_nearest_multiple(4.9, 2.0)
+    rp2 = dp.util.round_to_nearest_multiple(5.1, 2.0)
+    assert abs(rp1 - 4) < kTolerance
+    assert abs(rp2 - 6) < kTolerance
+
+
+def test_round_negative():
+    kTolerance = 1e-5
+    rn1 = dp.util.round_to_nearest_multiple(-4.9, 2.0)
+    rn2 = dp.util.round_to_nearest_multiple(-5.1, 2.0)
+    assert abs(rn1 + 4) < kTolerance
+    assert abs(rn2 + 6) < kTolerance
+
+
+def test_round_positive_ties():
+    kTolerance = 1e-5
+    rpt = dp.util.round_to_nearest_multiple(5.0, 2.0)
+    assert abs(rpt - 6.0) < kTolerance
+
+
+def test_round_negative_ties():
+    kTolerance = 1e-5
+    rnt = dp.util.round_to_nearest_multiple(-5.0, 2.0)
+    assert abs(rnt + 4.0) < kTolerance
+
+
+def test_statistics():
+    a = [1.0, 5.0, 7.0, 9.0, 13.0]
+    assert dp.util.mean(a) == 7.0
+    assert dp.util.variance(a) == 16.0
+    assert dp.util.standard_deviation(a) == 4.0
+    assert dp.util.order_statistics(0.6, a) == 8.0
+    assert dp.util.order_statistics(0, a) == 1.0
+    assert dp.util.order_statistics(1, a) == 13.0
+
+
+def test_vector_filter():
+    v = [1.0, 2.0, 2.0, 3.0]
+    selection = [False, True, True, False]
+    expected = [2.0, 2.0]
+    assert expected == dp.util.vector_filter(v, selection)
+
+
+def test_vector_to_string():
+    v = [1.0, 2.0, 2.0, 3.0]
+    expected = "[1, 2, 2, 3]"
+    assert dp.util.vector_to_string(v) == expected