diff --git a/src/ForwardModeAD.chpl b/src/ForwardModeAD.chpl
index 6dac733..b1d2449 100644
--- a/src/ForwardModeAD.chpl
+++ b/src/ForwardModeAD.chpl
@@ -1,16 +1,9 @@
 module ForwardModeAD {
+  include module DualType;
+  include module ElementaryFunctions;
+  include module Differentiation;
 
-  public use Math;
-
-  public use dualtype;
-  
-  public use arithmetic;
-
-  public use trigonometric;
-
-  public use transcendental;
-
-  public use hyperbolic;
-
-  public use differentiation;
+  public use DualType;
+  public use ElementaryFunctions;
+  public use Differentiation;
 }
diff --git a/src/differentiation.chpl b/src/ForwardModeAD/Differentiation.chpl
similarity index 98%
rename from src/differentiation.chpl
rename to src/ForwardModeAD/Differentiation.chpl
index 252b34e..d1fdade 100644
--- a/src/differentiation.chpl
+++ b/src/ForwardModeAD/Differentiation.chpl
@@ -1,6 +1,5 @@
-module differentiation {
-  
-  use ForwardModeAD;
+module Differentiation {
+  use ForwardModeAD.DualType;
 
   /*
   Initializes the input to the appropriate dual number to evalute the derivative.
diff --git a/src/dualtype.chpl b/src/ForwardModeAD/DualType.chpl
similarity index 99%
rename from src/dualtype.chpl
rename to src/ForwardModeAD/DualType.chpl
index a2c6f91..9887876 100644
--- a/src/dualtype.chpl
+++ b/src/ForwardModeAD/DualType.chpl
@@ -1,4 +1,4 @@
-module dualtype {  
+module DualType {  
 
   /*
   A dual number is a number in the form :math:`a + b\epsilon`, for which :math:`\epsilon^2 = 0`.
diff --git a/src/ForwardModeAD/ElementaryFunctions.chpl b/src/ForwardModeAD/ElementaryFunctions.chpl
new file mode 100644
index 0000000..48bcaf7
--- /dev/null
+++ b/src/ForwardModeAD/ElementaryFunctions.chpl
@@ -0,0 +1,184 @@
+module ElementaryFunctions {
+    use ForwardModeAD.DualType;
+    use Math;
+
+    // Arithmetic operations
+    operator +(a) where isDualType(a.type) { return a; }
+
+    operator -(a) where isDualType(a.type) { 
+        var f = -primalPart(a),
+            df = -dualPart(a);
+        return todual(f, df); 
+    }
+
+    operator +(a, b) where isEitherDualType(a.type, b.type) {
+        var f = primalPart(a) + primalPart(b);
+        var df = dualPart(a) + dualPart(b);
+        return todual(f, df);
+    }
+
+    operator -(a, b) where isEitherDualType(a.type, b.type) {
+        var f = primalPart(a) - primalPart(b);
+        var df = dualPart(a) - dualPart(b);
+        return todual(f, df);
+    }
+
+    operator *(a, b) where isEitherDualType(a.type, b.type) {
+        var f = primalPart(a) * primalPart(b),
+            df = dualPart(a) * primalPart(b) + primalPart(a) * dualPart(b);
+        return todual(f, df);
+    }
+
+    operator /(a, b) where isEitherDualType(a.type, b.type) {
+        var f = primalPart(a) / primalPart(b),
+            df = (dualPart(a) * primalPart(b) - primalPart(a) * dualPart(b)) / primalPart(b) ** 2;
+        return todual(f, df); 
+    }
+
+    operator **(a, b : real) where isDualType(a.type) {
+        var f = primalPart(a) ** b,
+            df = b * (primalPart(a) ** (b - 1)) * dualPart(a);
+        return todual(f, df);
+    }
+
+    proc sqrt(a) where isDualType(a.type) {
+        var f = sqrt(primalPart(a)),
+            df = 0.5 * dualPart(a) / sqrt(primalPart(a));
+        return todual(f, df);
+    }
+
+    proc cbrt(a) where isDualType(a.type) {
+        var f = cbrt(primalPart(a)),
+            df =  1.0 / 3.0 * dualPart(a) / cbrt(primalPart(a) ** 2);
+        return todual(f, df);
+    }
+
+    // Trigonometric functions
+    proc sin(a) where isDualType(a.type) {
+        var f = sin(primalPart(a)),
+            df = cos(primalPart(a)) * dualPart(a);
+        return todual(f, df);
+    }
+
+    proc cos(a) where isDualType(a.type) {
+        var f = cos(primalPart(a)),
+            df = -sin(primalPart(a)) * dualPart(a);
+        return todual(f, df);
+    }
+
+    proc tan(a) where isDualType(a.type) {
+        var f = tan(primalPart(a)),
+            df = dualPart(a) / (cos(primalPart(a)) ** 2);
+        return todual(f, df);
+    }
+
+    proc asin(a) where isDualType(a.type) {
+        var f = asin(primalPart(a)),
+            df = dualPart(a) / sqrt(1 - primalPart(a)**2);
+        return todual(f, df);
+    }
+
+    proc acos(a) where isDualType(a.type) {
+        var f = acos(primalPart(a)),
+            df = -dualPart(a) / sqrt(1 - primalPart(a)**2);
+        return todual(f, df);
+    }
+
+    proc atan(a) where isDualType(a.type) {
+        var f = atan(primalPart(a)),
+            df = dualPart(a) / (1 + primalPart(a)**2);
+        return todual(f, df);
+    }
+
+    // Trascendental functions
+    operator **(a : real, b) where isDualType(b.type) {
+        var f = a ** primalPart(b),
+            df = log(a) * (a ** primalPart(b)) * dualPart(b);
+        return todual(f, df);
+    }
+
+    operator **(a, b) where isDualType(a.type) && a.type == b.type {
+        var f = primalPart(a) ** primalPart(b);
+        var df = f * (dualPart(b) * log(primalPart(a)) + primalPart(b) * dualPart(a) / primalPart(a)); 
+        return todual(f, df);
+    }
+
+    proc exp(a) where isDualType(a.type) {
+        var f = exp(primalPart(a)),
+            df = dualPart(a) * exp(primalPart(a));
+        return todual(f, df);
+    }
+
+    proc exp2(a) where isDualType(a.type) {
+        var f = exp2(primalPart(a)), 
+            df = ln2 * exp2(primalPart(a)) * dualPart(a);
+        return todual(f, df);
+    }
+
+    proc expm1(a) where isDualType(a.type) {
+        var f = expm1(primalPart(a)), 
+            df = exp(primalPart(a)) * dualPart(a);
+        return todual(f, df);
+    }
+
+    proc log(a) where isDualType(a.type) {
+        var f = log(primalPart(a)), 
+            df = dualPart(a) / primalPart(a);
+        return todual(f, df);
+    }
+    
+    proc log2(a) where isDualType(a.type) {
+        var f = log2(primalPart(a)), 
+            df = dualPart(a)/(primalPart(a) * ln2);
+        return todual(f, df);
+    }
+
+    proc log10(a) where isDualType(a.type) {
+        var f = log10(primalPart(a)), 
+            df = dualPart(a) / (primalPart(a) * ln10);
+        return todual(f, df);
+    }
+
+    proc log1p(a) where isDualType(a.type) {
+        var f = log1p(primalPart(a)), 
+            df = dualPart(a) / (primalPart(a) + 1);
+        return todual(f, df);
+    }
+
+    // Hyperbolic functions
+    proc sinh(a) where isDualType(a.type) {
+        var f = sinh(primalPart(a)), 
+        df = dualPart(a) * cosh(primalPart(a));
+        return todual(f, df);
+    }
+
+    proc cosh(a) where isDualType(a.type) {
+        var f = cosh(primalPart(a)), 
+            df = dualPart(a) * sinh(primalPart(a));
+        return todual(f, df);
+    }
+
+    proc tanh(a) where isDualType(a.type) {
+        var f = tanh(primalPart(a)), 
+            df = dualPart(a) * (1 - tanh(primalPart(a))**2);
+        return todual(f, df);
+    }
+
+    proc asinh(a) where isDualType(a.type) {
+        var f = asinh(primalPart(a)), 
+            df = dualPart(a) / sqrt(primalPart(a)**2 + 1);
+        return todual(f, df);
+    }
+
+    proc acosh(a) where isDualType(a.type) {
+        var f = acosh(primalPart(a)), 
+            df = dualPart(a) / sqrt(primalPart(a)**2 - 1);
+        return todual(f, df);
+    }
+
+    proc atanh(a) where isDualType(a.type) {
+        var f = atanh(primalPart(a)), 
+            df = dualPart(a) / (1 - primalPart(a) ** 2);
+        return todual(f, df);
+    }
+}
diff --git a/src/arithmetic.chpl b/src/arithmetic.chpl
deleted file mode 100644
index 28be41c..0000000
--- a/src/arithmetic.chpl
+++ /dev/null
@@ -1,53 +0,0 @@
-module arithmetic {
-    use ForwardModeAD;
-    
-    operator +(a) where isDualType(a.type) { return a; }
-
-    operator -(a) where isDualType(a.type) { 
-        var f = -primalPart(a),
-            df = -dualPart(a);
-        return todual(f, df); 
-    }
-
-    operator +(a, b) where isEitherDualType(a.type, b.type) {
-        var f = primalPart(a) + primalPart(b);
-        var df = dualPart(a) + dualPart(b);
-        return todual(f, df);
-    }
-
-    operator -(a, b) where isEitherDualType(a.type, b.type) {
-        var f = primalPart(a) - primalPart(b);
-        var df = dualPart(a) - dualPart(b);
-        return todual(f, df);
-    }
-
-    operator *(a, b) where isEitherDualType(a.type, b.type) {
-        var f = primalPart(a) * primalPart(b),
-            df = dualPart(a) * primalPart(b) + primalPart(a) * dualPart(b);
-        return todual(f, df);
-    }
-
-    operator /(a, b) where isEitherDualType(a.type, b.type) {
-        var f = primalPart(a) / primalPart(b),
-            df = (dualPart(a) * primalPart(b) - primalPart(a) * dualPart(b)) / primalPart(b) ** 2;
-        return todual(f, df); 
-    }
-
-    operator **(a, b : real) where isDualType(a.type) {
-        var f = primalPart(a) ** b,
-            df = b * (primalPart(a) ** (b - 1)) * dualPart(a);
-        return todual(f, df);
-    }
-
-    proc sqrt(a) where isDualType(a.type) {
-        var f = sqrt(primalPart(a)),
-            df = 0.5 * dualPart(a) / sqrt(primalPart(a));
-        return todual(f, df);
-    }
-
-    proc cbrt(a) where isDualType(a.type) {
-        var f = cbrt(primalPart(a)),
-            df =  1.0 / 3.0 * dualPart(a) / cbrt(primalPart(a) ** 2);
-        return todual(f, df);
-    }
-}
\ No newline at end of file
diff --git a/src/hyperbolic.chpl b/src/hyperbolic.chpl
deleted file mode 100644
index 5e9f833..0000000
--- a/src/hyperbolic.chpl
+++ /dev/null
@@ -1,39 +0,0 @@
-module hyperbolic {
-    use ForwardModeAD;
-
-    proc sinh(a) where isDualType(a.type) {
-        var f = sinh(primalPart(a)), 
-        df = dualPart(a) * cosh(primalPart(a));
-        return todual(f, df);
-    }
-
-    proc cosh(a) where isDualType(a.type) {
-        var f = cosh(primalPart(a)), 
-            df = dualPart(a) * sinh(primalPart(a));
-        return todual(f, df);
-    }
-
-    proc tanh(a) where isDualType(a.type) {
-        var f = tanh(primalPart(a)), 
-            df = dualPart(a) * (1 - tanh(primalPart(a))**2);
-        return todual(f, df);
-    }
-
-    proc asinh(a) where isDualType(a.type) {
-        var f = asinh(primalPart(a)), 
-            df = dualPart(a) / sqrt(primalPart(a)**2 + 1);
-        return todual(f, df);
-    }
-
-    proc acosh(a) where isDualType(a.type) {
-        var f = acosh(primalPart(a)), 
-            df = dualPart(a) / sqrt(primalPart(a)**2 - 1);
-        return todual(f, df);
-    }
-
-    proc atanh(a) where isDualType(a.type) {
-        var f = atanh(primalPart(a)), 
-            df = dualPart(a) / (1 - primalPart(a) ** 2);
-        return todual(f, df);
-    }
-}
\ No newline at end of file
diff --git a/src/transcendental.chpl b/src/transcendental.chpl
deleted file mode 100644
index 36d2941..0000000
--- a/src/transcendental.chpl
+++ /dev/null
@@ -1,58 +0,0 @@
-module transcendental {
-    use ForwardModeAD;
-
-    operator **(a : real, b) where isDualType(b.type) {
-        var f = a ** primalPart(b),
-            df = log(a) * (a ** primalPart(b)) * dualPart(b);
-        return todual(f, df);
-    }
-
-
-    operator **(a, b) where isDualType(a.type) && a.type == b.type {
-        var f = primalPart(a) ** primalPart(b);
-        var df = f * (dualPart(b) * log(primalPart(a)) + primalPart(b) * dualPart(a) / primalPart(a)); 
-        return todual(f, df);
-    }
-
-    proc exp(a) where isDualType(a.type) {
-        var f = exp(primalPart(a)),
-            df = dualPart(a) * exp(primalPart(a));
-        return todual(f, df);
-    }
-
-    proc exp2(a) where isDualType(a.type) {
-        var f = exp2(primalPart(a)), 
-            df = ln2 * exp2(primalPart(a)) * dualPart(a);
-        return todual(f, df);
-    }
-
-    proc expm1(a) where isDualType(a.type) {
-        var f = expm1(primalPart(a)), 
-            df = exp(primalPart(a)) * dualPart(a);
-        return todual(f, df);
-    }
-
-    proc log(a) where isDualType(a.type) {
-        var f = log(primalPart(a)), 
-            df = dualPart(a) / primalPart(a);
-        return todual(f, df);
-    }
-    
-    proc log2(a) where isDualType(a.type) {
-        var f = log2(primalPart(a)), 
-            df = dualPart(a)/(primalPart(a) * ln2);
-        return todual(f, df);
-    }
-
-    proc log10(a) where isDualType(a.type) {
-        var f = log10(primalPart(a)), 
-            df = dualPart(a) / (primalPart(a) * ln10);
-        return todual(f, df);
-    }
-
-    proc log1p(a) where isDualType(a.type) {
-        var f = log1p(primalPart(a)), 
-            df = dualPart(a) / (primalPart(a) + 1);
-        return todual(f, df);
-    }
-}
diff --git a/src/trigonometric.chpl b/src/trigonometric.chpl
deleted file mode 100644
index 088b273..0000000
--- a/src/trigonometric.chpl
+++ /dev/null
@@ -1,39 +0,0 @@
-module trigonometric {
-    use ForwardModeAD;
-    
-    proc sin(a) where isDualType(a.type) {
-        var f = sin(primalPart(a)),
-            df = cos(primalPart(a)) * dualPart(a);
-        return todual(f, df);
-    }
-
-    proc cos(a) where isDualType(a.type) {
-        var f = cos(primalPart(a)),
-            df = -sin(primalPart(a)) * dualPart(a);
-        return todual(f, df);
-    }
-
-    proc tan(a) where isDualType(a.type) {
-        var f = tan(primalPart(a)),
-            df = dualPart(a) / (cos(primalPart(a)) ** 2);
-        return todual(f, df);
-    }
-
-    proc asin(a) where isDualType(a.type) {
-        var f = asin(primalPart(a)),
-            df = dualPart(a) / sqrt(1 - primalPart(a)**2);
-        return todual(f, df);
-    }
-
-    proc acos(a) where isDualType(a.type) {
-        var f = acos(primalPart(a)),
-            df = -dualPart(a) / sqrt(1 - primalPart(a)**2);
-        return todual(f, df);
-    }
-
-    proc atan(a) where isDualType(a.type) {
-        var f = atan(primalPart(a)),
-            df = dualPart(a) / (1 + primalPart(a)**2);
-        return todual(f, df);
-    }
-}
\ No newline at end of file
diff --git a/test/test_dual_algebra.chpl b/test/test_dual_algebra.chpl
index 9b5e397..98a92b2 100644
--- a/test/test_dual_algebra.chpl
+++ b/test/test_dual_algebra.chpl
@@ -1,5 +1,6 @@
 use UnitTest;
 use ForwardModeAD;
+use Math;
 
 proc test_arithmetic_operations(test : borrowed Test) throws {
     var x = new dual(1, 2),
diff --git a/test/test_multidual_algebra.chpl b/test/test_multidual_algebra.chpl
index 8fafba9..e5d9561 100644
--- a/test/test_multidual_algebra.chpl
+++ b/test/test_multidual_algebra.chpl
@@ -1,6 +1,6 @@
 use UnitTest;
 use ForwardModeAD;
-
+use Math;
 
 proc test_arithmetic_operations(test : borrowed Test) throws {
     var x = todual(1, [2, 3]),