From be8c884029c9dbb01f4e716e1af03be6b3b63fec Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 7 Nov 2023 01:01:53 -0600
Subject: [PATCH 01/18] Update GradientTape.cs

- Add support for cancelling print node operation
---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 9d25078b..de8b8651 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -86,7 +86,7 @@ public GradientTape()
 
     /// <summary>Print the nodes of the gradient tape to the console</summary>
     /// <param name="limit">The total number of nodes to print</param>
-    public void PrintNodes(int limit = 100)
+    public void PrintNodes(CancellationToken cancellationToken, int limit = 100)
     {
         const string tab = "    ";
 
@@ -96,21 +96,35 @@ public void PrintNodes(int limit = 100)
         int i = 0;
         while (i < Math.Min(_variableCount, limit))
         {
+            CheckForCancellation(cancellationToken);
             node = nodeSpan[i];
             Console.WriteLine($"Root Node {i}:");
             Console.WriteLine($"{tab}Weights: [{node.DX}, {node.DY}]");
             Console.WriteLine($"{tab}Parents: [{node.PX}, {node.PY}]");
             i++;
         }
+
+        CheckForCancellation(cancellationToken);
         Console.WriteLine();
+
         while (i < Math.Min(nodeSpan.Length, limit))
         {
+            CheckForCancellation(cancellationToken);
             node = nodeSpan[i];
             Console.WriteLine($"Node {i}:");
             Console.WriteLine($"{tab}Weights: [{node.DX}, {node.DY}]");
             Console.WriteLine($"{tab}Parents: [{node.PX}, {node.PY}]");
             i++;
         }
+
+        static void CheckForCancellation(CancellationToken cancellationToken)
+        {
+            if (cancellationToken.IsCancellationRequested)
+            {
+                Console.WriteLine("Print node operation cancelled");
+                cancellationToken.ThrowIfCancellationRequested();
+            }
+        }
     }
 
     /// <summary>Perform reverse accumulation on the gradient tape and output the resulting gradients</summary>

From 5ef762e23fa0c6066c7ddb3b58d71b80ffcf61ab Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 7 Nov 2023 16:04:55 -0600
Subject: [PATCH 02/18] Make index field internal

