diff --git a/README.md b/README.md
index 50f9be0e..643a111d 100644
--- a/README.md
+++ b/README.md
@@ -13,6 +13,6 @@
 
 ## About
 
-Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors.[^1] Mathematics.NET also supports first-order, forward and reverse-mode automatic differentiation.
+Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors.[^1] Mathematics.NET also supports first and second-order, forward and reverse-mode automatic differentiation.
 
 [^1]: Please visit the [documentation site](https://mathematics.hamlettanyavong.com) for detailed information.
diff --git a/docs/guide/autodiff/first-order-forward-mode.md b/docs/guide/autodiff/first-order-forward-mode.md
index dc64100c..4bb672b0 100644
--- a/docs/guide/autodiff/first-order-forward-mode.md
+++ b/docs/guide/autodiff/first-order-forward-mode.md
@@ -1,6 +1,6 @@
 # First-Order, Forward Mode Automatic Differentiation
 
-Support for first-order, forward-mode autodiff is provided by the `Dual<T>` class.
+Support for first-order, forward-mode autodiff is provided by the `Dual<T>` type.
 
 ## Dual Numbers
 
@@ -9,8 +9,13 @@ Forward-mode autodiff can be performed through the use of dual numbers which kee
 Dual<Real> x = new(1.23, 1.0);
 Dual<Real> y = new(2.34);
 ```
-The primal part represents the point at which we want to compute our derivative while the tangent part holds the information about our derivative. When we create a dual number with a tangent part, we specify a value that will be used as the seed, and it is important to know that the value of the derivative changes proportionally with this value.
+The primal part represents the point at which we want to compute our derivative while the tangent part holds the information about our derivative. When we create a dual number with a tangent part, we specify a value that will be used as the seed, and it is important to know that the value of the derivative changes proportionally with this value. We could also write, equivalently,
+```csharp
+using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
+var x = CreateVariable(1.23, 1.0);
+var y = CreateVariable(2.34);
+```
 Suppose we want to compute the partial derivative of the function
 $$
     f(x,y) = \frac{\sin(x + y)e^{-y}}{x^2+y^2+1}
@@ -19,11 +24,12 @@ at the points $ x=1.23 $ and $ y=2.34 $.
 with respect to $ x $. We must write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 // Add the following line to avoid having to type Dual<Real> every time we want to use a function.
 using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
-Dual<Real> x = new(1.23, 1.0);
-Dual<Real> y = new(2.34);
+Dual<Real> x = CreateVariable(1.23, 1.0);
+Dual<Real> y = CreateVariable(2.34);
 
 var result = Sin(x + y) * Exp(-y) / (x * x + y * y + 1);
 
@@ -31,8 +37,8 @@ Console.WriteLine("∂f/∂x: {0}", result);
 ```
 Notice that we set the seed for the variable of interest, $ x $, to 1 while the seed for the variable we do not care about, $ y $, was set to 0. If we had set both to 1.0, then we would have computed the total derivative of the function instead. To compute the partial derivative of our function with respect to $ y $, we write
 ```csharp
-Dual<Real> x = new(1.23);
-Dual<Real> y = new(2.34, 1.0);
+Dual<Real> x = CreateVariable(1.23);
+Dual<Real> y = CreateVariable(2.34, 1.0);
 ```
 with the tangent part of the variable of interest set to 1 and the other to 0. Doing this will print the following to the console:
 ```
@@ -40,11 +46,20 @@ with the tangent part of the variable of interest set to 1 and the other to 0. D
 ∂f/∂y: (-0.005009285670379789, -0.003024626925238263)
 ```
 
+### Total Derivative
+
+To get the total derivate of a function, set the seeds for each variable to `1.0`;
+```csharp
+Dual<Real> x = CreateVariable(1.23, 1.0);
+Dual<Real> y = CreateVariable(2.34, 1.0);
+// Repeat for each variable present
+```
+
 ## Dual Vectors
 
 We can create dual vectors to help us keep track of multiple dual numbers.
 ```csharp
-DualVector3<Real> x = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+DualVector3<Real> x = new(CreateVariable(1.23), CreateVariable(0.66), CreateVariable(2.34));
 ```
 We can use this to compute the vector-Jacobian product of the vector functions
 $$
@@ -57,9 +72,10 @@ $$
 with the vector $ v=(0.23,1.57,-1.71) $ at our points $ x_1=1.23 $, $ x_2=0.66 $, and $ x_3=2.34 $ by writing
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
-DualVector3<Real> x = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+DualVector3<Real> x = new(CreateVariable(1.23), CreateVariable(0.66), CreateVariable(2.34));
 Vector3<Real> v = new(0.23, 1.57, -1.71);
 
 var result = DualVector3<Real>.VJP(
diff --git a/docs/guide/autodiff/first-order-reverse-mode.md b/docs/guide/autodiff/first-order-reverse-mode.md
index 6415f538..02f58251 100644
--- a/docs/guide/autodiff/first-order-reverse-mode.md
+++ b/docs/guide/autodiff/first-order-reverse-mode.md
@@ -7,6 +7,7 @@ Support for first-order, reverse-mode automatic differentiation (autodiff) is pr
 Gradient tapes keep track of operations for autodiff; unlike forward-mode autodiff, tracking is required since gradients have to be calculated in reverse order. To begin using reverse-mode autodiff, we must create a gradient tape and assign it variables to track. These variables will be passed into and returned from methods that will compute the local gradients for us and record them on the tape.
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -33,6 +34,7 @@ $$
 at the point $ x=1.23 $. We can write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -107,6 +109,7 @@ Console.WriteLine(gradient[0]);
 The correct value for the derivative should be `3.525753368769319`. The complete code looks as follows:
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -289,6 +292,7 @@ $$
 at the points $ z=1.23+i2.34 $ and $ w=-0.66+i0.23 $. We can write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Complex> tape = new();
 var z = tape.CreateVariable(new(1.23, 2.34));
diff --git a/docs/guide/autodiff/second-order-forward-mode.md b/docs/guide/autodiff/second-order-forward-mode.md
new file mode 100644
index 00000000..518a8fbf
--- /dev/null
+++ b/docs/guide/autodiff/second-order-forward-mode.md
@@ -0,0 +1,40 @@
+# Second-Order, Forward Mode Automatic Differentiation
+
+Support for first-order, forward-mode autodiff is provided by the `HyperDual<T>` type. Because this type is used in a similar manner to `Dual<T>`, please refer to that section for help.
+
+### Second-Order Derivatives
+
+Suppose we wanted to find the second derivative of the complex function
+$$
+    f(z,w) = \sin(\tan{z}*\log{w})
+$$
+with respect to $ z $, at the points $ z=1.23+i0.66 $ and $ w=2.34-i0.25 $. We can do so by writing
+```csharp
+using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
+using static Mathematics.NET.AutoDiff.HyperDual<Mathematics.NET.Core.Complex>;
+
+var z = CreateVariable(new(1.23, 0.66), 1.0, 1.0);
+var w = CreateVariable(new(2.34, -0.25));
+
+var result = Sin(Tan(z) * Ln(w));
+
+Console.WriteLine(result.D3);
+```
+which will give us `(-6.158582087985498, 6.391603674636932)`. Notice that we now have to provide two seed values. Each one, very loosely speaking, "represents" a first-order derivative with respect to the variable in which it appears. Since there are two seeds present with the variable $ z $, it means we want to take its derivative twice. Similarly, we can write the following if we wanted the second derivative of our function with respect to $ w $:
+```csharp
+var z = CreateVariable(new(1.23, 0.66));
+var w = CreateVariable(new(2.34, -0.25), 1.0, 1.0);
+```
+This will give us `(0.30998196902728725, -0.11498565892578178)`. To get our mixed derivative, $ \partial f/\partial{z}\partial{w} $, we must indicate with we want one of each derivative
+```csharp
+var z = CreateVariable(new(1.23, 0.66), 1.0, 0.0);
+var w = CreateVariable(new(2.34, -0.25), 0.0, 1.0);
+```
+keeping in mind that the seeds must not occupy the same "slot." This, for example, will not give us the correct answer:
+```csharp
+// Incorrect, seeds must not occupy the same "slot"
+var z = CreateVariable(new(1.23, 0.66), 1.0, 0.0);
+var w = CreateVariable(new(2.34, -0.25), 1.0, 0.0);
+```
+This will print `(0.6670456012622978, 2.2955143408553718)` to the console.
diff --git a/docs/guide/autodiff/second-order-reverse-mode.md b/docs/guide/autodiff/second-order-reverse-mode.md
index 90107123..399082bb 100644
--- a/docs/guide/autodiff/second-order-reverse-mode.md
+++ b/docs/guide/autodiff/second-order-reverse-mode.md
@@ -1,6 +1,6 @@
 # Second-Order, Reverse Mode Automatic Differentiation
 
-Support for first-order, reverse-mode automatic differentiation (autodiff) is provided by the `HessianTape` class.
+Support for first-order, reverse-mode autodiff is provided by the `HessianTape` class.
 
 ## Hessian Tapes
 
@@ -26,7 +26,7 @@ $$
 $$
 Note that, in the future, we will not have to do this manually since there will be a method made specifically to compute Laplacians in spherical coordinates. For now, if we wanted to compute the Laplacian of the function
 $$
-    f(x,y,z) = \frac{\cos(x)}{(x+y)\sin(z)}
+    f(r,\theta,\phi) = \frac{\cos(r)}{(r+\theta)\sin(\phi)}
 $$
 we can write
 ```csharp
@@ -36,7 +36,7 @@ using Mathematics.NET.Core;
 HessianTape<Real> tape = new();
 var x = tape.CreateVariableVector(1.23, 0.66, 0.23);
 
-// f(x, y, z) = cos(x) / ((x + y) * sin(z))
+// f(r, θ, ϕ) = cos(r) / ((r + θ) * sin(ϕ))
 _ = tape.Divide(
         tape.Cos(x.X1),
         tape.Multiply(
diff --git a/docs/guide/toc.yml b/docs/guide/toc.yml
index 1feffb99..8eb58c91 100644
--- a/docs/guide/toc.yml
+++ b/docs/guide/toc.yml
@@ -16,3 +16,5 @@
     href: autodiff/first-order-forward-mode.md
   - name: Second-Order, Reverse-Mode Automatic Differentiation
     href: autodiff/second-order-reverse-mode.md
+  - name: Second-Order, Forward-Mode Automatic Differentiation
+    href: autodiff/second-order-forward-mode.md
diff --git a/docs/index.md b/docs/index.md
index 71aa48fa..15147c00 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -54,4 +54,4 @@
 
 ## About
 
-Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors. Mathematics.NET also supports forward and reverse-mode automatic differentiation.
+Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors. Mathematics.NET also supports first and second-order, forward and reverse-mode automatic differentiation.
diff --git a/src/Mathematics.NET/AutoDiff/AutoDiffExtensions.cs b/src/Mathematics.NET/AutoDiff/AutoDiffExtensions.cs
index 570a0624..7b3e0503 100644
--- a/src/Mathematics.NET/AutoDiff/AutoDiffExtensions.cs
+++ b/src/Mathematics.NET/AutoDiff/AutoDiffExtensions.cs
@@ -33,19 +33,19 @@ namespace Mathematics.NET.AutoDiff;
 public static class AutoDiffExtensions
 {
     //
-    // Variable creation
+    // Variable vector creation
     //
 
-    /// <summary>Create a three-element vector from a seed vector of length three.</summary>
+    /// <summary>Create a dual vector from a seed vector.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A type that implements <see cref="ITape{T}"/></param>
-    /// <param name="x">A three-element vector of seed values</param>
+    /// <param name="x">A vector of seed values</param>
     /// <returns>A variable vector of length three</returns>
     public static VariableVector3<T> CreateVariableVector<T>(this ITape<T> tape, Vector3<T> x)
         where T : IComplex<T>, IDifferentiableFunctions<T>
         => new(tape.CreateVariable(x.X1), tape.CreateVariable(x.X2), tape.CreateVariable(x.X3));
 
-    /// <summary>Create a three-element vector from seed values.</summary>
+    /// <summary>Create a vector from seed values.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A type that implements <see cref="ITape{T}"/></param>
     /// <param name="x1Seed">The first seed value</param>
@@ -159,6 +159,25 @@ public static Vector3<T> Gradient<T>(this ITape<T> tape, Func<ITape<T>, Variable
         return new(gradients[0], gradients[1], gradients[2]);
     }
 
+    /// <summary>Compute the Hessian of a scaler function using reverse-mode automatic differentiation.</summary>
+    /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
+    /// <param name="tape">A gradient or Hessian tape</param>
+    /// <param name="f">A scalar function</param>
+    /// <param name="x">The point at which to compute the gradient</param>
+    /// <returns>The Hessian of the scalar function</returns>
+    public static Matrix3x3<T> Hessian<T>(this HessianTape<T> tape, Func<HessianTape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
+        where T : IComplex<T>, IDifferentiableFunctions<T>
+    {
+        _ = f(tape, x);
+
+        tape.ReverseAccumulation(out ReadOnlySpan2D<T> hessian);
+
+        return new(
+            hessian[0, 0], hessian[0, 1], hessian[0, 2],
+            hessian[1, 0], hessian[1, 1], hessian[1, 2],
+            hessian[2, 0], hessian[2, 1], hessian[2, 2]);
+    }
+
     /// <summary>Compute the Jacobian of a vector function using reverse-mode automatic differentiation: $ \nabla^\text{T}f_i(\textbf{x}) $ for $ i=\left\{1,2,3\right\} $.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A gradient or Hessian tape</param>
@@ -233,7 +252,7 @@ public static Vector3<T> JVP<T>(
     /// <param name="f">A scalar function</param>
     /// <param name="x">The point at which to compute the Laplacian</param>
     /// <returns>The Laplacian of the scalar function</returns>
-    public static T Laplacian<T>(this HessianTape<T> tape, Func<ITape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
+    public static T Laplacian<T>(this HessianTape<T> tape, Func<HessianTape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
         where T : IComplex<T>, IDifferentiableFunctions<T>
     {
         _ = f(tape, x);
diff --git a/src/Mathematics.NET/AutoDiff/Dual.cs b/src/Mathematics.NET/AutoDiff/Dual.cs
index d8f21292..abec18d9 100644
--- a/src/Mathematics.NET/AutoDiff/Dual.cs
+++ b/src/Mathematics.NET/AutoDiff/Dual.cs
@@ -27,8 +27,6 @@
 
 using System.Diagnostics.CodeAnalysis;
 using System.Runtime.InteropServices;
-using Mathematics.NET.Core.Operations;
-using Mathematics.NET.Core.Relations;
 
 namespace Mathematics.NET.AutoDiff;
 
@@ -38,15 +36,7 @@ namespace Mathematics.NET.AutoDiff;
 /// <param name="d1">The tangent part of the dual number</param>
 [Serializable]
 [StructLayout(LayoutKind.Sequential)]
-public readonly struct Dual<T>(T d0, T d1)
-    : IAdditionOperation<Dual<T>, Dual<T>>,
-      IDivisionOperation<Dual<T>, Dual<T>>,
-      IMultiplicationOperation<Dual<T>, Dual<T>>,
-      INegationOperation<Dual<T>, Dual<T>>,
-      ISubtractionOperation<Dual<T>, Dual<T>>,
-      IEqualityRelation<Dual<T>, bool>,
-      IEquatable<Dual<T>>,
-      IFormattable
+public readonly struct Dual<T>(T d0, T d1) : IDual<Dual<T>, T>
     where T : IComplex<T>, IDifferentiableFunctions<T>
 {
     private readonly T _d0 = d0;
@@ -54,7 +44,6 @@ public readonly struct Dual<T>(T d0, T d1)
 
     public Dual(T value) : this(value, T.Zero) { }
 
-    /// <summary>Represents the primal part of the dual number</summary>
     public readonly T D0 => _d0;
 
     /// <summary>Represents the tangent part of the dual number</summary>
@@ -89,7 +78,7 @@ public Dual(T value) : this(value, T.Zero) { }
 
     public static Dual<T> operator /(T c, Dual<T> x) => new(c / x._d0, -x._d1 * c / (x._d0 * x._d0));
 
-    public static Dual<T> operator /(Dual<T> x, T c) => new(x._d0 / c, x._d1 * c / (c * c));
+    public static Dual<T> operator /(Dual<T> x, T c) => new(x._d0 / c, x._d1 / c);
 
     public static Dual<Real> Modulo(Dual<Real> x, Dual<Real> y)
     {
@@ -124,23 +113,26 @@ public static Dual<Real> Modulo(Dual<Real> x, Dual<Real> y)
     // Other operations
     //
 
+    public static Dual<T> CreateVariable(T value) => new(value, T.Zero);
+
+    public static Dual<T> CreateVariable(T value, T seed) => new(value, seed);
+
+    public Dual<T> WithSeed(T seed) => new(_d0, seed);
+
     // Exponential functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp(T)"/>
     public static Dual<T> Exp(Dual<T> x)
     {
         var exp = T.Exp(x._d0);
         return new(exp, x._d1 * exp);
     }
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp2(T)"/>
     public static Dual<T> Exp2(Dual<T> x)
     {
         var exp2 = T.Exp2(x._d0);
         return new(exp2, Real.Ln2 * x._d1 * exp2);
     }
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp10(T)"/>
     public static Dual<T> Exp10(Dual<T> x)
     {
         var exp10 = T.Exp10(x._d0);
@@ -149,25 +141,19 @@ public static Dual<T> Exp10(Dual<T> x)
 
     // Hyperbolic functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Acosh(T)"/>
     public static Dual<T> Acosh(Dual<T> x)
         => new(T.Acosh(x._d0), x._d1 / (T.Sqrt(x._d0 - T.One) * T.Sqrt(x._d0 + T.One)));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Asinh(T)"/>
     public static Dual<T> Asinh(Dual<T> x)
         => new(T.Asinh(x._d0), x._d1 / T.Sqrt(x._d0 * x._d0 + T.One));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Atanh(T)"/>
     public static Dual<T> Atanh(Dual<T> x)
         => new(T.Atanh(x._d0), x._d1 / (T.One - x._d0 * x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cosh(T)"/>
     public static Dual<T> Cosh(Dual<T> x) => new(T.Cosh(x._d0), x._d1 * T.Sinh(x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sinh(T)"/>
     public static Dual<T> Sinh(Dual<T> x) => new(T.Sinh(x._d0), x._d1 * T.Cosh(x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Tanh(T)"/>
     public static Dual<T> Tanh(Dual<T> x)
     {
         var u = T.One / T.Cosh(x._d0);
@@ -176,36 +162,28 @@ public static Dual<T> Tanh(Dual<T> x)
 
     // Logarithmic functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Ln(T)"/>
     public static Dual<T> Ln(Dual<T> x) => new(T.Ln(x._d0), x._d1 / x._d0);
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log(T, T)"/>
     public static Dual<T> Log(Dual<T> x, Dual<T> b) => Ln(x) / Ln(b);
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log2(T)"/>
     public static Dual<T> Log2(Dual<T> x) => new(T.Log2(x._d0), x._d1 / (Real.Ln2 * x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log10(T)"/>
     public static Dual<T> Log10(Dual<T> x) => new(T.Log10(x._d0), x._d1 / (Real.Ln10 * x._d0));
 
     // Power functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Pow(T, T)"/>
     public static Dual<T> Pow(Dual<T> x, Dual<T> y) => Exp(y * Ln(x));
 
     // Root functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cbrt(T)"/>
     public static Dual<T> Cbrt(Dual<T> x)
     {
         var cbrt = T.Cbrt(x._d0);
         return new(cbrt, x._d1 / (3.0 * cbrt * cbrt));
     }
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Root(T, T)"/>
     public static Dual<T> Root(Dual<T> x, Dual<T> n) => Exp(Ln(x) / n);
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sqrt(T)"/>
     public static Dual<T> Sqrt(Dual<T> x)
     {
         var sqrt = T.Sqrt(x._d0);
@@ -214,13 +192,10 @@ public static Dual<T> Sqrt(Dual<T> x)
 
     // Trigonometric functions
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Acos(T)"/>
     public static Dual<T> Acos(Dual<T> x) => new(T.Acos(x._d0), -x._d1 / T.Sqrt(T.One - x._d0 * x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Asin(T)"/>
     public static Dual<T> Asin(Dual<T> x) => new(T.Asin(x._d0), x._d1 / T.Sqrt(T.One - x._d0 * x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Atan(T)"/>
     public static Dual<T> Atan(Dual<T> x) => new(T.Atan(x._d0), x._d1 / (T.One + x._d0 * x._d0));
 
     /// <inheritdoc cref="IReal{T}.Atan2(T, T)"/>
@@ -230,13 +205,10 @@ public static Dual<Real> Atan2(Dual<Real> y, Dual<Real> x)
         return new(Real.Atan2(y._d0, x._d0), (y._d1 * x._d0 - x._d1 * y._d0) * u);
     }
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cos(T)"/>
     public static Dual<T> Cos(Dual<T> x) => new(T.Cos(x._d0), -x._d1 * T.Sin(x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sin(T)"/>
     public static Dual<T> Sin(Dual<T> x) => new(T.Sin(x._d0), x._d1 * T.Cos(x._d0));
 
-    /// <inheritdoc cref="IDifferentiableFunctions{T}.Tan(T)"/>
     public static Dual<T> Tan(Dual<T> x)
     {
         var sec = T.One / T.Cos(x._d0);
@@ -265,11 +237,9 @@ public static Dual<T> CustomOperation(Dual<T> x, Dual<T> y, Func<T, T, T> f, Fun
         => new(f(x._d0, y._d0), dfy(x._d0, y._d0) * x._d1 + dfx(x._d0, y._d1) * y._d1);
 
     //
-    // Explicit operators
+    // Dual vector creation
     //
 
-    /// <summary>Convert a value of type <see cref="double"/> to one of type <see cref="Dual{T}"/></summary>
-    /// <remarks>The tangent part of the converted value will be zero.</remarks>
-    /// <param name="x">The value to convert</param>
-    public static explicit operator Dual<T>(double x) => new(x);
+    public static DualVector3<Dual<T>, T> CreateDualVector(Dual<T> x1Seed, Dual<T> x2Seed, Dual<T> x3Seed)
+        => new(x1Seed, x2Seed, x3Seed);
 }
diff --git a/src/Mathematics.NET/AutoDiff/DualVector3.cs b/src/Mathematics.NET/AutoDiff/DualVector3.cs
index c36b04f7..39c92a79 100644
--- a/src/Mathematics.NET/AutoDiff/DualVector3.cs
+++ b/src/Mathematics.NET/AutoDiff/DualVector3.cs
@@ -34,20 +34,22 @@ namespace Mathematics.NET.AutoDiff;
 
 /// <summary>Represents a vector of three dual numbers for use in forward-mode automatic differentiation</summary>
 /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/></typeparam>
+/// <typeparam name="U">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
 [StructLayout(LayoutKind.Sequential)]
-public record struct DualVector3<T>
-    where T : IComplex<T>, IDifferentiableFunctions<T>
+public record struct DualVector3<T, U>
+    where T : IDual<T, U>
+    where U : IComplex<U>, IDifferentiableFunctions<U>
 {
     /// <summary>The first element of the vector</summary>
-    public Dual<T> X1;
+    public T X1;
 
     /// <summary>The second element of the vector</summary>
-    public Dual<T> X2;
+    public T X2;
 
     /// <summary>The third element of the vector</summary>
-    public Dual<T> X3;
+    public T X3;
 
-    public DualVector3(Dual<T> x1, Dual<T> x2, Dual<T> x3)
+    public DualVector3(T x1, T x2, T x3)
     {
         X1 = x1;
         X2 = x2;
@@ -66,7 +68,7 @@ public T this[int index]
 
     // Get
 
-    internal static T GetElement(DualVector3<T> vector, int index)
+    internal static T GetElement(DualVector3<T, U> vector, int index)
     {
         if ((uint)index >= 3)
         {
@@ -77,31 +79,31 @@ internal static T GetElement(DualVector3<T> vector, int index)
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static T GetElementUnsafe(ref DualVector3<T> vector, int index)
+    private static T GetElementUnsafe(ref DualVector3<T, U> vector, int index)
     {
         Debug.Assert(index is >= 0 and < 3);
-        return Unsafe.Add(ref Unsafe.As<DualVector3<T>, T>(ref vector), index);
+        return Unsafe.Add(ref Unsafe.As<DualVector3<T, U>, T>(ref vector), index);
     }
 
     // Set
 
-    internal static DualVector3<T> WithElement(DualVector3<T> vector, int index, T value)
+    internal static DualVector3<T, U> WithElement(DualVector3<T, U> vector, int index, T value)
     {
         if ((uint)index >= 3)
         {
             throw new IndexOutOfRangeException();
         }
 
-        DualVector3<T> result = vector;
+        DualVector3<T, U> result = vector;
         SetElementUnsafe(ref result, index, value);
         return result;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static void SetElementUnsafe(ref DualVector3<T> vector, int index, T value)
+    private static void SetElementUnsafe(ref DualVector3<T, U> vector, int index, T value)
     {
         Debug.Assert(index is >= 0 and < 3);
-        Unsafe.Add(ref Unsafe.As<DualVector3<T>, T>(ref vector), index) = value;
+        Unsafe.Add(ref Unsafe.As<DualVector3<T, U>, T>(ref vector), index) = value;
     }
 
     //
@@ -109,7 +111,7 @@ private static void SetElementUnsafe(ref DualVector3<T> vector, int index, T val
     //
 
     /// <inheritdoc cref="Vector3{T}.Cross(Vector3{T}, Vector3{T})"/>
-    public static DualVector3<T> Cross(DualVector3<T> left, DualVector3<T> right)
+    public static DualVector3<T, U> Cross(DualVector3<T, U> left, DualVector3<T, U> right)
     {
         return new(
             left.X2 * right.X3 - left.X3 * right.X2,
@@ -127,20 +129,20 @@ public static DualVector3<T> Cross(DualVector3<T> left, DualVector3<T> right)
     /// <param name="fz">The z-component of the vector field</param>
     /// <param name="x">The point at which to compute the curl</param>
     /// <returns>The curl of the vector field</returns>
-    public static Vector3<T> Curl(
-        Func<DualVector3<T>, Dual<T>> fx,
-        Func<DualVector3<T>, Dual<T>> fy,
-        Func<DualVector3<T>, Dual<T>> fz,
-        DualVector3<T> x)
+    public static Vector3<U> Curl(
+        Func<DualVector3<T, U>, T> fx,
+        Func<DualVector3<T, U>, T> fy,
+        Func<DualVector3<T, U>, T> fz,
+        DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, Real.One), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.One), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.One))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
 
-        ReadOnlySpan<T> dfx = [fx(seeds[0]).D1, fx(seeds[1]).D1, fx(seeds[2]).D1];
-        ReadOnlySpan<T> dfy = [fy(seeds[0]).D1, fy(seeds[1]).D1, fy(seeds[2]).D1];
-        ReadOnlySpan<T> dfz = [fz(seeds[0]).D1, fz(seeds[1]).D1, fz(seeds[2]).D1];
+        ReadOnlySpan<U> dfx = [fx(seeds[0]).D1, fx(seeds[1]).D1, fx(seeds[2]).D1];
+        ReadOnlySpan<U> dfy = [fy(seeds[0]).D1, fy(seeds[1]).D1, fy(seeds[2]).D1];
+        ReadOnlySpan<U> dfz = [fz(seeds[0]).D1, fz(seeds[1]).D1, fz(seeds[2]).D1];
 
         return new(
             dfz[1] - dfy[2],
@@ -153,15 +155,15 @@ public static Vector3<T> Curl(
     /// <param name="f">A scalar function</param>
     /// <param name="x">The point at which to compute the directional derivative</param>
     /// <returns>A directional derivative</returns>
-    public static T DirectionalDerivative(
-        Vector3<T> v,
-        Func<DualVector3<T>, Dual<T>> f,
-        DualVector3<T> x)
+    public static U DirectionalDerivative(
+        Vector3<U> v,
+        Func<DualVector3<T, U>, T> f,
+        DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, v.X1), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, v.X2), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, v.X3))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(v.X1), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(v.X2), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(v.X3))];
 
         return f(seeds[0]).D1 + f(seeds[1]).D1 + f(seeds[2]).D1;
     }
@@ -172,16 +174,16 @@ public static T DirectionalDerivative(
     /// <param name="fz">The z-component of the vector field</param>
     /// <param name="x">The point at which to compute the divergence</param>
     /// <returns>The divergence of the vector field</returns>
-    public static T Divergence(
-        Func<DualVector3<T>, Dual<T>> fx,
-        Func<DualVector3<T>, Dual<T>> fy,
-        Func<DualVector3<T>, Dual<T>> fz,
-        DualVector3<T> x)
+    public static U Divergence(
+        Func<DualVector3<T, U>, T> fx,
+        Func<DualVector3<T, U>, T> fy,
+        Func<DualVector3<T, U>, T> fz,
+        DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, Real.One), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.One), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.One))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
 
         return fx(seeds[0]).D1 + fy(seeds[1]).D1 + fz(seeds[2]).D1;
     }
@@ -190,34 +192,63 @@ public static T Divergence(
     /// <param name="f">A scalar function</param>
     /// <param name="x">The point at which to compute the gradient</param>
     /// <returns>The gradient of the scalar function</returns>
-    public static Vector3<T> Gradient(Func<DualVector3<T>, Dual<T>> f, DualVector3<T> x)
+    public static Vector3<U> Gradient(Func<DualVector3<T, U>, T> f, DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, Real.One), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.One), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.One))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
 
         return new(f(seeds[0]).D1, f(seeds[1]).D1, f(seeds[2]).D1);
     }
 
+    /// <summary>Compute the Hessian of a scalar function using forward-mode automatic differentiation.</summary>
+    /// <param name="f">A scalar function</param>
+    /// <param name="x">The point at which to compute the gradient</param>
+    /// <returns>The Hessian of the scalar function</returns>
+    public static Matrix3x3<U> Hessian(Func<DualVector3<HyperDual<U>, U>, HyperDual<U>> f, DualVector3<HyperDual<U>, U> x)
+    {
+        ReadOnlySpan<DualVector3<HyperDual<U>, U>> seeds = [
+            new(x.X1.WithSeed(U.One, U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
+
+        Matrix3x3<U> result = new();
+
+        result[0, 0] = f(new(x.X1.WithSeed(U.One, U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero))).D3;
+        result[1, 1] = f(new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One, U.One), x.X3.WithSeed(U.Zero))).D3;
+        result[2, 2] = f(new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One, U.One))).D3;
+
+        var e12 = f(new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero, U.One), x.X3.WithSeed(U.Zero)));
+        result[0, 1] = e12.D3; result[1, 0] = e12.D3;
+
+        var e13 = f(new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero, U.One)));
+        result[0, 2] = e13.D3; result[2, 0] = e13.D3;
+
+        var e23 = f(new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero, U.One)));
+        result[1, 2] = e23.D3; result[2, 1] = e23.D3;
+
+        return result;
+    }
+
     /// <summary>Compute the Jacobian of a vector function using forward-mode automatic differentiation: $ \nabla^\text{T}f_i(\textbf{x}) $ for $ i=\left\{1,2,3\right\} $.</summary>
     /// <param name="fx">The first function</param>
     /// <param name="fy">The second function</param>
     /// <param name="fz">The third function</param>
     /// <param name="x">The point at which to compute the Jacobian</param>
     /// <returns>The Jacobian of the vector function</returns>
-    public static Matrix3x3<T> Jacobian(
-        Func<DualVector3<T>, Dual<T>> fx,
-        Func<DualVector3<T>, Dual<T>> fy,
-        Func<DualVector3<T>, Dual<T>> fz,
-        DualVector3<T> x)
+    public static Matrix3x3<U> Jacobian(
+        Func<DualVector3<T, U>, T> fx,
+        Func<DualVector3<T, U>, T> fy,
+        Func<DualVector3<T, U>, T> fz,
+        DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, Real.One), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.One), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.One))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
 
-        Matrix3x3<T> result = new();
+        Matrix3x3<U> result = new();
 
         result[0, 0] = fx(seeds[0]).D1; result[0, 1] = fx(seeds[1]).D1; result[0, 2] = fx(seeds[2]).D1;
         result[1, 0] = fy(seeds[0]).D1; result[1, 1] = fy(seeds[1]).D1; result[1, 2] = fy(seeds[2]).D1;
@@ -233,17 +264,31 @@ public static Matrix3x3<T> Jacobian(
     /// <param name="x">The point at which to compute the Jacobian-vector product</param>
     /// <param name="v">A vector</param>
     /// <returns>The Jacobian-vector product of the vector function and vector</returns>
-    public static Vector3<T> JVP(
-        Func<DualVector3<T>, Dual<T>> fx,
-        Func<DualVector3<T>, Dual<T>> fy,
-        Func<DualVector3<T>, Dual<T>> fz,
-        DualVector3<T> x,
-        Vector3<T> v)
+    public static Vector3<U> JVP(
+        Func<DualVector3<T, U>, T> fx,
+        Func<DualVector3<T, U>, T> fy,
+        Func<DualVector3<T, U>, T> fz,
+        DualVector3<T, U> x,
+        Vector3<U> v)
     {
-        DualVector3<T> seed = new(new(x.X1.D0, v.X1), new(x.X2.D0, v.X2), new(x.X3.D0, v.X3));
+        DualVector3<T, U> seed = new(x.X1.WithSeed(v.X1), x.X2.WithSeed(v.X2), x.X3.WithSeed(v.X3));
         return new(fx(seed).D1, fy(seed).D1, fz(seed).D1);
     }
 
+    /// <summary>Compute the Laplacian of a scalar function using forward-mode automatic differentiation:  $ \nabla^2f $.</summary>
+    /// <param name="f">A scalar function</param>
+    /// <param name="x">The point at which to compute the gradient</param>
+    /// <returns>The gradient of the scalar function</returns>
+    public static U Laplacian(Func<DualVector3<HyperDual<U>, U>, HyperDual<U>> f, DualVector3<HyperDual<U>, U> x)
+    {
+        ReadOnlySpan<DualVector3<HyperDual<U>, U>> seeds = [
+            new(x.X1.WithSeed(U.One, U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One, U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One, U.One))];
+
+        return f(seeds[0]).D3 + f(seeds[1]).D3 + f(seeds[2]).D3;
+    }
+
     /// <summary>Compute the vector-Jacobian product of a vector and a vector function using forward-mode automatic differentiation.</summary>
     /// <param name="v">A vector</param>
     /// <param name="fx">The first function</param>
@@ -251,19 +296,19 @@ public static Vector3<T> JVP(
     /// <param name="fz">The third function</param>
     /// <param name="x">The point at which to compute the vector-Jacobian product</param>
     /// <returns>The vector-Jacobian product of the vector and vector-function</returns>
-    public static Vector3<T> VJP(
-        Vector3<T> v,
-        Func<DualVector3<T>, Dual<T>> fx,
-        Func<DualVector3<T>, Dual<T>> fy,
-        Func<DualVector3<T>, Dual<T>> fz,
-        DualVector3<T> x)
+    public static Vector3<U> VJP(
+        Vector3<U> v,
+        Func<DualVector3<T, U>, T> fx,
+        Func<DualVector3<T, U>, T> fy,
+        Func<DualVector3<T, U>, T> fz,
+        DualVector3<T, U> x)
     {
-        ReadOnlySpan<DualVector3<T>> seeds = [
-            new(new(x.X1.D0, Real.One), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.One), new(x.X3.D0, Real.Zero)),
-            new(new(x.X1.D0, Real.Zero), new(x.X2.D0, Real.Zero), new(x.X3.D0, Real.One))];
+        ReadOnlySpan<DualVector3<T, U>> seeds = [
+            new(x.X1.WithSeed(U.One), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.One), x.X3.WithSeed(U.Zero)),
+            new(x.X1.WithSeed(U.Zero), x.X2.WithSeed(U.Zero), x.X3.WithSeed(U.One))];
 
-        Vector3<T> result = new();
+        Vector3<U> result = new();
 
         result[0] = v.X1 * fx(seeds[0]).D1 + v.X2 * fy(seeds[0]).D1 + v.X3 * fz(seeds[0]).D1;
         result[1] = v.X1 * fx(seeds[1]).D1 + v.X2 * fy(seeds[1]).D1 + v.X3 * fz(seeds[1]).D1;
diff --git a/src/Mathematics.NET/AutoDiff/HessianTape.cs b/src/Mathematics.NET/AutoDiff/HessianTape.cs
index 47846a51..6d665124 100644
--- a/src/Mathematics.NET/AutoDiff/HessianTape.cs
+++ b/src/Mathematics.NET/AutoDiff/HessianTape.cs
@@ -173,8 +173,8 @@ public void ReverseAccumulation(out ReadOnlySpan2D<T> hessian, T seed)
             var node = Unsafe.Add(ref start, i);
             var gradientElement = gradientSpan[i];
 
-            EdgePush(hessianSpan, ref node, i);
-            Accumulate(hessianSpan, ref node, gradientElement);
+            EdgePush(hessianSpan, in node, i);
+            Accumulate(hessianSpan, in node, gradientElement);
 
             gradientSpan[node.PX] += gradientElement * node.DX;
             gradientSpan[node.PY] += gradientElement * node.DY;
@@ -214,8 +214,8 @@ public void ReverseAccumulation(out ReadOnlySpan<T> gradient, out ReadOnlySpan2D
             var node = Unsafe.Add(ref start, i);
             var gradientElement = gradientSpan[i];
 
-            EdgePush(hessianSpan, ref node, i);
-            Accumulate(hessianSpan, ref node, gradientElement);
+            EdgePush(hessianSpan, in node, i);
+            Accumulate(hessianSpan, in node, gradientElement);
 
             gradientSpan[node.PX] += gradientElement * node.DX;
             gradientSpan[node.PY] += gradientElement * node.DY;
@@ -226,7 +226,7 @@ public void ReverseAccumulation(out ReadOnlySpan<T> gradient, out ReadOnlySpan2D
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
-    private static void EdgePush(Span2D<T> weight, ref HessianNode<T> node, int i)
+    private static void EdgePush(Span2D<T> weight, in HessianNode<T> node, int i)
     {
         for (int p = 0; p <= i; p++)
         {
@@ -270,7 +270,7 @@ private static void EdgePush(Span2D<T> weight, ref HessianNode<T> node, int i)
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
-    private static void Accumulate(Span2D<T> weight, ref HessianNode<T> node, T v)
+    private static void Accumulate(Span2D<T> weight, in HessianNode<T> node, T v)
     {
         weight[node.PX, node.PX] += v * node.DXX;
         weight[node.PX, node.PY] += v * node.DXY;
diff --git a/src/Mathematics.NET/AutoDiff/HyperDual.cs b/src/Mathematics.NET/AutoDiff/HyperDual.cs
new file mode 100644
index 00000000..5e6cd7cc
--- /dev/null
+++ b/src/Mathematics.NET/AutoDiff/HyperDual.cs
@@ -0,0 +1,259 @@
+// <copyright file="HyperDual.cs" company="Mathematics.NET">
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+// </copyright>
+
+using System.Diagnostics.CodeAnalysis;
+using System.Runtime.InteropServices;
+
+namespace Mathematics.NET.AutoDiff;
+
+/// <summary>Represents a hyper-dual number</summary>
+/// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
+/// <param name="d0">The primal part of a hyper-dual number</param>
+/// <param name="d1">The tangent part of a hyper-dual number</param>
+[Serializable]
+[StructLayout(LayoutKind.Sequential)]
+public readonly struct HyperDual<T>(Dual<T> d0, Dual<T> d1) : IDual<HyperDual<T>, T>
+    where T : IComplex<T>, IDifferentiableFunctions<T>
+{
+    private readonly Dual<T> _d0 = d0;
+    private readonly Dual<T> _d1 = d1;
+
+    public HyperDual(Dual<T> value) : this(value, new(T.Zero)) { }
+
+    /// <summary>Represents the primal part of the hyper-dual number</summary>
+    public T D0 => _d0.D0;
+
+    /// <summary>Represents the first tangent part of the hyper-dual number</summary>
+    public T D1 => _d0.D1;
+
+    /// <summary>Represents the second tangent part of the hyper-dual number</summary>
+    public T D2 => _d1.D0;
+
+    /// <summary>Represents the third tangent part, the cross term, of a hyper-dual number</summary>
+    public T D3 => _d1.D1;
+
+    //
+    // Operators
+    //
+
+    public static HyperDual<T> operator -(HyperDual<T> x) => new(-x._d0, -x._d1);
+
+    public static HyperDual<T> operator +(HyperDual<T> x, HyperDual<T> y) => new(x._d0 + y._d0, x._d1 + y._d1);
+
+    public static HyperDual<T> operator +(T c, HyperDual<T> x) => new(c + x._d0, x._d1);
+
+    public static HyperDual<T> operator +(HyperDual<T> x, T c) => new(x._d0 + c, x._d1);
+
+    public static HyperDual<T> operator -(HyperDual<T> x, HyperDual<T> y) => new(x._d0 - y._d0, x._d1 - y._d1);
+
+    public static HyperDual<T> operator -(T c, HyperDual<T> x) => new(c - x._d0, -x._d1);
+
+    public static HyperDual<T> operator -(HyperDual<T> x, T c) => new(x._d0 - c, x._d1);
+
+    public static HyperDual<T> operator *(HyperDual<T> x, HyperDual<T> y) => new(x._d0 * y._d0, x._d0 * y._d1 + y._d0 * x._d1);
+
+    public static HyperDual<T> operator *(T c, HyperDual<T> x) => new(c * x._d0, c * x._d1);
+
+    public static HyperDual<T> operator *(HyperDual<T> x, T c) => new(x._d0 * c, c * x._d1);
+
+    public static HyperDual<T> operator /(HyperDual<T> x, HyperDual<T> y)
+        => new(x._d0 / y._d0, (x._d1 * y._d0 - y._d1 * x._d0) / (y._d0 * y._d0));
+
+    public static HyperDual<T> operator /(T c, HyperDual<T> x) => new(c / x._d0, -x._d1 * c / (x._d0 * x._d0));
+
+    public static HyperDual<T> operator /(HyperDual<T> x, T c) => new(x._d0 / c, x._d1 / c);
+
+    //
+    // Equality
+    //
+
+    public static bool operator ==(HyperDual<T> left, HyperDual<T> right) => left._d0 == right._d0 && left._d1 == right._d1;
+
+    public static bool operator !=(HyperDual<T> left, HyperDual<T> right) => left._d0 != right._d0 && left._d1 != right._d1;
+
+    public override bool Equals([NotNullWhen(true)] object? obj) => obj is HyperDual<T> other && Equals(other);
+
+    public bool Equals(HyperDual<T> value) => _d0.Equals(value._d0) && _d1.Equals(value._d1);
+
+    public override int GetHashCode() => HashCode.Combine(_d0, _d0);
+
+    //
+    // Formatting
+    //
+
+    public string ToString(string? format, IFormatProvider? formatProvider) => $"({_d0.D0}, {_d0.D1}, {_d1.D0} {_d1.D1})";
+
+    //
+    // Other operations
+    //
+
+    public static HyperDual<T> CreateVariable(T value) => new(new(value, T.Zero));
+
+    public static HyperDual<T> CreateVariable(T value, T seed) => new(new(value, T.One));
+
+    /// <summary>Create an instance of the type with a specified value and two seed values.</summary>
+    /// <remarks>Use this to compute mixed second derivatives.</remarks>
+    /// <param name="value">A value</param>
+    /// <param name="e1Seed">A seed value for the first variable of interest</param>
+    /// <param name="e2Seed">A seec value for the second variable of interest</param>
+    /// <returns>An instance of the type</returns>
+    public static HyperDual<T> CreateVariable(T value, T e1Seed, T e2Seed) => new(new(value, e1Seed), new(e2Seed));
+
+    public HyperDual<T> WithSeed(T seed) => new(new(_d0.D0, seed));
+
+    public HyperDual<T> WithSeed(T e1Seed, T e2Seed) => new(new(_d0.D0, e1Seed), new(e2Seed));
+
+    // Exponential functions
+
+    public static HyperDual<T> Exp(HyperDual<T> x)
+    {
+        var exp = Dual<T>.Exp(x._d0);
+        return new(exp, x._d1 * exp);
+    }
+
+    public static HyperDual<T> Exp2(HyperDual<T> x)
+    {
+        var exp2 = Dual<T>.Exp2(x._d0);
+        return new(exp2, Real.Ln2 * x._d1 * exp2);
+    }
+
+    public static HyperDual<T> Exp10(HyperDual<T> x)
+    {
+        var exp10 = Dual<T>.Exp10(x._d0);
+        return new(exp10, Real.Ln10 * x._d1 * exp10);
+    }
+
+    // Hyperbolic functions
+
+    public static HyperDual<T> Acosh(HyperDual<T> x)
+        => new(Dual<T>.Acosh(x._d0), x._d1 / (Dual<T>.Sqrt(x._d0 - T.One) * Dual<T>.Sqrt(x._d0 + T.One)));
+
+    public static HyperDual<T> Asinh(HyperDual<T> x)
+        => new(Dual<T>.Asinh(x._d0), x._d1 / Dual<T>.Sqrt(x._d0 * x._d0 + T.One));
+
+    public static HyperDual<T> Atanh(HyperDual<T> x)
+        => new(Dual<T>.Atanh(x._d0), x._d1 / (T.One - x._d0 * x._d0));
+
+    public static HyperDual<T> Cosh(HyperDual<T> x) => new(Dual<T>.Cosh(x._d0), x._d1 * Dual<T>.Sinh(x._d0));
+
+    public static HyperDual<T> Sinh(HyperDual<T> x) => new(Dual<T>.Sinh(x._d0), x._d1 * Dual<T>.Cosh(x._d0));
+
+    public static HyperDual<T> Tanh(HyperDual<T> x)
+    {
+        var u = T.One / Dual<T>.Cosh(x._d0);
+        return new(Dual<T>.Tanh(x._d0), x._d1 * u * u);
+    }
+
+    // Logarithmic functions
+
+    public static HyperDual<T> Ln(HyperDual<T> x) => new(Dual<T>.Ln(x._d0), x._d1 / x._d0);
+
+    public static HyperDual<T> Log(HyperDual<T> x, HyperDual<T> b) => Ln(x) / Ln(b);
+
+    public static HyperDual<T> Log2(HyperDual<T> x) => new(Dual<T>.Log2(x._d0), x._d1 / (Real.Ln2 * x._d0));
+
+    public static HyperDual<T> Log10(HyperDual<T> x) => new(Dual<T>.Log10(x._d0), x._d1 / (Real.Ln10 * x._d0));
+
+    // Power functions
+
+    public static HyperDual<T> Pow(HyperDual<T> x, HyperDual<T> y) => Exp(y * Ln(x));
+
+    // Root functions
+
+    public static HyperDual<T> Cbrt(HyperDual<T> x)
+    {
+        var cbrt = Dual<T>.Cbrt(x._d0);
+        return new(cbrt, x._d1 / (3.0 * cbrt * cbrt));
+    }
+
+    public static HyperDual<T> Root(HyperDual<T> x, HyperDual<T> n) => Exp(Ln(x) / n);
+
+    public static HyperDual<T> Sqrt(HyperDual<T> x)
+    {
+        var sqrt = Dual<T>.Sqrt(x._d0);
+        return new(sqrt, 0.5 * x._d1 / sqrt);
+    }
+
+    // Trigonometric functions
+
+    public static HyperDual<T> Acos(HyperDual<T> x) => new(Dual<T>.Acos(x._d0), -x._d1 / Dual<T>.Sqrt(T.One - x._d0 * x._d0));
+
+    public static HyperDual<T> Asin(HyperDual<T> x) => new(Dual<T>.Asin(x._d0), x._d1 / Dual<T>.Sqrt(T.One - x._d0 * x._d0));
+
+    public static HyperDual<T> Atan(HyperDual<T> x) => new(Dual<T>.Atan(x._d0), x._d1 / (T.One + x._d0 * x._d0));
+
+    /// <inheritdoc cref="IReal{T}.Atan2(T, T)"/>
+    public static HyperDual<Real> Atan2(HyperDual<Real> y, HyperDual<Real> x)
+    {
+        var u = Real.One / (x._d0 * x._d0 + y._d0 * y._d0);
+        return new(Dual<T>.Atan2(y._d0, x._d0), (y._d1 * x._d0 - x._d1 * y._d0) * u);
+    }
+
+    public static HyperDual<T> Cos(HyperDual<T> x) => new(Dual<T>.Cos(x._d0), -x._d1 * Dual<T>.Sin(x._d0));
+
+    public static HyperDual<T> Sin(HyperDual<T> x) => new(Dual<T>.Sin(x._d0), x._d1 * Dual<T>.Cos(x._d0));
+
+    public static HyperDual<T> Tan(HyperDual<T> x)
+    {
+        var sec = T.One / Dual<T>.Cos(x._d0);
+        return new(Dual<T>.Tan(x._d0), x._d1 * sec * sec);
+    }
+
+    //
+    // Custom operations
+    //
+
+    /// <summary>Perform forward-mode autodiff using a custom unary operation.</summary>
+    /// <param name="x">A hyper-dual number</param>
+    /// <param name="f">A function</param>
+    /// <param name="df">The derivative of the function</param>
+    /// <returns>A hyper-dual number</returns>
+    public static HyperDual<T> CustomOperation(HyperDual<T> x, Func<Dual<T>, Dual<T>> f, Func<Dual<T>, Dual<T>> df)
+        => new(f(x._d0), x._d1 * df(x._d0));
+
+    /// <summary>Perform forward-mode autodiff using a custom binary operation.</summary>
+    /// <param name="x">A variable</param>
+    /// <param name="y">A variable</param>
+    /// <param name="f">A function</param>
+    /// <param name="dfx">The derivative of the function with respect to the left variable</param>
+    /// <param name="dfy">The derivative of the function with respect to the right variable</param>
+    /// <returns>A hyper-dual number</returns>
+    public static HyperDual<T> CustomOperation(
+        HyperDual<T> x,
+        HyperDual<T> y,
+        Func<Dual<T>, Dual<T>, Dual<T>> f,
+        Func<Dual<T>, Dual<T>, Dual<T>> dfx,
+        Func<Dual<T>, Dual<T>, Dual<T>> dfy)
+        => new(f(x._d0, y._d0), dfy(x._d0, y._d0) * x._d1 + dfx(x._d0, y._d1) * y._d1);
+
+    //
+    // Dual vector creation
+    //
+
+    public static DualVector3<HyperDual<T>, T> CreateDualVector(HyperDual<T> x1Seed, HyperDual<T> x2Seed, HyperDual<T> x3Seed)
+        => new(x1Seed, x2Seed, x3Seed);
+}
diff --git a/src/Mathematics.NET/AutoDiff/IDual.cs b/src/Mathematics.NET/AutoDiff/IDual.cs
new file mode 100644
index 00000000..fc194327
--- /dev/null
+++ b/src/Mathematics.NET/AutoDiff/IDual.cs
@@ -0,0 +1,77 @@
+// <copyright file="IDual.cs" company="Mathematics.NET">
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+// </copyright>
+
+using Mathematics.NET.Core.Operations;
+using Mathematics.NET.Core.Relations;
+
+namespace Mathematics.NET.AutoDiff;
+
+/// <summary>Defines support for dual numbers</summary>
+/// <typeparam name="T">The type that implements the interface</typeparam>
+/// <typeparam name="U">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
+public interface IDual<T, U>
+    : IAdditionOperation<T, T>,
+      IDivisionOperation<T, T>,
+      IMultiplicationOperation<T, T>,
+      INegationOperation<T, T>,
+      ISubtractionOperation<T, T>,
+      IFormattable,
+      IEqualityRelation<T, bool>,
+      IEquatable<T>,
+      IDifferentiableFunctions<T>
+    where T : IDual<T, U>
+    where U : IComplex<U>, IDifferentiableFunctions<U>
+{
+    /// <summary>Represents the primal part of the dual number</summary>
+    U D0 { get; }
+
+    /// <summary>Represents the tangent part of the dual number</summary>
+    U D1 { get; }
+
+    /// <summary>Create a vector from seed values.</summary>
+    /// <param name="x1Seed">The first seed value</param>
+    /// <param name="x2Seed">The second seed value</param>
+    /// <param name="x3Seed">The third seed value</param>
+    /// <returns>A dual vector of length three</returns>
+    static abstract DualVector3<T, U> CreateDualVector(T x1Seed, T x2Seed, T x3Seed);
+
+    /// <summary>Create an instance of the type with a specified value.</summary>
+    /// <param name="value">A value</param>
+    /// <returns>An instance of the type</returns>
+    static abstract T CreateVariable(U value);
+
+    /// <summary>Create an instance of the type with a specified value and seed.</summary>
+    /// <param name="value">A value</param>
+    /// <param name="seed">A seed</param>
+    /// <returns>An instance of the type</returns>
+    static abstract T CreateVariable(U value, U seed);
+
+    /// <summary>Create a seeded instance of this type</summary>
+    /// <param name="seed">The seed value</param>
+    /// <returns>A seeded value</returns>
+    T WithSeed(U seed);
+}
diff --git a/src/Mathematics.NET/Core/IDifferentiableFunctions.cs b/src/Mathematics.NET/Core/IDifferentiableFunctions.cs
index 5c2a10bc..f1e12817 100644
--- a/src/Mathematics.NET/Core/IDifferentiableFunctions.cs
+++ b/src/Mathematics.NET/Core/IDifferentiableFunctions.cs
@@ -29,7 +29,7 @@ namespace Mathematics.NET.Core;
 
 /// <summary>Defines support for common differentiable functions</summary>
 /// <typeparam name="T">The type that implements the interface</typeparam>
-public interface IDifferentiableFunctions<T> : IComplex<T>
+public interface IDifferentiableFunctions<T>
     where T : IDifferentiableFunctions<T>
 {
     //
diff --git a/src/Mathematics.NET/Mathematics.NET.csproj b/src/Mathematics.NET/Mathematics.NET.csproj
index db00fb91..643d7acd 100644
--- a/src/Mathematics.NET/Mathematics.NET.csproj
+++ b/src/Mathematics.NET/Mathematics.NET.csproj
@@ -7,7 +7,7 @@
     <Nullable>enable</Nullable>
     <Platforms>x64</Platforms>
     <Title>Mathematics.NET</Title>
-    <Version>0.1.0-alpha.8</Version>
+    <Version>0.1.0-alpha.9</Version>
     <PackageIcon>mathematics.net.png</PackageIcon>
     <Authors>Hamlet Tanyavong</Authors>
     <Description>Mathematics.NET is a C# class library that provides tools for solving mathematical problems. Included are custom types for real, complex, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors. The library also contains methods for performing forward and reverse-mode automatic differentiation.</Description>
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfDualOfRealTests.cs
similarity index 65%
rename from tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs
rename to tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfDualOfRealTests.cs
index 7ab90b59..5db3d7d6 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfDualOfRealTests.cs
@@ -27,23 +27,22 @@
 
 using Mathematics.NET.AutoDiff;
 using Mathematics.NET.LinearAlgebra;
-using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
 namespace Mathematics.NET.Tests.AutoDiff;
 
 [TestClass]
 [TestCategory("AutoDiff"), TestCategory("Real Number")]
-public sealed class DualVector3OfRealTests
+public sealed class DualVector3OfDualOfRealTests
 {
     [TestMethod]
     [TestCategory("Vector Calculus")]
     [DataRow(1.23, 0.66, 2.34, 1.954144178335244, -1.142124546272508, 0.820964086423733)]
     public void Curl_VectorField_ReturnsCurl(double x, double y, double z, double expectedX, double expectedY, double expectedZ)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
 
-        var actual = DualVector3<Real>.Curl(FX, FY, FZ, u);
+        var actual = DualVector3<Dual<Real>, Real>.Curl(FX, FY, FZ, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -53,10 +52,10 @@ public void Curl_VectorField_ReturnsCurl(double x, double y, double z, double ex
     [DataRow(0.23, 1.57, -1.71, 1.23, 0.66, 2.34, -0.801549048972843)]
     public void DirectionalDerivative_ScalarFunctionAndDirection_ReturnsDirectionalDerivative(double vx, double vy, double vz, double x, double y, double z, double expected)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Vector3<Real> v = new(vx, vy, vz);
 
-        var actual = DualVector3<Real>.DirectionalDerivative(v, F, u);
+        var actual = DualVector3<Dual<Real>, Real>.DirectionalDerivative(v, F, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -66,9 +65,9 @@ public void DirectionalDerivative_ScalarFunctionAndDirection_ReturnsDirectionalD
     [DataRow(1.23, 0.66, 2.34, 0.3987010509910668)]
     public void Divergence_VectorField_ReturnsDivergence(double x, double y, double z, double expected)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
 
-        var actual = DualVector3<Real>.Divergence(FX, FY, FZ, u);
+        var actual = DualVector3<Dual<Real>, Real>.Divergence(FX, FY, FZ, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -78,25 +77,26 @@ public void Divergence_VectorField_ReturnsDivergence(double x, double y, double
     [DataRow(1.23, 0.66, 2.34, -0.824313594924351, -0.1302345967828155, 0.2382974299363869)]
     public void Gradient_ScalarFunction_ReturnsGradient(double x, double y, double z, double expectedX, double expectedY, double expectedZ)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
 
-        var actual = DualVector3<Real>.Gradient(F, u);
+        var actual = DualVector3<Dual<Real>, Real>.Gradient(F, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
 
     [TestMethod]
     [TestCategory("Vector Calculus")]
-    public void Jacobian_R3VectorFunction_ReturnsJacobian()
+    [DataRow(1.23, 0.66, 2.34)]
+    public void Jacobian_R3VectorFunction_ReturnsJacobian(double x, double y, double z)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Matrix3x3<Real> expected = new(
             0.775330615737715, -0.5778557672605755, 0.3080621020764366,
             0.2431083191631576, 0.2431083191631576, 0.2431083191631576,
             1.450186648348945, 2.197252497498402, -0.6197378839098056);
 
-        var actual = DualVector3<Real>.Jacobian(FX, FY, FZ, u);
+        var actual = DualVector3<Dual<Real>, Real>.Jacobian(FX, FY, FZ, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -106,11 +106,11 @@ public void Jacobian_R3VectorFunction_ReturnsJacobian()
     [DataRow(1.23, 0.66, 2.34, 0.23, 1.57, -1.71, -1.255693707530136, 0.0218797487246842, 4.842981131678516)]
     public void JVP_R3VectorFunctionAndVector_ReturnsJVP(double x, double y, double z, double vx, double vy, double vz, double expectedX, double expectedY, double expectedZ)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Vector3<Real> v = new(vx, vy, vz);
         Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
 
-        var actual = DualVector3<Real>.JVP(FX, FY, FZ, u, v);
+        var actual = DualVector3<Dual<Real>, Real>.JVP(FX, FY, FZ, u, v);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -120,11 +120,11 @@ public void JVP_R3VectorFunctionAndVector_ReturnsJVP(double x, double y, double
     [DataRow(0.23, 1.57, -1.71, 1.23, 0.66, 2.34, -1.919813065970865, -3.508528536106042, 1.512286126049506)]
     public void VJP_VectorAndR3VectorFunction_ReturnsVJP(double vx, double vy, double vz, double x, double y, double z, double expectedX, double expectedY, double expectedZ)
     {
-        DualVector3<Real> u = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+        DualVector3<Dual<Real>, Real> u = new(Dual<Real>.CreateVariable(x), Dual<Real>.CreateVariable(y), Dual<Real>.CreateVariable(z));
         Vector3<Real> v = new(vx, vy, vz);
         Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
 
-        var actual = DualVector3<Real>.VJP(v, FX, FY, FZ, u);
+        var actual = DualVector3<Dual<Real>, Real>.VJP(v, FX, FY, FZ, u);
 
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
@@ -133,15 +133,23 @@ public void VJP_VectorAndR3VectorFunction_ReturnsVJP(double vx, double vy, doubl
     // Helpers
     //
 
-    private static Dual<Real> F(DualVector3<Real> x)
-        => Cos(x.X1) / ((x.X1 + x.X2) * Sin(x.X3));
-
-    private static Dual<Real> FX(DualVector3<Real> x)
-        => Sin(x.X1) * (Cos(x.X2) + Sqrt(x.X3));
-
-    private static Dual<Real> FY(DualVector3<Real> x)
-        => Sqrt(x.X1 + x.X2 + x.X3);
-
-    private static Dual<Real> FZ(DualVector3<Real> x)
-        => Sinh(Exp(x.X1) * x.X2 / x.X3);
+    private static T F<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Cos(x.X1) / ((x.X1 + x.X2) * T.Sin(x.X3));
+
+    private static T FX<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sin(x.X1) * (T.Cos(x.X2) + T.Sqrt(x.X3));
+
+    private static T FY<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sqrt(x.X1 + x.X2 + x.X3);
+
+    private static T FZ<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sinh(T.Exp(x.X1) * x.X2 / x.X3);
 }
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfHyperDualOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfHyperDualOfRealTests.cs
new file mode 100644
index 00000000..fc2e7030
--- /dev/null
+++ b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfHyperDualOfRealTests.cs
@@ -0,0 +1,167 @@
+// <copyright file="DualVector3OfHyperDualOfRealTests.cs" company="Mathematics.NET">
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+// </copyright>
+
+using Mathematics.NET.AutoDiff;
+using Mathematics.NET.LinearAlgebra;
+
+namespace Mathematics.NET.Tests.AutoDiff;
+
+[TestClass]
+[TestCategory("AutoDiff"), TestCategory("Real Number")]
+public sealed class DualVector3OfHyperDualOfRealTests
+{
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34, 1.954144178335244, -1.142124546272508, 0.820964086423733)]
+    public void Curl_VectorField_ReturnsCurl(double x, double y, double z, double expectedX, double expectedY, double expectedZ)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.Curl(FX, FY, FZ, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(0.23, 1.57, -1.71, 1.23, 0.66, 2.34, -0.801549048972843)]
+    public void DirectionalDerivative_ScalarFunctionAndDirection_ReturnsDirectionalDerivative(double vx, double vy, double vz, double x, double y, double z, double expected)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Vector3<Real> v = new(vx, vy, vz);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.DirectionalDerivative(v, F, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34, 0.3987010509910668)]
+    public void Divergence_VectorField_ReturnsDivergence(double x, double y, double z, double expected)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+
+        var actual = DualVector3<HyperDual<Real>, Real>.Divergence(FX, FY, FZ, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34, -0.824313594924351, -0.1302345967828155, 0.2382974299363869)]
+    public void Gradient_ScalarFunction_ReturnsGradient(double x, double y, double z, double expectedX, double expectedY, double expectedZ)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.Gradient(F, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34)]
+    public void Jacobian_R3VectorFunction_ReturnsJacobian(double x, double y, double z)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Matrix3x3<Real> expected = new(
+            0.775330615737715, -0.5778557672605755, 0.3080621020764366,
+            0.2431083191631576, 0.2431083191631576, 0.2431083191631576,
+            1.450186648348945, 2.197252497498402, -0.6197378839098056);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.Jacobian(FX, FY, FZ, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34, 0.23, 1.57, -1.71, -1.255693707530136, 0.0218797487246842, 4.842981131678516)]
+    public void JVP_R3VectorFunctionAndVector_ReturnsJVP(double x, double y, double z, double vx, double vy, double vz, double expectedX, double expectedY, double expectedZ)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Vector3<Real> v = new(vx, vy, vz);
+        Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.JVP(FX, FY, FZ, u, v);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(1.23, 0.66, 2.34, 1.471507039061705)]
+    public void Laplacian_ScalarFunction_ReturnsLaplacian(double x, double y, double z, double expected)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+
+        var actual = DualVector3<HyperDual<Real>, Real>.Laplacian(F, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    [DataRow(0.23, 1.57, -1.71, 1.23, 0.66, 2.34, -1.919813065970865, -3.508528536106042, 1.512286126049506)]
+    public void VJP_VectorAndR3VectorFunction_ReturnsVJP(double vx, double vy, double vz, double x, double y, double z, double expectedX, double expectedY, double expectedZ)
+    {
+        DualVector3<HyperDual<Real>, Real> u = new(HyperDual<Real>.CreateVariable(x), HyperDual<Real>.CreateVariable(y), HyperDual<Real>.CreateVariable(z));
+        Vector3<Real> v = new(vx, vy, vz);
+        Vector3<Real> expected = new(expectedX, expectedY, expectedZ);
+
+        var actual = DualVector3<HyperDual<Real>, Real>.VJP(v, FX, FY, FZ, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    //
+    // Helpers
+    //
+
+    private static T F<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Cos(x.X1) / ((x.X1 + x.X2) * T.Sin(x.X3));
+
+    private static T FX<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sin(x.X1) * (T.Cos(x.X2) + T.Sqrt(x.X3));
+
+    private static T FY<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sqrt(x.X1 + x.X2 + x.X3);
+
+    private static T FZ<T, U>(DualVector3<T, U> x)
+        where T : IDual<T, U>
+        where U : IComplex<U>, IDifferentiableFunctions<U>
+        => T.Sinh(T.Exp(x.X1) * x.X2 / x.X3);
+}
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
index 258038d9..82ac2a21 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
@@ -346,10 +346,10 @@ public void CustomOperation_Binary_ReturnsGradient(double left, double right, do
             (y, x) => -y.Value * u,
             (y, x) => Real.Zero); // Not of interest
 
-        _tape.ReverseAccumulation(out ReadOnlySpan<Real> actual);
-
         Real[] expected = [expectedLeft, expectedRight];
 
+        _tape.ReverseAccumulation(out ReadOnlySpan<Real> actual);
+
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
 
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/TapeExtensionsOfRealTests.cs
similarity index 90%
rename from tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs
rename to tests/Mathematics.NET.Tests/AutoDiff/TapeExtensionsOfRealTests.cs
index cf9a4ea1..770bbd71 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/TapeExtensionsOfRealTests.cs
@@ -1,4 +1,4 @@
-// <copyright file="GradientTapeExtensionsOfRealTests.cs" company="Mathematics.NET">
+// <copyright file="TapeExtensionsOfRealTests.cs" company="Mathematics.NET">
 // Mathematics.NET
 // https://github.com/HamletTanyavong/Mathematics.NET
 //
@@ -32,11 +32,11 @@ namespace Mathematics.NET.Tests.AutoDiff;
 
 [TestClass]
 [TestCategory("AutoDiff"), TestCategory("Real Number")]
-public sealed class GradientTapeExtensionsOfRealTests
+public sealed class TapeExtensionsOfRealTests
 {
     private GradientTape<Real> _tape;
 
-    public GradientTapeExtensionsOfRealTests()
+    public TapeExtensionsOfRealTests()
     {
         _tape = new();
     }
@@ -92,6 +92,22 @@ public void Gradient_ScalarFunction_ReturnsGradient(double x, double y, double z
         Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
     }
 
+    [TestMethod]
+    [TestCategory("Vector Calculus")]
+    public void Hessian_ScalarFunction_ReturnsHessian()
+    {
+        HessianTape<Real> tape = new();
+        var u = tape.CreateVariableVector(1.23, 0.66, 2.34);
+        Matrix3x3<Real> expected = new(
+            0.6261461305189455, 0.5050519532842152, -0.7980381386329245,
+            0.5050519532842152, 0.1378143881299635, -0.1260832962626385,
+            -0.7980381386329245, -0.1260832962626385, 0.707546520412796);
+
+        var actual = tape.Hessian(F, u);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
     [TestMethod]
     [TestCategory("Vector Calculus")]
     public void Jacobian_R3VectorFunction_ReturnsJacobian()