diff --git a/AutomaticDifferentiation.hpp b/AutomaticDifferentiation.hpp
index a939ddb..c104968 100755
--- a/AutomaticDifferentiation.hpp
+++ b/AutomaticDifferentiation.hpp
@@ -7,13 +7,13 @@
 
 #include <valarray>
 
-template<typename T>
-std::ostream & operator<<(std::ostream & out, std::valarray<T> const& v)
-{
-	for(size_t i = 0 ; i < v.size() ; i++)
-		out << v[i] << " ";
-	return out;
-}
+// template<typename T>
+// std::ostream & operator<<(std::ostream & out, std::valarray<T> const& v)
+// {
+// 	for(size_t i = 0 ; i < v.size() ; i++)
+// 		out << v[i] << " ";
+// 	return out;
+// }
 
 #define MINAB(a, b) (((a) < (b)) ? (a) : (b))
 
@@ -64,7 +64,7 @@ class Dual
 				// copy old data into new b vector
 				VectorT b_old = b;
 				b.resize(N);
-				for(size_t j = 0 ; j < MINAB(b.size(), b_old.size()) ; j++)
+				for(size_t j = 0 ; j < min(b.size(), b_old.size()) ; j++)
 					b[j] = b_old[j];
 			}
 			b[i] = Scalar(1.);
@@ -433,6 +433,10 @@ template<typename Scalar> Dual<Scalar> sqrt(const Dual<Scalar> & x) {
 	return Dual<Scalar>(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a)));
 }
 
+template<typename Scalar> Dual<Scalar> cbrt(const Dual<Scalar> & x) {
+	return Dual<Scalar>(cbrt(x.a), x.b/(Scalar(3.)*pow(x.a, Scalar(2.)/Scalar(3.))));
+}
+
 template<typename Scalar> Dual<Scalar> sign(const Dual<Scalar> & x) {
 	return Dual<Scalar>(sign(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
 }
diff --git a/test/Makefile b/test/Makefile
new file mode 100755
index 0000000..a1efaa2
--- /dev/null
+++ b/test/Makefile
@@ -0,0 +1,37 @@
+# Declaration of variables
+C = clang
+COMMON_FLAGS = -Wall -MMD
+C_FLAGS = $(COMMON_FLAGS)
+CC = clang++
+CC_FLAGS = $(COMMON_FLAGS) -std=c++17 -g -O0#-O2
+LD_FLAGS = -lgmp -lmpfr
+INCLUDES = -I../../Catch2-master/single_include
+
+# File names
+EXEC = run
+CSOURCES = $(wildcard *.c)
+COBJECTS = $(CSOURCES:.c=.o)
+SOURCES = $(wildcard *.cpp)
+OBJECTS = $(SOURCES:.cpp=.o)
+
+# Main target
+$(EXEC): $(COBJECTS) $(OBJECTS)
+	$(CC) $(COBJECTS) $(OBJECTS) -o $(EXEC) $(LD_FLAGS)
+
+# To obtain object files
+%.o: %.cpp
+	$(CC) $(INCLUDES) $(CC_FLAGS) -o $@ -c $<
+
+# To obtain object files
+%.o: %.c
+	$(C) $(INCLUDES) $(C_FLAGS) -o $@ -c $<
+
+-include $(SOURCES:%.cpp=%.d)
+-include $(CSOURCES:%.c=%.d)
+
+# To remove generated files
+clean:
+	rm -f $(COBJECTS) $(OBJECTS) $(SOURCES:%.cpp=%.d) $(CSOURCES:%.c=%.d)
+
+cleaner: clean
+	rm -f $(EXEC)
diff --git a/test/dual_test.cpp b/test/dual_test.cpp
new file mode 100755
index 0000000..f73803a
--- /dev/null
+++ b/test/dual_test.cpp
@@ -0,0 +1,40 @@
+#include <catch2/catch.hpp>
+#include <iostream>
+
+#include "../utils.hpp"
+#include "../AutomaticDifferentiation.hpp"
+
+#define print(x) PRINT_VAR(x);
+#define printstr(x) PRINT_STR(x);
+
+template<typename T>
+bool check_almost_equal(const std::vector<T> & v1, const std::vector<T> & v2, T tol)
+{
+    bool ok = true;
+    for(size_t i = 0 ; i < v1.size() ; i++)
+        if(fabs(v1[i] - v2[i]) > tol)
+        {
+            ok = false;
+            break;
+        }
+    return ok;
+}
+
+TEST_CASE( "Dual numbers : Constructors", "[dual]" )
+{
+    std::cout.precision(16);
+    using S = double;
+    using D = Dual<S>;
+    using V = std::vector<S>;
+
+    SECTION( "Default Constructor" ) {
+
+    }
+
+    SECTION( "Constructor with int" ) {
+
+    }
+}
+
+// CHECK()
+// REQUIRE()
diff --git a/test/main_test.cpp b/test/main_test.cpp
new file mode 100755
index 0000000..0465d9b
--- /dev/null
+++ b/test/main_test.cpp
@@ -0,0 +1,4 @@
+// Let Catch provide main()
+#define CATCH_CONFIG_MAIN
+
+#include <catch2/catch.hpp>
diff --git a/tests/compile_time_derivative/main.cpp b/tests/compile_time_derivative/main.cpp
index daa535f..d5db6b7 100755
--- a/tests/compile_time_derivative/main.cpp
+++ b/tests/compile_time_derivative/main.cpp
@@ -21,7 +21,7 @@ T dfct(const T & x) {
 	return (6)*x + (2) - (1)/(x*x);
 }
 
-// The file must be compiled with -O0 because with -O2 and up
+// The file must be compiled with -O0 because with -O1 and up the compiler simplifies func and dfunc out...
 int main()
 {
 	auto func = (3*x*x + 2*x - 3 + 1./x);
diff --git a/utils.hpp b/utils.hpp
new file mode 100755
index 0000000..b419ae6
--- /dev/null
+++ b/utils.hpp
@@ -0,0 +1,35 @@
+#ifndef DEF_utils
+#define DEF_utils
+
+#include <ostream>
+
+#define PRINT_VAR(x) std::cout << #x << "\t= " << (x) << "\n"
+#define PRINT_VEC(x);	{std::cout << #x << "\t= "; \
+						for(unsigned int i_print_vec = 0 ; i_print_vec < (x).size() ; i_print_vec++) \
+							std::cout << (x)[i_print_vec] << "\t"; \
+						std::cout << "\n";}
+#define PRINT_STR(x) std::cout << (x) << "\n"
+
+template<typename T, typename T2> auto min(const T & a, const T2 & b) -> decltype(a | b) { return (a < b) ? a : b; }
+template<typename T, typename T2> auto max(const T & a, const T2 & b) -> decltype(a | b) { return (a < b) ? b : a; }
+
+template<typename T, template<class,class...> class C, class... Args>
+std::ostream& operator<<(std::ostream& os, const C<T,Args...>& objs)
+{
+    for (auto const& obj : objs)
+        os << obj << ' ';
+    return os;
+}
+
+/// Convenience template class to do StdPairValueCatcher(a, b) = someFunctionThatReturnsAnStdPair<A,B>()
+/// Instead of doing a = std::get<0>(result_from_function), b = std::get<1>(result_from_function)
+template<typename A, typename B>
+struct StdPairValueCatcher
+{
+    A & a;
+    B & b;
+    StdPairValueCatcher(A & a_, B & b_) : a(a_), b(b_) {}
+    std::pair<A, B> const& operator=(std::pair<A, B> const& p) {a = std::get<0>(p);  b = std::get<1>(p); return p; }
+};
+
+#endif