---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 78 ++++++++++----------
 src/Mathematics.NET/AutoDiff/Variable.cs     | 14 ++--
 2 files changed, 46 insertions(+), 46 deletions(-)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index de8b8651..3512a853 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -169,94 +169,94 @@ public Variable CreateVariable(Real seed)
 
     public Variable Add(Variable x, Variable y)
     {
-        _nodes.Add(new(Real.One, Real.One, x.Index, y.Index));
+        _nodes.Add(new(Real.One, Real.One, x._index, y._index));
         return new(_nodes.Count - 1, x.Value + y.Value);
     }
 
     public Variable Add(Real c, Variable x)
     {
-        _nodes.Add(new(Real.One, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c + x.Value);
     }
 
     public Variable Add(Variable x, Real c)
     {
-        _nodes.Add(new(Real.One, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value + c);
     }
 
     public Variable Divide(Variable x, Variable y)
     {
         var u = Real.One / y.Value;
-        _nodes.Add(new(Real.One / y.Value, -x.Value * u * u, x.Index, y.Index));
+        _nodes.Add(new(Real.One / y.Value, -x.Value * u * u, x._index, y._index));
         return new(_nodes.Count - 1, x.Value * u);
     }
 
     public Variable Divide(Real c, Variable x)
     {
         var u = Real.One / x.Value;
-        _nodes.Add(new(-c.Value * u * u, x.Index, _nodes.Count));
+        _nodes.Add(new(-c.Value * u * u, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value * u);
     }
 
     public Variable Divide(Variable x, Real c)
     {
         var u = Real.One / c.Value;
-        _nodes.Add(new(u, x.Index, _nodes.Count));
+        _nodes.Add(new(u, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value * u);
     }
 
     public Variable Modulo(Variable x, Variable y)
     {
-        _nodes.Add(new(Real.One, x.Value * Real.Floor(x.Value / y.Value), x.Index, y.Index));
+        _nodes.Add(new(Real.One, x.Value * Real.Floor(x.Value / y.Value), x._index, y._index));
         return new(_nodes.Count - 1, x.Value % y.Value);
     }
 
     public Variable Modulo(Real c, Variable x)
     {
-        _nodes.Add(new(c * Real.Floor(c / x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(c * Real.Floor(c / x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, c % x.Value);
     }
 
     public Variable Modulo(Variable x, Real c)
     {
-        _nodes.Add(new(Real.One, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value % c);
     }
 
     public Variable Multiply(Variable x, Variable y)
     {
-        _nodes.Add(new(y.Value, x.Value, x.Index, y.Index));
+        _nodes.Add(new(y.Value, x.Value, x._index, y._index));
         return new(_nodes.Count - 1, x.Value * y.Value);
     }
 
     public Variable Multiply(Real c, Variable x)
     {
-        _nodes.Add(new(c, x.Index, _nodes.Count));
+        _nodes.Add(new(c, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c * x.Value);
     }
 
     public Variable Multiply(Variable x, Real c)
     {
-        _nodes.Add(new(c, x.Index, _nodes.Count));
+        _nodes.Add(new(c, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value * c);
     }
 
     public Variable Subtract(Variable x, Variable y)
     {
-        _nodes.Add(new(Real.One, -Real.One, x.Index, y.Index));
+        _nodes.Add(new(Real.One, -Real.One, x._index, y._index));
         return new(_nodes.Count - 1, x.Value - y.Value);
     }
 
     public Variable Subtract(Real c, Variable x)
     {
-        _nodes.Add(new(-Real.One, x.Index, _nodes.Count));
+        _nodes.Add(new(-Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c - x.Value);
     }
 
     public Variable Subtract(Variable x, Real c)
     {
-        _nodes.Add(new(Real.One, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value - c);
     }
 
@@ -269,21 +269,21 @@ public Variable Subtract(Variable x, Real c)
     public Variable Exp(Variable x)
     {
         var exp = Real.Exp(x.Value);
-        _nodes.Add(new(exp, x.Index, _nodes.Count));
+        _nodes.Add(new(exp, x._index, _nodes.Count));
         return new(_nodes.Count - 1, exp);
     }
 
     public Variable Exp2(Variable x)
     {
         var exp2 = Real.Exp2(x.Value);
-        _nodes.Add(new(Real.Ln2 * exp2, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.Ln2 * exp2, x._index, _nodes.Count));
         return new(_nodes.Count - 1, exp2);
     }
 
     public Variable Exp10(Variable x)
     {
         var exp10 = Real.Exp10(x.Value);
-        _nodes.Add(new(Real.Ln10 * exp10, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.Ln10 * exp10, x._index, _nodes.Count));
         return new(_nodes.Count - 1, exp10);
     }
 
@@ -291,38 +291,38 @@ public Variable Exp10(Variable x)
 
     public Variable Acosh(Variable x)
     {
-        _nodes.Add(new(Real.One / (Complex.Sqrt(x.Value - Real.One) * Complex.Sqrt(x.Value + Real.One)).Re, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (Complex.Sqrt(x.Value - Real.One) * Complex.Sqrt(x.Value + Real.One)).Re, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Acosh(x.Value));
     }
 
     public Variable Asinh(Variable x)
     {
-        _nodes.Add(new(Real.One / Real.Sqrt(x.Value * x.Value + Real.One), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / Real.Sqrt(x.Value * x.Value + Real.One), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Asinh(x.Value));
     }
 
     public Variable Atanh(Variable x)
     {
-        _nodes.Add(new(Real.One / (Real.One - x.Value * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Atanh(x.Value));
     }
 
     public Variable Cosh(Variable x)
     {
-        _nodes.Add(new(Real.Sinh(x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.Sinh(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Cosh(x.Value));
     }
 
     public Variable Sinh(Variable x)
     {
-        _nodes.Add(new(Real.Cosh(x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.Cosh(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Sinh(x.Value));
     }
 
     public Variable Tanh(Variable x)
     {
         var u = Real.One / Real.Cosh(x.Value);
-        _nodes.Add(new(u * u, x.Index, _nodes.Count));
+        _nodes.Add(new(u * u, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Tanh(x.Value));
     }
 
@@ -330,26 +330,26 @@ public Variable Tanh(Variable x)
 
     public Variable Ln(Variable x)
     {
-        _nodes.Add(new(Real.One / x.Value, x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / x.Value, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Ln(x.Value));
     }
 
     public Variable Log(Variable x, Variable b)
     {
         var lnB = Real.Ln(b.Value);
-        _nodes.Add(new(Real.One / (x.Value * lnB), -Real.Ln(x.Value) / (b.Value * lnB * lnB), x.Index, b.Index));
+        _nodes.Add(new(Real.One / (x.Value * lnB), -Real.Ln(x.Value) / (b.Value * lnB * lnB), x._index, b._index));
         return new(_nodes.Count - 1, Real.Log(x.Value, b.Value));
     }
 
     public Variable Log2(Variable x)
     {
-        _nodes.Add(new(Real.One / (Real.Ln2 * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (Real.Ln2 * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Log2(x.Value));
     }
 
     public Variable Log10(Variable x)
     {
-        _nodes.Add(new(Real.One / (Real.Ln10 * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (Real.Ln10 * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Log10(x.Value));
     }
 
@@ -358,7 +358,7 @@ public Variable Log10(Variable x)
     public Variable Pow(Variable x, Variable y)
     {
         var pow = Real.Pow(x.Value, y.Value);
-        _nodes.Add(new(y.Value * Real.Pow(x.Value, y.Value - Real.One), Real.Ln(x.Value) * pow, x.Index, y.Index));
+        _nodes.Add(new(y.Value * Real.Pow(x.Value, y.Value - Real.One), Real.Ln(x.Value) * pow, x._index, y._index));
         return new(_nodes.Count - 1, pow);
     }
 
@@ -367,21 +367,21 @@ public Variable Pow(Variable x, Variable y)
     public Variable Cbrt(Variable x)
     {
         var cbrt = Real.Cbrt(x.Value);
-        _nodes.Add(new(Real.One / (3.0 * cbrt * cbrt), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (3.0 * cbrt * cbrt), x._index, _nodes.Count));
         return new(_nodes.Count - 1, cbrt);
     }
 
     public Variable Root(Variable x, Variable n)
     {
         var root = Real.Root(x.Value, n.Value);
-        _nodes.Add(new(root / (n.Value * x.Value), -Real.Ln(x.Value) * root / (n.Value * n.Value), x.Index, n.Index));
+        _nodes.Add(new(root / (n.Value * x.Value), -Real.Ln(x.Value) * root / (n.Value * n.Value), x._index, n._index));
         return new(_nodes.Count - 1, root);
     }
 
     public Variable Sqrt(Variable x)
     {
         var sqrt = Real.Sqrt(x.Value);
-        _nodes.Add(new(0.5 / sqrt, x.Index, _nodes.Count));
+        _nodes.Add(new(0.5 / sqrt, x._index, _nodes.Count));
         return new(_nodes.Count - 1, sqrt);
     }
 
@@ -389,45 +389,45 @@ public Variable Sqrt(Variable x)
 
     public Variable Acos(Variable x)
     {
-        _nodes.Add(new(-Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(-Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Acos(x.Value));
     }
 
     public Variable Asin(Variable x)
     {
-        _nodes.Add(new(Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Asin(x.Value));
     }
 
     public Variable Atan(Variable x)
     {
-        _nodes.Add(new(Real.One / (Real.One + x.Value * x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.One / (Real.One + x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Atan(x.Value));
     }
 
     public Variable Atan2(Variable y, Variable x)
     {
         var u = Real.One / (x.Value * x.Value + y.Value * y.Value);
-        _nodes.Add(new(-x.Value * u, y.Value * u, y.Index, x.Index));
+        _nodes.Add(new(-x.Value * u, y.Value * u, y._index, x._index));
         return new(_nodes.Count - 1, Real.Atan2(y.Value, x.Value));
     }
 
     public Variable Cos(Variable x)
     {
-        _nodes.Add(new(-Real.Sin(x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(-Real.Sin(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Cos(x.Value));
     }
 
     public Variable Sin(Variable x)
     {
-        _nodes.Add(new(Real.Cos(x.Value), x.Index, _nodes.Count));
+        _nodes.Add(new(Real.Cos(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Sin(x.Value));
     }
 
     public Variable Tan(Variable x)
     {
         var sec = Real.One / Real.Cos(x.Value);
-        _nodes.Add(new(sec * sec, x.Index, _nodes.Count));
+        _nodes.Add(new(sec * sec, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Tan(x.Value));
     }
 }
diff --git a/src/Mathematics.NET/AutoDiff/Variable.cs b/src/Mathematics.NET/AutoDiff/Variable.cs
index 18f68c65..87a783c6 100644
--- a/src/Mathematics.NET/AutoDiff/Variable.cs
+++ b/src/Mathematics.NET/AutoDiff/Variable.cs
@@ -30,15 +30,15 @@ namespace Mathematics.NET.AutoDiff;
 /// <summary>Represents a variable used in reverse-mode automatic differentiation</summary>
 public readonly record struct Variable
 {
-    internal Variable(int index, Real value)
-    {
-        Index = index;
-        Value = value;
-    }
-
     /// <summary>The index of the variable</summary>
-    public readonly int Index;
+    internal readonly int _index;
 
     /// <summary>The value of the variable</summary>
     public readonly Real Value;
+
+    internal Variable(int index, Real value)
+    {
+        _index = index;
+        Value = value;
+    }
 }

From c29e05e8eabffc625f038748b6c79f3dbe69b4d0 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 7 Nov 2023 16:05:25 -0600
Subject: [PATCH 03/18] Update MatrixMultiplyByScalarBenchmarks.cs

---
 .../LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs  | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/Mathematics.NET.Benchmarks/LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs b/tests/Mathematics.NET.Benchmarks/LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs
index f8fbc156..41be8eb4 100644
--- a/tests/Mathematics.NET.Benchmarks/LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs
+++ b/tests/Mathematics.NET.Benchmarks/LinearAlgebra/MatrixMultiplyByScalarBenchmarks.cs
@@ -38,8 +38,8 @@ public class MatrixMultiplyByScalarBenchmarks
 {
     public int Rows { get; set; }
     public int Cols { get; set; }
-    public required NET.Core.Complex[,] MatrixOne { get; set; }
-    public required NET.Core.Complex[,] MatrixTwo { get; set; }
+    public required Complex[,] MatrixOne { get; set; }
+    public required Complex[,] MatrixTwo { get; set; }
 
     [GlobalSetup]
     public void GlobalSetup()
@@ -47,7 +47,7 @@ public void GlobalSetup()
         Rows = 100;
         Cols = 100;
 
-        MatrixOne = new NET.Core.Complex[Rows, Cols];
+        MatrixOne = new Complex[Rows, Cols];
 
         for (int i = 0; i < Rows; i++)
         {
@@ -75,7 +75,7 @@ public void MultiplyByScalarNaive()
     [Benchmark]
     public void MultiplyByScalarParallel()
     {
-        Memory2D<NET.Core.Complex> matrixAsMemory = MatrixTwo;
-        ParallelHelper.ForEach(matrixAsMemory, new MultiplyByScalarAction<NET.Core.Complex>(Real.Pi));
+        Memory2D<Complex> matrixAsMemory = MatrixTwo;
+        ParallelHelper.ForEach(matrixAsMemory, new MultiplyByScalarAction<Complex>(Real.Pi));
     }
 }

From 030907c717b88103e28a10d100220044d658e458 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 8 Nov 2023 02:07:57 -0600
Subject: [PATCH 04/18] Update GradientTape.cs

---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 3512a853..4fcc1837 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -133,6 +133,11 @@ static void CheckForCancellation(CancellationToken cancellationToken)
     [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
     public void ReverseAccumulation(out ReadOnlySpan<Real> gradients, double seed = 1.0)
     {
+        if (_variableCount == 0)
+        {
+            throw new Exception("Gradient tape contains no nodes");
+        }
+
         ReadOnlySpan<Node> nodesAsSpan = CollectionsMarshal.AsSpan(_nodes);
         ref var start = ref MemoryMarshal.GetReference(nodesAsSpan);
 

From a3776410e4002a6259c192c5eb46fbcd03c03032 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 8 Nov 2023 17:13:48 -0600
Subject: [PATCH 05/18] Update GradientTape.cs

- Minor changes
---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 40 ++++++++++----------
 1 file changed, 21 insertions(+), 19 deletions(-)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 4fcc1837..75784bb1 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -66,6 +66,8 @@ namespace Mathematics.NET.AutoDiff;
 /// <summary>Represents a gradient tape</summary>
 public record class GradientTape
 {
+    // TODO: Measure performance with Stack<Node> instead of List<Node>
+    // TODO: Consider using array pools or something similar
     private List<Node> _nodes;
     private int _variableCount;
 
@@ -74,17 +76,28 @@ public GradientTape()
         _nodes = new();
     }
 
-    /// <summary>Get the number of nodes on the gradient tape</summary>
+    /// <summary>Get the number of nodes on the gradient tape.</summary>
     public int NodeCount => _nodes.Count;
 
-    /// <summary>Get the number of variables that are being tracked</summary>
+    /// <summary>Get the number of variables that are being tracked.</summary>
     public int VariableCount => _variableCount;
 
     //
     // Methods
     //
 
-    /// <summary>Print the nodes of the gradient tape to the console</summary>
+    /// <summary>Create a variable for the gradient tape to track.</summary>
+    /// <param name="seed">A seed value</param>
+    /// <returns>A variable</returns>
+    public Variable CreateVariable(Real seed)
+    {
+        _nodes.Add(new(_variableCount));
+        Variable variable = new(_variableCount++, seed);
+        return variable;
+    }
+
+    /// <summary>Print the nodes of the gradient tape to the console.</summary>
+    /// <param name="cancellationToken">A cancellation token</param>
     /// <param name="limit">The total number of nodes to print</param>
     public void PrintNodes(CancellationToken cancellationToken, int limit = 100)
     {
@@ -127,7 +140,7 @@ static void CheckForCancellation(CancellationToken cancellationToken)
         }
     }
 
-    /// <summary>Perform reverse accumulation on the gradient tape and output the resulting gradients</summary>
+    /// <summary>Perform reverse accumulation on the gradient tape and output the resulting gradients.</summary>
     /// <param name="gradients">The gradients</param>
     /// <param name="seed">A seed value</param>
     [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
@@ -135,13 +148,13 @@ public void ReverseAccumulation(out ReadOnlySpan<Real> gradients, double seed =
     {
         if (_variableCount == 0)
         {
-            throw new Exception("Gradient tape contains no nodes");
+            throw new Exception("Gradient tape contains no root nodes");
         }
 
-        ReadOnlySpan<Node> nodesAsSpan = CollectionsMarshal.AsSpan(_nodes);
-        ref var start = ref MemoryMarshal.GetReference(nodesAsSpan);
+        ReadOnlySpan<Node> nodes = CollectionsMarshal.AsSpan(_nodes);
+        ref var start = ref MemoryMarshal.GetReference(nodes);
 
-        var length = nodesAsSpan.Length;
+        var length = nodes.Length;
         Span<Real> partialGradients = new Real[length];
         partialGradients[length - 1] = seed;
 
@@ -157,17 +170,6 @@ public void ReverseAccumulation(out ReadOnlySpan<Real> gradients, double seed =
         gradients = partialGradients[.._variableCount];
     }
 
-    /// <summary>Create a variable for the gradient tape to track</summary>
-    /// <param name="seed">A seed value</param>
-    /// <returns>A variable</returns>
-    public Variable CreateVariable(Real seed)
-    {
-        _nodes.Add(new(_variableCount));
-        Variable variable = new(_variableCount, seed);
-        _variableCount++;
-        return variable;
-    }
-
     //
     // Basic operations
     //

From 48378345b5326e368b8ae49b1c9e9de59ce0d734 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 13:39:27 -0600
Subject: [PATCH 06/18] Update Variable.cs

- Override ToString
---
 src/Mathematics.NET/AutoDiff/Variable.cs | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/Mathematics.NET/AutoDiff/Variable.cs b/src/Mathematics.NET/AutoDiff/Variable.cs
index 87a783c6..a95c5ab0 100644
--- a/src/Mathematics.NET/AutoDiff/Variable.cs
+++ b/src/Mathematics.NET/AutoDiff/Variable.cs
@@ -41,4 +41,6 @@ internal Variable(int index, Real value)
         _index = index;
         Value = value;
     }
+
+    public override string ToString() => Value.ToString();
 }

From a201a297f3133ccd88561ec30386b0ae42e4aa55 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 14:15:16 -0600
Subject: [PATCH 07/18] Update GradientTape.cs

- Add documentation comments
---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 83 ++++++++++++++++++++
 1 file changed, 83 insertions(+)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 75784bb1..62b3a3b3 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -174,24 +174,40 @@ public void ReverseAccumulation(out ReadOnlySpan<Real> gradients, double seed =
     // Basic operations
     //
 
+    /// <summary>Add two variables</summary>
+    /// <param name="x">The first variable</param>
+    /// <param name="y">The second variable</param>
+    /// <returns>A variable</returns>
     public Variable Add(Variable x, Variable y)
     {
         _nodes.Add(new(Real.One, Real.One, x._index, y._index));
         return new(_nodes.Count - 1, x.Value + y.Value);
     }
 
+    /// <summary>Add a real value and a variable</summary>
+    /// <param name="c">A real value</param>
+    /// <param name="x">A variable</param>
+    /// <returns>A variable</returns>
     public Variable Add(Real c, Variable x)
     {
         _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c + x.Value);
     }
 
+    /// <summary>Add a variable and a real value</summary>
+    /// <param name="x">A variable</param>
+    /// <param name="c">A real value</param>
+    /// <returns>A variable</returns>
     public Variable Add(Variable x, Real c)
     {
         _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value + c);
     }
 
+    /// <summary>Divide two variables</summary>
+    /// <param name="x">A dividend</param>
+    /// <param name="y">A divisor</param>
+    /// <returns>A variable</returns>
     public Variable Divide(Variable x, Variable y)
     {
         var u = Real.One / y.Value;
@@ -199,6 +215,10 @@ public Variable Divide(Variable x, Variable y)
         return new(_nodes.Count - 1, x.Value * u);
     }
 
+    /// <summary>Divide a real value by a variable</summary>
+    /// <param name="c">A real dividend</param>
+    /// <param name="x">A variable divisor</param>
+    /// <returns>A variable</returns>
     public Variable Divide(Real c, Variable x)
     {
         var u = Real.One / x.Value;
@@ -206,6 +226,10 @@ public Variable Divide(Real c, Variable x)
         return new(_nodes.Count - 1, x.Value * u);
     }
 
+    /// <summary>Divide a variable by a real value</summary>
+    /// <param name="x">A variable dividend</param>
+    /// <param name="c">A real divisor</param>
+    /// <returns>A variable</returns>
     public Variable Divide(Variable x, Real c)
     {
         var u = Real.One / c.Value;
@@ -213,54 +237,90 @@ public Variable Divide(Variable x, Real c)
         return new(_nodes.Count - 1, x.Value * u);
     }
 
+    /// <summary>Compute the modulo of a variable given a divisor</summary>
+    /// <param name="x">A dividend</param>
+    /// <param name="y">A divisor</param>
+    /// <returns><paramref name="x"/> mod <paramref name="y"/></returns>
     public Variable Modulo(Variable x, Variable y)
     {
         _nodes.Add(new(Real.One, x.Value * Real.Floor(x.Value / y.Value), x._index, y._index));
         return new(_nodes.Count - 1, x.Value % y.Value);
     }
 
+    /// <summary>Compute the modulo of a real value given a divisor</summary>
+    /// <param name="c">A real dividend</param>
+    /// <param name="x">A variable divisor</param>
+    /// <returns><paramref name="c"/> mod <paramref name="x"/></returns>
     public Variable Modulo(Real c, Variable x)
     {
         _nodes.Add(new(c * Real.Floor(c / x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, c % x.Value);
     }
 
+    /// <summary>Compute the modulo of a variable given a divisor</summary>
+    /// <param name="x">A variable dividend</param>
+    /// <param name="c">A real divisor</param>
+    /// <returns><paramref name="x"/> mod <paramref name="c"/></returns>
     public Variable Modulo(Variable x, Real c)
     {
         _nodes.Add(new(Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value % c);
     }
 
+    /// <summary>Multiply two variables</summary>
+    /// <param name="x">The first variable</param>
+    /// <param name="y">The second variable</param>
+    /// <returns>A variable</returns>
     public Variable Multiply(Variable x, Variable y)
     {
         _nodes.Add(new(y.Value, x.Value, x._index, y._index));
         return new(_nodes.Count - 1, x.Value * y.Value);
     }
 
+    /// <summary>Multiply a real value by a variable</summary>
+    /// <param name="c">A real value</param>
+    /// <param name="x">A variable</param>
+    /// <returns>A variable</returns>
     public Variable Multiply(Real c, Variable x)
     {
         _nodes.Add(new(c, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c * x.Value);
     }
 
+    /// <summary>Multiply a variable by a real value</summary>
+    /// <param name="x">A variable</param>
+    /// <param name="c">A real value</param>
+    /// <returns>A variable</returns>
     public Variable Multiply(Variable x, Real c)
     {
         _nodes.Add(new(c, x._index, _nodes.Count));
         return new(_nodes.Count - 1, x.Value * c);
     }
 
+    /// <summary>Subract two variables</summary>
+    /// <param name="x">The first variable</param>
+    /// <param name="y">The second variable</param>
+    /// <returns>A variable</returns>
     public Variable Subtract(Variable x, Variable y)
     {
         _nodes.Add(new(Real.One, -Real.One, x._index, y._index));
         return new(_nodes.Count - 1, x.Value - y.Value);
     }
 
+    /// <summary>Subtract a variable from a real value</summary>
+    /// <param name="c">A real value</param>
+    /// <param name="x">A variable</param>
+    /// <returns>A variable</returns>
     public Variable Subtract(Real c, Variable x)
     {
         _nodes.Add(new(-Real.One, x._index, _nodes.Count));
         return new(_nodes.Count - 1, c - x.Value);
     }
 
+    /// <summary>Subtract a real value from a variable</summary>
+    /// <param name="x">A variable</param>
+    /// <param name="c">A real value</param>
+    /// <returns>A variable</returns>
     public Variable Subtract(Variable x, Real c)
     {
         _nodes.Add(new(Real.One, x._index, _nodes.Count));
@@ -273,6 +333,7 @@ public Variable Subtract(Variable x, Real c)
 
     // Exponential functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp(T)"/>
     public Variable Exp(Variable x)
     {
         var exp = Real.Exp(x.Value);
@@ -280,6 +341,7 @@ public Variable Exp(Variable x)
         return new(_nodes.Count - 1, exp);
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp2(T)"/>
     public Variable Exp2(Variable x)
     {
         var exp2 = Real.Exp2(x.Value);
@@ -287,6 +349,7 @@ public Variable Exp2(Variable x)
         return new(_nodes.Count - 1, exp2);
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Exp10(T)"/>
     public Variable Exp10(Variable x)
     {
         var exp10 = Real.Exp10(x.Value);
@@ -296,36 +359,42 @@ public Variable Exp10(Variable x)
 
     // Hyperbolic functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Acosh(T)"/>
     public Variable Acosh(Variable x)
     {
         _nodes.Add(new(Real.One / (Complex.Sqrt(x.Value - Real.One) * Complex.Sqrt(x.Value + Real.One)).Re, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Acosh(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Asinh(T)"/>
     public Variable Asinh(Variable x)
     {
         _nodes.Add(new(Real.One / Real.Sqrt(x.Value * x.Value + Real.One), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Asinh(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Atanh(T)"/>
     public Variable Atanh(Variable x)
     {
         _nodes.Add(new(Real.One / (Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Atanh(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cosh(T)"/>
     public Variable Cosh(Variable x)
     {
         _nodes.Add(new(Real.Sinh(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Cosh(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sinh(T)"/>
     public Variable Sinh(Variable x)
     {
         _nodes.Add(new(Real.Cosh(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Sinh(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Tanh(T)"/>
     public Variable Tanh(Variable x)
     {
         var u = Real.One / Real.Cosh(x.Value);
@@ -335,12 +404,14 @@ public Variable Tanh(Variable x)
 
     // Logarithmic functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Ln(T)"/>
     public Variable Ln(Variable x)
     {
         _nodes.Add(new(Real.One / x.Value, x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Ln(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log(T, T)"/>
     public Variable Log(Variable x, Variable b)
     {
         var lnB = Real.Ln(b.Value);
@@ -348,12 +419,14 @@ public Variable Log(Variable x, Variable b)
         return new(_nodes.Count - 1, Real.Log(x.Value, b.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log2(T)"/>
     public Variable Log2(Variable x)
     {
         _nodes.Add(new(Real.One / (Real.Ln2 * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Log2(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Log10(T)"/>
     public Variable Log10(Variable x)
     {
         _nodes.Add(new(Real.One / (Real.Ln10 * x.Value), x._index, _nodes.Count));
@@ -362,6 +435,7 @@ public Variable Log10(Variable x)
 
     // Power functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Pow(T, T)"/>
     public Variable Pow(Variable x, Variable y)
     {
         var pow = Real.Pow(x.Value, y.Value);
@@ -371,6 +445,7 @@ public Variable Pow(Variable x, Variable y)
 
     // Root functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cbrt(T)"/>
     public Variable Cbrt(Variable x)
     {
         var cbrt = Real.Cbrt(x.Value);
@@ -378,6 +453,7 @@ public Variable Cbrt(Variable x)
         return new(_nodes.Count - 1, cbrt);
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Root(T, T)"/>
     public Variable Root(Variable x, Variable n)
     {
         var root = Real.Root(x.Value, n.Value);
@@ -385,6 +461,7 @@ public Variable Root(Variable x, Variable n)
         return new(_nodes.Count - 1, root);
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sqrt(T)"/>
     public Variable Sqrt(Variable x)
     {
         var sqrt = Real.Sqrt(x.Value);
@@ -394,18 +471,21 @@ public Variable Sqrt(Variable x)
 
     // Trigonometric functions
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Acos(T)"/>
     public Variable Acos(Variable x)
     {
         _nodes.Add(new(-Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Acos(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Asin(T)"/>
     public Variable Asin(Variable x)
     {
         _nodes.Add(new(Real.One / Real.Sqrt(Real.One - x.Value * x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Asin(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Atan(T)"/>
     public Variable Atan(Variable x)
     {
         _nodes.Add(new(Real.One / (Real.One + x.Value * x.Value), x._index, _nodes.Count));
@@ -419,18 +499,21 @@ public Variable Atan2(Variable y, Variable x)
         return new(_nodes.Count - 1, Real.Atan2(y.Value, x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Cos(T)"/>
     public Variable Cos(Variable x)
     {
         _nodes.Add(new(-Real.Sin(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Cos(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Sin(T)"/>
     public Variable Sin(Variable x)
     {
         _nodes.Add(new(Real.Cos(x.Value), x._index, _nodes.Count));
         return new(_nodes.Count - 1, Real.Sin(x.Value));
     }
 
+    /// <inheritdoc cref="IDifferentiableFunctions{T}.Tan(T)"/>
     public Variable Tan(Variable x)
     {
         var sec = Real.One / Real.Cos(x.Value);

From e223e2f1ba533b882bc0e248bc4119e271a6bbad Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 15:05:51 -0600
Subject: [PATCH 08/18] Add ToReal method

---
 src/Mathematics.NET/Core/IReal.cs    | 6 ++++++
 src/Mathematics.NET/Core/Rational.cs | 5 ++++-
 src/Mathematics.NET/Core/Real.cs     | 2 ++
 3 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/Mathematics.NET/Core/IReal.cs b/src/Mathematics.NET/Core/IReal.cs
index 3b9a1bbb..67299ca3 100644
--- a/src/Mathematics.NET/Core/IReal.cs
+++ b/src/Mathematics.NET/Core/IReal.cs
@@ -91,4 +91,10 @@ public interface IReal<T>
     /// </returns>
     /// <exception cref="ArithmeticException">An overflow or underflow has occurred</exception>
     static abstract int Sign(T x);
+
+    /// <summary>Create a real number from a type that implements <typeparamref name="T"/></summary>
+    /// <param name="x">A type that implements <typeparamref name="T"/></param>
+    /// <returns>A real value</returns>
+    /// <exception cref="OverflowException">Thrown when the value cannot be converted to the type <see cref="Real"/></exception>
+    static abstract Real ToReal(T x);
 }
diff --git a/src/Mathematics.NET/Core/Rational.cs b/src/Mathematics.NET/Core/Rational.cs
index 0f7bfa6a..0a32bbd7 100644
--- a/src/Mathematics.NET/Core/Rational.cs
+++ b/src/Mathematics.NET/Core/Rational.cs
@@ -612,6 +612,9 @@ public static int Sign(Rational<T> x)
         return x._numerator > T.Zero ? 1 : -1;
     }
 
+    public static Real ToReal(Rational<T> x)
+        => checked(double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator));
+
     public static bool TryConvertFromChecked<U>(U value, out Rational<T> result)
         where U : INumberBase<U>
     {
@@ -715,7 +718,7 @@ public static explicit operator Rational<T>(double x)
 
     public static explicit operator checked Real(Rational<T> x) => checked(double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator));
 
-    public static explicit operator Real(Rational<T> x) => double.CreateSaturating(x._numerator) / double.CreateSaturating(x._denominator);
+    public static explicit operator Real(Rational<T> x) => ToReal(x);
 
     public static explicit operator checked double(Rational<T> x) => double.CreateChecked(x._numerator) / double.CreateChecked(x._denominator);
 
diff --git a/src/Mathematics.NET/Core/Real.cs b/src/Mathematics.NET/Core/Real.cs
index 8865dfbe..ba462f14 100644
--- a/src/Mathematics.NET/Core/Real.cs
+++ b/src/Mathematics.NET/Core/Real.cs
@@ -300,6 +300,8 @@ public static Real Reciprocate(Real x)
 
     public static int Sign(Real x) => Math.Sign(x._value);
 
+    public static Real ToReal(Real x) => x;
+
     public static bool TryConvertFromChecked<U>(U value, out Real result)
         where U : INumberBase<U>
     {

From 8f8bf28f1d759ae2aae51a46bda8f6dc250242f3 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 15:12:27 -0600
Subject: [PATCH 09/18] Add Atan2

---
 src/Mathematics.NET/Core/IReal.cs    | 6 ++++++
 src/Mathematics.NET/Core/Rational.cs | 2 ++
 2 files changed, 8 insertions(+)

diff --git a/src/Mathematics.NET/Core/IReal.cs b/src/Mathematics.NET/Core/IReal.cs
index 67299ca3..c663f184 100644
--- a/src/Mathematics.NET/Core/IReal.cs
+++ b/src/Mathematics.NET/Core/IReal.cs
@@ -47,6 +47,12 @@ public interface IReal<T>
     /// <summary>The backing value of the type</summary>
     double Value { get; }
 
+    /// <summary>Compute a quadrant-aware arctangent given two values</summary>
+    /// <param name="y">The first value</param>
+    /// <param name="x">The second value</param>
+    /// <returns>An angle</returns>
+    static abstract Real Atan2(T y, T x);
+
     /// <summary>Compute the ceiling function of a value</summary>
     /// <param name="x">A value</param>
     /// <returns>The smallest integer greater than or equal to the value</returns>
diff --git a/src/Mathematics.NET/Core/Rational.cs b/src/Mathematics.NET/Core/Rational.cs
index 0a32bbd7..210955e0 100644
--- a/src/Mathematics.NET/Core/Rational.cs
+++ b/src/Mathematics.NET/Core/Rational.cs
@@ -501,6 +501,8 @@ public static bool TryParse(ReadOnlySpan<char> s, NumberStyles style, IFormatPro
 
     public static Rational<T> Abs(Rational<T> x) => new(T.Abs(x._numerator), T.Abs(x._denominator));
 
+    public static Real Atan2(Rational<T> y, Rational<T> x) => Real.Atan2(ToReal(y), ToReal(x));
+
     public static Rational<T> Ceiling(Rational<T> x)
     {
         var (quotient, remainder) = T.DivRem(x._numerator, x._denominator);

From 1966e5370be34cdbae47cfcfea5328f9708297cc Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 15:12:36 -0600
Subject: [PATCH 10/18] Update GradientTape.cs

- Add documentation comment
---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 62b3a3b3..5eced7a6 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -492,6 +492,7 @@ public Variable Atan(Variable x)
         return new(_nodes.Count - 1, Real.Atan(x.Value));
     }
 
+    /// <inheritdoc cref="IReal{T}.Atan2(T, T)"/>
     public Variable Atan2(Variable y, Variable x)
     {
         var u = Real.One / (x.Value * x.Value + y.Value * y.Value);

From 8e88ad6641aef382d5d8ad051a5cb33ab54dbaec Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 17:14:06 -0600
Subject: [PATCH 11/18] Update Vector3.cs

- Add cross product
---
 src/Mathematics.NET/LinearAlgebra/Vector3.cs | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/Mathematics.NET/LinearAlgebra/Vector3.cs b/src/Mathematics.NET/LinearAlgebra/Vector3.cs
index ec103a13..5822a026 100644
--- a/src/Mathematics.NET/LinearAlgebra/Vector3.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Vector3.cs
@@ -178,6 +178,18 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Compute the cross product of two vectors.</summary>
+    /// <param name="left">The first vector</param>
+    /// <param name="right">The second vector</param>
+    /// <returns>The cross product</returns>
+    public static Vector3<T> Cross(Vector3<T> left, Vector3<T> right)
+    {
+        return new(
+            left.X2 * right.X3 - left.X3 * right.X2,
+            left.X3 * right.X1 - left.X1 * right.X3,
+            left.X1 * right.X2 - left.X2 * right.X1);
+    }
+
     public static T InnerProduct(Vector3<T> left, Vector3<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2 + T.Conjugate(left.X3) * right.X3;
 

From 39cd675d2a0b7614ddda53578b93a9b9bac78099 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 17:20:36 -0600
Subject: [PATCH 12/18] Return readonly values

- Add documentation comments
---
 src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs | 10 +++++++---
 src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs | 11 ++++++++---
 src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs | 12 +++++++++---
 src/Mathematics.NET/LinearAlgebra/Vector2.cs   |  4 ++--
 src/Mathematics.NET/LinearAlgebra/Vector3.cs   |  4 ++--
 src/Mathematics.NET/LinearAlgebra/Vector4.cs   |  4 ++--
 6 files changed, 30 insertions(+), 15 deletions(-)

diff --git a/src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs b/src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs
index 69904fdf..c3791f5f 100644
--- a/src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Matrix2x2.cs
@@ -194,12 +194,16 @@ public string ToString(string? format, IFormatProvider? provider)
         return string.Format(provider, builder.ToString());
     }
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix2x2<T> CreateDiagonal(T e11, T e22)
         => new(e11, T.Zero, T.Zero, e22);
 
     public readonly T Determinant() => E11 * E22 - E12 * E21;
 
-    public Matrix2x2<T> Inverse()
+    public readonly Matrix2x2<T> Inverse()
     {
         var det = Determinant();
         if (det == T.Zero)
@@ -214,7 +218,7 @@ public Matrix2x2<T> Inverse()
     public static bool IsNaM(Matrix2x2<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22);
 
-    public T Trace() => E11 + E22;
+    public readonly T Trace() => E11 + E22;
 
-    public Matrix2x2<T> Transpose() => new(E11, E21, E12, E22);
+    public readonly Matrix2x2<T> Transpose() => new(E11, E21, E12, E22);
 }
diff --git a/src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs b/src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs
index f82f4c38..e78df3cd 100644
--- a/src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Matrix3x3.cs
@@ -247,6 +247,11 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <param name="e33">The $ e_{33} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix3x3<T> CreateDiagonal(T e11, T e22, T e33)
     {
         return new(
@@ -268,7 +273,7 @@ public readonly T Determinant()
         return a * ei_fh - b * di_fg + c * dh_eg;
     }
 
-    public Matrix3x3<T> Inverse()
+    public readonly Matrix3x3<T> Inverse()
     {
         T a = E11, b = E12, c = E13;
         T d = E21, e = E22, f = E23;
@@ -313,9 +318,9 @@ public Matrix3x3<T> Inverse()
     public static bool IsNaM(Matrix3x3<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22) && T.IsNaN(matrix.E33);
 
-    public T Trace() => E11 + E22 + E33;
+    public readonly T Trace() => E11 + E22 + E33;
 
-    public Matrix3x3<T> Transpose()
+    public readonly Matrix3x3<T> Transpose()
     {
         return new(
             E11, E21, E31,
diff --git a/src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs b/src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs
index 7c51b2ce..e3a90f53 100644
--- a/src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Matrix4x4.cs
@@ -285,6 +285,12 @@ public string ToString(string? format, IFormatProvider? provider)
     // Methods
     //
 
+    /// <summary>Create a diagonal matrix from specified values along the diagonal.</summary>
+    /// <param name="e11">The $ e_{11} $ component</param>
+    /// <param name="e22">The $ e_{22} $ component</param>
+    /// <param name="e33">The $ e_{33} $ component</param>
+    /// <param name="e44">The $ e_{44} $ component</param>
+    /// <returns>A diagonal matrix</returns>
     public static Matrix4x4<T> CreateDiagonal(T e11, T e22, T e33, T e44)
     {
         return new(
@@ -315,7 +321,7 @@ public readonly T Determinant()
     }
 
     // TODO: Optimize
-    public Matrix4x4<T> Inverse()
+    public readonly Matrix4x4<T> Inverse()
     {
         T a = E11, b = E12, c = E13, d = E14;
         T e = E21, f = E22, g = E23, h = E24;
@@ -383,9 +389,9 @@ public Matrix4x4<T> Inverse()
     public static bool IsNaM(Matrix4x4<T> matrix)
         => T.IsNaN(matrix.E11) && T.IsNaN(matrix.E22) && T.IsNaN(matrix.E33) && T.IsNaN(matrix.E44);
 
-    public T Trace() => E11 + E22 + E33 + E44;
+    public readonly T Trace() => E11 + E22 + E33 + E44;
 
-    public Matrix4x4<T> Transpose()
+    public readonly Matrix4x4<T> Transpose()
     {
         return new(
             E11, E21, E31, E41,
diff --git a/src/Mathematics.NET/LinearAlgebra/Vector2.cs b/src/Mathematics.NET/LinearAlgebra/Vector2.cs
index 3bc2b5e3..d3f04928 100644
--- a/src/Mathematics.NET/LinearAlgebra/Vector2.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Vector2.cs
@@ -153,7 +153,7 @@ public string ToString(string? format, IFormatProvider? provider)
     public static T InnerProduct(Vector2<T> left, Vector2<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         var x0 = (X1 * T.Conjugate(X1)).Re;
         var x1 = (X2 * T.Conjugate(X2)).Re;
@@ -164,7 +164,7 @@ public Real Norm()
         return max * Real.Sqrt(x0 / maxSquared + x1 / maxSquared);
     }
 
-    public Vector2<T> Normalize()
+    public readonly Vector2<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm);
diff --git a/src/Mathematics.NET/LinearAlgebra/Vector3.cs b/src/Mathematics.NET/LinearAlgebra/Vector3.cs
index 5822a026..07900288 100644
--- a/src/Mathematics.NET/LinearAlgebra/Vector3.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Vector3.cs
@@ -193,7 +193,7 @@ public static Vector3<T> Cross(Vector3<T> left, Vector3<T> right)
     public static T InnerProduct(Vector3<T> left, Vector3<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2 + T.Conjugate(left.X3) * right.X3;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         Span<Real> components = stackalloc Real[3];
 
@@ -219,7 +219,7 @@ public Real Norm()
         return max * Real.Sqrt(partialSum);
     }
 
-    public Vector3<T> Normalize()
+    public readonly Vector3<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm, X3 / norm);
diff --git a/src/Mathematics.NET/LinearAlgebra/Vector4.cs b/src/Mathematics.NET/LinearAlgebra/Vector4.cs
index 304ea0a4..aa870c99 100644
--- a/src/Mathematics.NET/LinearAlgebra/Vector4.cs
+++ b/src/Mathematics.NET/LinearAlgebra/Vector4.cs
@@ -192,7 +192,7 @@ public string ToString(string? format, IFormatProvider? provider)
     public static T InnerProduct(Vector4<T> left, Vector4<T> right)
         => T.Conjugate(left.X1) * right.X1 + T.Conjugate(left.X2) * right.X2 + T.Conjugate(left.X3) * right.X3 + T.Conjugate(left.X4) * right.X4;
 
-    public Real Norm()
+    public readonly Real Norm()
     {
         Span<Real> components = stackalloc Real[4];
 
@@ -220,7 +220,7 @@ public Real Norm()
         return max * Real.Sqrt(partialSum);
     }
 
-    public Vector4<T> Normalize()
+    public readonly Vector4<T> Normalize()
     {
         var norm = Norm();
         return new(X1 / norm, X2 / norm, X3 / norm, X4 / norm);

From 317e8c94553fb9d5e21d214e994651621e5ae3d4 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 17:46:37 -0600
Subject: [PATCH 13/18] Update Mathematics.NET.Tests.csproj

---
 tests/Mathematics.NET.Tests/Mathematics.NET.Tests.csproj | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/Mathematics.NET.Tests/Mathematics.NET.Tests.csproj b/tests/Mathematics.NET.Tests/Mathematics.NET.Tests.csproj
index 5829ef8a..87e2a518 100644
--- a/tests/Mathematics.NET.Tests/Mathematics.NET.Tests.csproj
+++ b/tests/Mathematics.NET.Tests/Mathematics.NET.Tests.csproj
@@ -11,7 +11,7 @@
   </PropertyGroup>
 
   <ItemGroup>
-    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.7.2" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.8.0" />
     <PackageReference Include="MSTest.TestAdapter" Version="3.1.1" />
     <PackageReference Include="MSTest.TestFramework" Version="3.1.1" />
     <PackageReference Include="coverlet.collector" Version="6.0.0">

From e7f3763e1316fc3537b4329ab4c731dcceb84271 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 18:04:29 -0600
Subject: [PATCH 14/18] Update GradientTape.cs

- Rename some items
---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 5eced7a6..9e6d44a3 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -155,19 +155,19 @@ public void ReverseAccumulation(out ReadOnlySpan<Real> gradients, double seed =
         ref var start = ref MemoryMarshal.GetReference(nodes);
 
         var length = nodes.Length;
-        Span<Real> partialGradients = new Real[length];
-        partialGradients[length - 1] = seed;
+        Span<Real> gradientSpan = new Real[length];
+        gradientSpan[length - 1] = seed;
 
         for (int i = length - 1; i >= _variableCount; i--)
         {
             var node = Unsafe.Add(ref start, i);
-            var partialGradient = partialGradients[i];
+            var gradient = gradientSpan[i];
 
-            partialGradients[node.PX] += partialGradient * node.DX;
-            partialGradients[node.PY] += partialGradient * node.DY;
+            gradientSpan[node.PX] += gradient * node.DX;
+            gradientSpan[node.PY] += gradient * node.DY;
         }
 
-        gradients = partialGradients[.._variableCount];
+        gradients = gradientSpan[.._variableCount];
     }
 
     //

From 26b81128411885e4b632ba05e4f3ffde147a1e4c Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 23:27:31 -0600
Subject: [PATCH 15/18] Update Assert.cs

- Use ReadOnlySpan2D
- Add method for asserting ReadOnlySpan equality
- Rename parameters
---
 tests/Mathematics.NET.Tests/Assert.cs | 52 +++++++++++++++++++++++----
 1 file changed, 46 insertions(+), 6 deletions(-)

diff --git a/tests/Mathematics.NET.Tests/Assert.cs b/tests/Mathematics.NET.Tests/Assert.cs
index e81483eb..c27a9b0d 100644
--- a/tests/Mathematics.NET.Tests/Assert.cs
+++ b/tests/Mathematics.NET.Tests/Assert.cs
@@ -29,20 +29,26 @@
 
 namespace Mathematics.NET.Tests;
 
+/// <summary>Assert helpers for Mathemtics.NET</summary>
+/// <typeparam name="T">A type that implements <see cref="IComplex{T}"/></typeparam>
 public sealed class Assert<T>
     where T : IComplex<T>
 {
-    public static void AreApproximatelyEqual(T expected, T actual, Real delta)
+    /// <summary>Assert that two values are approximately equal.</summary>
+    /// <param name="expected">The expected value</param>
+    /// <param name="actual">The actual value</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void AreApproximatelyEqual(T expected, T actual, Real epsilon)
     {
         if (T.IsNaN(expected) && T.IsNaN(actual) || T.IsInfinity(expected) && T.IsInfinity(actual))
         {
             return;
         }
 
-        if (!Precision.AreApproximatelyEqual(expected, actual, delta))
+        if (!Precision.AreApproximatelyEqual(expected, actual, epsilon))
         {
             Assert.Fail($$"""
-                Actual value does not fall within the specifed margin of error, {{delta}}:
+                Actual value does not fall within the specifed margin of error, {{epsilon}}:
 
                 Expected: {{expected}}
                 Actual: {{actual}}
@@ -50,7 +56,41 @@ public static void AreApproximatelyEqual(T expected, T actual, Real delta)
         }
     }
 
-    public static void ElementsAreApproximatelyEqual(Span2D<T> expected, Span2D<T> actual, Real delta)
+    /// <summary>Assert that the elements in two read-only spans are approximately equal.</summary>
+    /// <param name="expected">A read-only span of expected values</param>
+    /// <param name="actual">A read-only span of actual values</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void ElementsAreApproximatelyEqual(ReadOnlySpan<T> expected, ReadOnlySpan<T> actual, Real epsilon)
+    {
+        if (expected.Length != actual.Length)
+        {
+            Assert.Fail($"Length of the actual array, {actual.Length}, does not match the length of the expected array, {expected.Length}");
+        }
+
+        for (int i = 0; i < expected.Length; i++)
+        {
+            if (T.IsNaN(expected[i]) && T.IsNaN(actual[i]) || T.IsInfinity(expected[i]) && T.IsInfinity(actual[i]))
+            {
+                continue;
+            }
+
+            if (!Precision.AreApproximatelyEqual(expected[i], actual[i], epsilon))
+            {
+                Assert.Fail($$"""
+                    Actual value at index {{i}} does not fall within the specified margin of error, {{epsilon}}
+
+                    Expected: {{expected[i]}}
+                    Actual: {{actual[i]}}
+                    """);
+            }
+        }
+    }
+
+    /// <summary>Assert that the elements in two 2D, read-only spans are approximately equal.</summary>
+    /// <param name="expected">A 2D, read-only span of expected values</param>
+    /// <param name="actual">A 2D, read-only span of actual values</param>
+    /// <param name="epsilon">A margin of error</param>
+    public static void ElementsAreApproximatelyEqual(ReadOnlySpan2D<T> expected, ReadOnlySpan2D<T> actual, Real epsilon)
     {
         if (expected.Height != actual.Height || expected.Width != actual.Width)
         {
@@ -66,10 +106,10 @@ public static void ElementsAreApproximatelyEqual(Span2D<T> expected, Span2D<T> a
                     continue;
                 }
 
-                if (!Precision.AreApproximatelyEqual(expected[i, j], actual[i, j], delta))
+                if (!Precision.AreApproximatelyEqual(expected[i, j], actual[i, j], epsilon))
                 {
                     Assert.Fail($$"""
-                        Actual value at row {{i}} and column {{j}} does not fall within the specifed margin of error, {{delta}}:
+                        Actual value at row {{i}} and column {{j}} does not fall within the specifed margin of error, {{epsilon}}:
 
                         Expected: {{expected[i, j]}}
                         Actual: {{actual[i, j]}}

From 6dde2562da0a079e789293d40f5d9974925add74 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Thu, 9 Nov 2023 23:51:30 -0600
Subject: [PATCH 16/18] Create GradientTapeTests.cs

---
 .../AutoDiff/GradientTapeTests.cs             | 561 ++++++++++++++++++
 1 file changed, 561 insertions(+)
 create mode 100644 tests/Mathematics.NET.Tests/AutoDiff/GradientTapeTests.cs

diff --git a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeTests.cs
new file mode 100644
index 00000000..477b8fdd
--- /dev/null
+++ b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeTests.cs
@@ -0,0 +1,561 @@
+// <copyright file="GradientTapeTests.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.Data;
+using Mathematics.NET.AutoDiff;
+
+namespace Mathematics.NET.Tests.AutoDiff;
+
+[TestClass]
+[TestCategory("AutoDiff"), TestCategory("Gradient Tape")]
+public sealed class GradientTapeTests
+{
+    [TestMethod]
+    [DataRow(0.123, -1.007651429146436)]
+    public void Acos_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Acos, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 1.396315794095838)]
+    public void Acosh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Acosh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1, 1)]
+    public void Add_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Add, left, right);
+
+        Assert<Real>.ElementsAreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1)]
+    public void Add_RealAndVariable_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(right);
+        _ = tape.Add(left, x);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1)]
+    public void Add_VariableAndReal_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(left);
+        _ = tape.Add(x, right);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(0.123, 1.007651429146436)]
+    public void Asin_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Asin, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.6308300845448597)]
+    public void Asinh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Asinh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.3979465955668749)]
+    public void Atan_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Atan, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0.334835801674179, -0.1760034342133505)]
+    public void Atan2_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Atan2, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, -1.94969779684149)]
+    public void Atanh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Atanh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.2903634877210767)]
+    public void Cbrt_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Cbrt, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, -0.942488801931697)]
+    public void Cos_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Cos, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 1.564468479304407)]
+    public void Cosh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Cosh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(2)]
+    [DataRow(3)]
+    [DataRow(4)]
+    public void CreateVariable_Multiple_TracksCorrectNumberOfVariables(int amount)
+    {
+        GradientTape tape = new();
+        for (int i = 0; i < amount; i++)
+        {
+            _ = tape.CreateVariable(0);
+        }
+
+        var actual = tape.VariableCount;
+
+        Assert.AreEqual(amount, actual);
+    }
+
+    [TestMethod]
+    [DataRow(1.23)]
+    public void CreateVariable_WithSeedValue_ReturnsVariable(double value)
+    {
+        GradientTape tape = new();
+
+        var actual = tape.CreateVariable(value).Value;
+
+        Assert.AreEqual(value, actual);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0.4273504273504274, -0.2246329169406093)]
+    public void Divide_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Divide, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, -0.2246329169406093)]
+    public void Divide_RealAndVariable_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(right);
+        _ = tape.Divide(left, x);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0.4273504273504274)]
+    public void Divide_VariableAndReal_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(left);
+        _ = tape.Divide(x, right);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 3.421229536289673)]
+    public void Exp_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Exp, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 1.625894476644487)]
+    public void Exp2_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Exp2, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 39.10350518430174)]
+    public void Exp10_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Exp10, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.813008130081301)]
+    public void Ln_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Ln, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0.9563103467806, -0.1224030239537303)]
+    public void Log_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Log, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 1.172922797470702)]
+    public void Log2_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Log2, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.35308494463679)]
+    public void Log10_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Log10, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1, 0)]
+    public void Modulo_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Modulo, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0)]
+    public void Modulo_RealAndVariable_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(right);
+        _ = tape.Modulo(left, x);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1)]
+    public void Modulo_VariableAndReal_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(left);
+        _ = tape.Modulo(x, right);
+        tape.ReverseAccumulation(out var gradient);
+
+        var actual = gradient[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 2.34, 1.23)]
+    public void Multiply_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Multiply, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-16);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1.23)]
+    public void Multiply_RealAndVariable_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(right);
+        _ = tape.Multiply(left, x);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 2.34)]
+    public void Multiply_VariableAndReal_ReturnsGradient(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(left);
+        _ = tape.Multiply(x, right);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 3.088081166620949, 0.3360299854573856)]
+    public void Pow_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Pow, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 0.3795771135606888, -0.04130373687338086)]
+    public void Root_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Root, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-15);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.3342377271245026)]
+    public void Sin_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Sin, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 1.856761056985266)]
+    public void Sinh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Sinh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.4508348173337161)]
+    public void Sqrt_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Sqrt, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1, -1)]
+    public void Subtract_TwoVariables_ReturnsGradients(double left, double right, double expectedLeft, double expectedRight)
+    {
+        GradientTape tape = new();
+        var expected = new Real[2] { expectedLeft, expectedRight };
+
+        var actual = ComputeGradients(tape, tape.Subtract, left, right);
+
+        Assert<Real>.AreApproximatelyEqual(expectedLeft, actual[0], 1e-16);
+        Assert<Real>.AreApproximatelyEqual(expectedRight, actual[1], 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, -1)]
+    public void Subtract_RealAndVariable_ReturnsGradients(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(right);
+        _ = tape.Subtract(left, x);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 2.34, 1)]
+    public void Subtract_VariableAndReal_ReturnsGradients(double left, double right, double expected)
+    {
+        GradientTape tape = new();
+        var x = tape.CreateVariable(left);
+        _ = tape.Subtract(x, right);
+        tape.ReverseAccumulation(out var gradients);
+
+        var actual = gradients[0];
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-16);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 8.95136077522624)]
+    public void Tan_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Tan, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    [TestMethod]
+    [DataRow(1.23, 0.2900600799721436)]
+    public void Tanh_Variable_ReturnsGradient(double input, double expected)
+    {
+        GradientTape tape = new();
+
+        var actual = ComputeGradient(tape, tape.Tanh, input);
+
+        Assert<Real>.AreApproximatelyEqual(expected, actual, 1e-15);
+    }
+
+    //
+    // Helpers
+    //
+
+    private static Real ComputeGradient(GradientTape tape, Func<Variable, Variable> function, Real input)
+    {
+        var x = tape.CreateVariable(input);
+        _ = function(x);
+        tape.ReverseAccumulation(out var gradients);
+        return gradients[0];
+    }
+
+    private static Real[] ComputeGradients(GradientTape tape, Func<Variable, Variable, Variable> function, Real left, Real right)
+    {
+        var x = tape.CreateVariable(left);
+        var y = tape.CreateVariable(right);
+        _ = function(x, y);
+        tape.ReverseAccumulation(out var gradients);
+        return gradients.ToArray();
+    }
+}

From 6373583a855c8509d3b7eb41321d1ef1c559fb4d Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Fri, 10 Nov 2023 01:10:04 -0600
Subject: [PATCH 17/18] Fix Atan2

---
 src/Mathematics.NET/AutoDiff/GradientTape.cs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 9e6d44a3..f8f94b48 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -496,7 +496,7 @@ public Variable Atan(Variable x)
     public Variable Atan2(Variable y, Variable x)
     {
         var u = Real.One / (x.Value * x.Value + y.Value * y.Value);
-        _nodes.Add(new(-x.Value * u, y.Value * u, y._index, x._index));
+        _nodes.Add(new(x.Value * u, -y.Value * u, y._index, x._index));
         return new(_nodes.Count - 1, Real.Atan2(y.Value, x.Value));
     }
 

From 71ae241e8bd8609bc138e748d74cfb1bf895fdf5 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Fri, 10 Nov 2023 01:13:58 -0600
Subject: [PATCH 18/18] Update Mathematics.NET.csproj

- Update version
---
 src/Mathematics.NET/Mathematics.NET.csproj | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Mathematics.NET/Mathematics.NET.csproj b/src/Mathematics.NET/Mathematics.NET.csproj
index a4f3e7ad..82d918dc 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.3</Version>
+    <Version>0.1.0-alpha.4</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.</Description>