diff --git a/Project.toml b/Project.toml
index e4b603d..4c1d4c3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,4 +1,4 @@
-name = "ArrayAD"
+name = "ArrayDiff"
 uuid = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
 authors = ["Benoît Legat <benoit.legat@gmail.com>"]
 version = "0.1.0"
diff --git a/README.md b/README.md
index d04a3ba..441231a 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# ArrayAD
+# ArrayDiff
 
 | **Build Status** |
 |:----------------:|
@@ -9,7 +9,7 @@ Experimental addition of array support to `MOI.Nonlinear.ReverseAD`
 > [!WARNING]
 > This code is still very much experimental
 
-[build-img]: https://github.com/blegat/ArrayAD.jl/actions/workflows/ci.yml/badge.svg?branch=main
-[build-url]: https://github.com/blegat/ArrayAD.jl/actions?query=workflow%3ACI
-[codecov-img]: https://codecov.io/gh/blegat/ArrayAD.jl/branch/main/graph/badge.svg
-[codecov-url]: https://codecov.io/gh/blegat/ArrayAD.jl/branch/main
+[build-img]: https://github.com/blegat/ArrayDiff.jl/actions/workflows/ci.yml/badge.svg?branch=main
+[build-url]: https://github.com/blegat/ArrayDiff.jl/actions?query=workflow%3ACI
+[codecov-img]: https://codecov.io/gh/blegat/ArrayDiff.jl/branch/main/graph/badge.svg
+[codecov-url]: https://codecov.io/gh/blegat/ArrayDiff.jl/branch/main
diff --git a/src/ArrayAD.jl b/src/ArrayDiff.jl
similarity index 98%
rename from src/ArrayAD.jl
rename to src/ArrayDiff.jl
index dabf51b..1224763 100644
--- a/src/ArrayAD.jl
+++ b/src/ArrayDiff.jl
@@ -4,7 +4,7 @@
 # Use of this source code is governed by an MIT-style license that can be found
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
-module ArrayAD
+module ArrayDiff
 
 import ForwardDiff
 import MathOptInterface as MOI
diff --git a/src/utils.jl b/src/utils.jl
index 9bea2fd..d1d69e4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -188,8 +188,8 @@ julia> x = [(1, 2, 3), (4, 5, 6), (7, 8, 9)]
  (4, 5, 6)
  (7, 8, 9)
 
-julia> MOI.Nonlinear.ArrayAD._reinterpret_unsafe(NTuple{2,Int}, x)
-4-element MathOptInterface.Nonlinear.ArrayAD._UnsafeVectorView{Tuple{Int64, Int64}}:
+julia> MOI.Nonlinear.ArrayDiff._reinterpret_unsafe(NTuple{2,Int}, x)
+4-element MathOptInterface.Nonlinear.ArrayDiff._UnsafeVectorView{Tuple{Int64, Int64}}:
  (1, 2)
  (3, 4)
  (5, 6)
diff --git a/test/Project.toml b/test/Project.toml
index 810d8d5..7ed96af 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,5 @@
 [deps]
-ArrayAD = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
+ArrayDiff = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
diff --git a/test/ReverseAD.jl b/test/ReverseAD.jl
index cf9e9a9..552ff3f 100644
--- a/test/ReverseAD.jl
+++ b/test/ReverseAD.jl
@@ -13,8 +13,8 @@ import SparseArrays
 import MathOptInterface as MOI
 const Nonlinear = MOI.Nonlinear
 
-import ArrayAD
-const Coloring = ArrayAD.Coloring
+import ArrayDiff
+const Coloring = ArrayDiff.Coloring
 
 function runtests()
     for name in names(@__MODULE__; all = true)
@@ -31,7 +31,7 @@ function test_objective_quadratic_univariate()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :($x^2 + 1))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2]) == 1.2^2 + 1
     g = [NaN]
@@ -62,7 +62,7 @@ function test_objective_and_constraints_quadratic_univariate()
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :($x^2 + 1))
     Nonlinear.add_constraint(model, :($x^2), MOI.LessThan(2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2]) == 1.2^2 + 1
     g = [NaN]
@@ -98,7 +98,7 @@ function test_objective_quadratic_multivariate()
     y = MOI.VariableIndex(2)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :($x^2 + $x * $y + $y^2))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2 + 1.2 * 2.3 + 2.3^2
     g = [NaN, NaN]
@@ -135,7 +135,7 @@ function test_objective_quadratic_multivariate_subexpressions()
     ey = Nonlinear.add_expression(model, :($y^2))
     exy = Nonlinear.add_expression(model, :($ex + $x * $y))
     Nonlinear.set_objective(model, :($exy + $ey))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     val = [1.2, 2.3]
     @test MOI.eval_objective(evaluator, val) == 1.2^2 + 1.2 * 2.3 + 2.3^2
@@ -177,7 +177,7 @@ function test_objective_ifelse_comparison()
     y = MOI.VariableIndex(2)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(ifelse(1 <= $x <= 2, $x^2, $y^2)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2
     @test MOI.eval_objective(evaluator, [2.2, 2.3]) == 2.3^2
@@ -194,7 +194,7 @@ function test_objective_ifelse_logic()
     y = MOI.VariableIndex(2)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(ifelse(1 <= $x && $x <= 2, $x^2, $y^2)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2
     @test MOI.eval_objective(evaluator, [2.2, 2.3]) == 2.3^2
@@ -211,7 +211,7 @@ function test_objective_parameter()
     model = Nonlinear.Model()
     p = Nonlinear.add_parameter(model, 1.2)
     Nonlinear.set_objective(model, :($p * $x))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.3]) == 1.2 * 1.3
     g = [NaN]
@@ -227,7 +227,7 @@ function test_objective_subexpression()
     expr = Nonlinear.add_expression(model, input)
     expr_2 = Nonlinear.add_expression(model, :($expr^2))
     Nonlinear.set_objective(model, :($expr_2^2))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad])
     @test MOI.eval_objective(evaluator, [1.3]) == ((1.3^2 + 1)^2)^2
     g = [NaN]
@@ -240,7 +240,7 @@ function test_constraint_quadratic_univariate()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.add_constraint(model, :($x^2), MOI.LessThan(2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     g = [NaN]
     x_val = [1.2]
@@ -266,7 +266,7 @@ function test_constraint_quadratic_multivariate()
     y = MOI.VariableIndex(2)
     model = Nonlinear.Model()
     Nonlinear.add_constraint(model, :($x^2 + $x * $y + $y^2), MOI.LessThan(2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     g = [NaN]
     x_val = [1.2, 2.3]
@@ -292,7 +292,7 @@ function test_constraint_quadratic_multivariate_subexpressions()
     ey = Nonlinear.add_expression(model, :($y^2))
     exy = Nonlinear.add_expression(model, :($ex + $x * $y))
     Nonlinear.add_constraint(model, :($exy + $ey), MOI.LessThan(2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     g = [NaN]
     x_val = [1.2, 2.3]
@@ -339,7 +339,7 @@ function test_hessian_sparsity_registered_function()
     model = Nonlinear.Model()
     Nonlinear.register_operator(model, :f, 2, f, ∇f, ∇²f)
     Nonlinear.set_objective(model, :(f($x, $z) + $y^2))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y, z])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y, z])
     @test :Hess in MOI.features_available(evaluator)
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.hessian_lagrangian_structure(evaluator) ==
@@ -369,7 +369,7 @@ function test_hessian_sparsity_registered_rosenbrock()
     model = Nonlinear.Model()
     Nonlinear.register_operator(model, :rosenbrock, 2, f, ∇f, ∇²f)
     Nonlinear.set_objective(model, :(rosenbrock($x, $y)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     @test :Hess in MOI.features_available(evaluator)
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.hessian_lagrangian_structure(evaluator) ==
@@ -399,7 +399,7 @@ function test_hessian_registered_error()
     model = Nonlinear.Model()
     Nonlinear.register_operator(model, :rosenbrock, 2, f, ∇f, ∇²f)
     Nonlinear.set_objective(model, :(rosenbrock($x, $y)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     H = fill(NaN, 3)
     @test_throws(
@@ -496,7 +496,7 @@ function test_derivatives()
     b = MOI.VariableIndex(2)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(sin($a^2) + cos($b * 4) / 5 - 2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [a, b])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [a, b])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     x = [2.0, 3.0]
     @test MOI.eval_objective(evaluator, x) ==
@@ -522,7 +522,7 @@ function test_NLPBlockData()
     Nonlinear.add_constraint(model, :($x - 2), MOI.GreaterThan(0.0))
     Nonlinear.add_constraint(model, :($x - 3), MOI.EqualTo(0.0))
     Nonlinear.add_constraint(model, :($x), MOI.Interval(4.0, 5.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     block = MOI.NLPBlockData(evaluator)
     @test block.constraint_bounds == [
         MOI.NLPBoundsPair(-Inf, 0.0),
@@ -544,11 +544,11 @@ function test_linearity()
         ex = Nonlinear.add_expression(model, input)
         expr = model[ex]
         adj = Nonlinear.adjacency_matrix(expr.nodes)
-        nodes = ArrayAD._replace_moi_variables(expr.nodes, variables)
-        ret = ArrayAD._classify_linearity(nodes, adj, ArrayAD.Linearity[])
+        nodes = ArrayDiff._replace_moi_variables(expr.nodes, variables)
+        ret = ArrayDiff._classify_linearity(nodes, adj, ArrayDiff.Linearity[])
         @test ret[1] == test_value
         indexed_set = Coloring.IndexedSet(100)
-        edge_list = ArrayAD._compute_hessian_sparsity(
+        edge_list = ArrayDiff._compute_hessian_sparsity(
             nodes,
             adj,
             ret,
@@ -556,14 +556,14 @@ function test_linearity()
             Set{Tuple{Int,Int}}[],
             Vector{Int}[],
         )
-        if ret[1] != ArrayAD.NONLINEAR
+        if ret[1] != ArrayDiff.NONLINEAR
             @test length(edge_list) == 0
         elseif length(IJ) > 0
             @test IJ == edge_list
         end
         if length(indices) > 0
             empty!(indexed_set)
-            ArrayAD._compute_gradient_sparsity!(indexed_set, nodes)
+            ArrayDiff._compute_gradient_sparsity!(indexed_set, nodes)
             ix = sort(collect(indexed_set))
             @test indices == ix
         end
@@ -571,58 +571,58 @@ function test_linearity()
     end
     _test_linearity(
         :(sin($x^2) + cos($y * 4) - 2.0),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(2, 2), (1, 1)]),
         [1, 2],
     )
-    _test_linearity(:(3 * 4 * ($x + $y)), ArrayAD.LINEAR)
+    _test_linearity(:(3 * 4 * ($x + $y)), ArrayDiff.LINEAR)
     _test_linearity(
         :($z * $y),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(3, 2), (3, 3), (2, 2)]),
         [2, 3],
     )
-    _test_linearity(:(3 + 4), ArrayAD.CONSTANT)
-    _test_linearity(:(sin(3) + $x), ArrayAD.LINEAR)
+    _test_linearity(:(3 + 4), ArrayDiff.CONSTANT)
+    _test_linearity(:(sin(3) + $x), ArrayDiff.LINEAR)
     _test_linearity(
         :(cos($z) * sin(3) + $x),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(3, 3)]),
         [1, 3],
     )
-    _test_linearity(:($x - 3 * $y), ArrayAD.LINEAR)
-    _test_linearity(:(-$x), ArrayAD.LINEAR)
-    _test_linearity(:(+$x), ArrayAD.LINEAR)
+    _test_linearity(:($x - 3 * $y), ArrayDiff.LINEAR)
+    _test_linearity(:(-$x), ArrayDiff.LINEAR)
+    _test_linearity(:(+$x), ArrayDiff.LINEAR)
     _test_linearity(
         :($x^$y),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(2, 2), (1, 1), (2, 1)]),
         [1, 2],
     )
-    _test_linearity(:($x / 3 + $y), ArrayAD.LINEAR)
+    _test_linearity(:($x / 3 + $y), ArrayDiff.LINEAR)
     _test_linearity(
         :(3 / ($x * $y)),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(2, 2), (1, 1), (2, 1)]),
         [1, 2],
     )
-    _test_linearity(:(1 / ($x + 3)), ArrayAD.NONLINEAR, Set([(1, 1)]), [1])
+    _test_linearity(:(1 / ($x + 3)), ArrayDiff.NONLINEAR, Set([(1, 1)]), [1])
     _test_linearity(
         :(ifelse($x <= 1, $x, $y)),
-        ArrayAD.PIECEWISE_LINEAR,
+        ArrayDiff.PIECEWISE_LINEAR,
         Set([]),
         [],
     )
     _test_linearity(
         :(ifelse($x <= 1, $x^2, $y)),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(1, 1)]),
         [1, 2],
     )
-    _test_linearity(:(ifelse(1 <= 1, 2, 3)), ArrayAD.CONSTANT)
+    _test_linearity(:(ifelse(1 <= 1, 2, 3)), ArrayDiff.CONSTANT)
     _test_linearity(
         :(1 / ifelse($x < 1, $x, 0)),
-        ArrayAD.NONLINEAR,
+        ArrayDiff.NONLINEAR,
         Set([(1, 1)]),
         [1],
     )
@@ -634,11 +634,12 @@ function test_linearity_no_hess()
     model = Nonlinear.Model()
     ex = Nonlinear.add_expression(model, :($x + 1))
     Nonlinear.set_objective(model, ex)
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac])
     # We initialized without the need for the hessian so
     # the linearity shouldn't be computed.
-    @test only(evaluator.backend.subexpressions).linearity == ArrayAD.NONLINEAR
+    @test only(evaluator.backend.subexpressions).linearity ==
+          ArrayDiff.NONLINEAR
     return
 end
 
@@ -648,8 +649,11 @@ function test_dual_forward()
     function _test_dual_forward(input, x_input, test_value)
         model = Nonlinear.Model()
         Nonlinear.set_objective(model, input)
-        evaluator =
-            Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y])
+        evaluator = Nonlinear.Evaluator(
+            model,
+            ArrayDiff.Mode(),
+            MOI.VariableIndex[x, y],
+        )
         MOI.initialize(evaluator, [:Grad])
         ∇f = fill(NaN, 2)
         MOI.eval_objective_gradient(evaluator, ∇f, x_input)
@@ -694,7 +698,7 @@ function test_gradient_registered_function()
     Nonlinear.register_operator(model, :c, 1, cos, x -> -sin(x), x -> -cos(x))
     Nonlinear.set_objective(model, :(Φ($y, $x - 1) * c($z)))
     evaluator =
-        Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y, z])
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x, y, z])
     MOI.initialize(evaluator, [:Grad])
     ∇f = fill(NaN, 3)
     x_input = [2.0, 3.0, 4.0]
@@ -715,7 +719,8 @@ function test_gradient_jump_855()
         model,
         :(ifelse($x <= 3.0, ($x - 2.0)^2, 2 * log($x - 2.0) + 1.0)),
     )
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     ∇f = fill(NaN, 1)
     MOI.eval_objective_gradient(evaluator, ∇f, [-1.0])
@@ -729,7 +734,8 @@ function test_gradient_abs()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(abs($x)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     ∇f = fill(NaN, 1)
     MOI.eval_objective_gradient(evaluator, ∇f, [2.0])
@@ -745,7 +751,7 @@ function test_gradient_trig()
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(sin($x^2) + cos($y * 4) / 5 - 2.0))
     evaluator =
-        Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y])
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x, y])
     MOI.initialize(evaluator, [:Grad])
     ∇f = fill(NaN, 2)
     MOI.eval_objective_gradient(evaluator, ∇f, [2.0, 3.0])
@@ -757,7 +763,8 @@ function test_gradient_logical()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :($x > 0.5 && $x < 0.9))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     @test MOI.eval_objective(evaluator, [1.5]) == 0.0
     ∇f = fill(NaN, 1)
@@ -770,7 +777,8 @@ function test_gradient_ifelse()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(ifelse($x >= 0.5 || $x < 0.1, $x, 5)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     @test MOI.eval_objective(evaluator, [1.5]) == 1.5
     ∇f = fill(NaN, 1)
@@ -789,7 +797,8 @@ function test_gradient_sqrt_nan()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(sqrt($x)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     @test isnan(MOI.eval_objective(evaluator, [-1.5]))
     ∇f = fill(Inf, 1)
@@ -805,7 +814,7 @@ function test_gradient_variable_power()
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :((1 / $x)^$y - $z))
     evaluator =
-        Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y, z])
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x, y, z])
     MOI.initialize(evaluator, [:Grad])
     x_input = [2.5, 3.5, 1.0]
     @test MOI.eval_objective(evaluator, x_input) == (1 / 2.5)^3.5 - 1.0
@@ -824,7 +833,8 @@ function test_single_parameter()
     model = Nonlinear.Model()
     p = Nonlinear.add_parameter(model, 105.2)
     Nonlinear.set_objective(model, :($p))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x])
+    evaluator =
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x])
     MOI.initialize(evaluator, [:Grad])
     @test MOI.eval_objective(evaluator, [-0.1]) == 105.2
     return
@@ -838,7 +848,7 @@ function test_gradient_nested_subexpressions()
     ex2 = Nonlinear.add_expression(model, :($ex1))
     Nonlinear.set_objective(model, ex2)
     evaluator =
-        Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y])
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x, y])
     MOI.initialize(evaluator, [:Grad])
     ∇f = fill(NaN, 2)
     MOI.eval_objective_gradient(evaluator, ∇f, [2.0, 3.0])
@@ -852,7 +862,7 @@ function test_gradient_view()
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(($x - 1)^2 + 4 * ($y - $x^2)^2))
     evaluator =
-        Nonlinear.Evaluator(model, ArrayAD.Mode(), MOI.VariableIndex[x, y])
+        Nonlinear.Evaluator(model, ArrayDiff.Mode(), MOI.VariableIndex[x, y])
     MOI.initialize(evaluator, [:Grad])
     X_input = [0.0; 1.0; 2.0; 3.0]
     for idx in [1:2, 3:4, [1, 3], 4:-2:2]
@@ -895,7 +905,7 @@ function _test_odd_chunks_Hessian_products(N)
     model = Nonlinear.Model()
     x = MOI.VariableIndex.(1:N)
     Nonlinear.set_objective(model, Expr(:call, :*, x...))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(evaluator, [:Hess])
     hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator)
     V = zeros(length(hessian_sparsity))
@@ -925,7 +935,7 @@ function test_jacobians_and_jacvec()
     Nonlinear.set_objective(model, :($a * $b + $c^2))
     Nonlinear.add_constraint(model, :($c * $b), MOI.LessThan(1.0))
     Nonlinear.add_constraint(model, :($a^2 / 2), MOI.LessThan(1.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(evaluator, [:Jac, :JacVec])
     values = [1.0, 2.0, 3.0] # For a, b, c.
     jacobian_sparsity = MOI.jacobian_structure(evaluator)
@@ -957,7 +967,7 @@ function test_jacobians_and_jacvec_with_subexpressions()
     Nonlinear.set_objective(model, :($a * $b + $c^2))
     Nonlinear.add_constraint(model, :($bc), MOI.LessThan(1.0))
     Nonlinear.add_constraint(model, :($a^2 / 2), MOI.LessThan(1.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(evaluator, [:Jac, :JacVec])
     values = [1.0, 2.0, 3.0] # For a, b, c.
     jacobian_sparsity = MOI.jacobian_structure(evaluator)
@@ -985,7 +995,7 @@ function test_pow_complex_result()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(ifelse($x > 0, $x^1.5, -(-$x)^1.5)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     x = [-2.0]
     @test MOI.eval_objective(evaluator, x) ≈ -(2^1.5)
@@ -1005,7 +1015,7 @@ function test_constraint_gradient()
     model = Nonlinear.Model()
     Nonlinear.add_constraint(model, :($x^2 + $x * $y + $y^2), MOI.LessThan(2.0))
     Nonlinear.add_constraint(model, :(cos($y)), MOI.LessThan(2.0))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac])
     @test MOI.constraint_gradient_structure(evaluator, 1) == [1, 2]
     @test MOI.constraint_gradient_structure(evaluator, 2) == [2]
@@ -1024,7 +1034,7 @@ function test_hessian_length()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(log($x)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Hess])
     H = Float64[]
     got, want = 0, 1
@@ -1042,7 +1052,7 @@ function test_jacobian_length()
     x = MOI.VariableIndex(1)
     model = Nonlinear.Model()
     Nonlinear.add_constraint(model, :(sin($x)), MOI.LessThan(0.5))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Jac])
     J = Float64[]
     @test_throws BoundsError MOI.eval_constraint_jacobian(evaluator, J, [1.0])
@@ -1054,7 +1064,7 @@ function test_timers()
     model = Nonlinear.Model()
     Nonlinear.set_objective(model, :(log($x)))
     Nonlinear.add_constraint(model, :(sin($x)), MOI.LessThan(0.5))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     y = [1.2]
     g = [NaN]
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
@@ -1094,7 +1104,7 @@ function test_varying_length_x()
     model = MOI.Nonlinear.Model()
     x = MOI.VariableIndex(1)
     MOI.Nonlinear.set_objective(model, :(sin($x)))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, Symbol[:Jac])
     ∇f = [NaN]
     MOI.eval_objective_gradient(evaluator, ∇f, [1.0, 2.0])
@@ -1110,7 +1120,7 @@ function test_univariate_operator_with_no_second_order()
     MOI.Nonlinear.register_operator(model, :op_f, 1, f, df)
     x = MOI.VariableIndex(1)
     MOI.Nonlinear.add_constraint(model, :(op_f($x)), MOI.LessThan(2.0))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     @test !(:Hess in MOI.features_available(evaluator))
     MOI.initialize(evaluator, [:Grad, :Jac])
     J = zeros(length(MOI.jacobian_structure(evaluator)))
@@ -1122,7 +1132,7 @@ end
 function test_no_objective()
     model = Nonlinear.Model()
     x = MOI.VariableIndex(1)
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad])
     @test_throws(
         ErrorException("No nonlinear objective."),
@@ -1140,7 +1150,7 @@ function test_x_power_1()
     model = Nonlinear.Model()
     x = MOI.VariableIndex(1)
     MOI.Nonlinear.set_objective(model, :($x^1))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Hess])
     @test MOI.eval_objective(evaluator, [2.0]) ≈ 2.0
     H = [NaN]
@@ -1154,7 +1164,7 @@ function test_variable_first_node_in_tape()
     x = MOI.VariableIndex(1)
     expr = MOI.Nonlinear.add_expression(model, :($x))
     MOI.Nonlinear.set_objective(model, :(sin($expr)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     H = [NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [2.0], 1.5, [])
@@ -1168,7 +1178,7 @@ function test_subexpression_first_node_in_tape()
     expr = MOI.Nonlinear.add_expression(model, :($x))
     expr2 = MOI.Nonlinear.add_expression(model, :($expr))
     MOI.Nonlinear.set_objective(model, :(sin($expr2)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     H = [NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [2.0], 1.5, [])
@@ -1181,7 +1191,7 @@ function test_parameter_in_hessian()
     x = MOI.VariableIndex(1)
     p = MOI.Nonlinear.add_parameter(model, 3.0)
     MOI.Nonlinear.set_objective(model, :(sin($x + $p)))
-    evaluator = Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     H = [NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [2.0], 1.5, [])
@@ -1192,7 +1202,7 @@ end
 function test_unsafe_vector_view()
     x = Float64[]
     GC.@preserve x begin
-        view = ArrayAD._UnsafeVectorView(x, 3)
+        view = ArrayDiff._UnsafeVectorView(x, 3)
         @test length(x) == 3
         view[2] = 1.0
         @test x[2] == 1.0
@@ -1205,7 +1215,7 @@ function test_classify_linearity_ifelse()
     y = MOI.VariableIndex(2)
     model = MOI.Nonlinear.Model()
     MOI.Nonlinear.set_objective(model, :(ifelse($y, $x, 1)))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.2, 1.0]) == 1.2
     @test MOI.eval_objective(evaluator, [1.2, 0.0]) == 1.0
@@ -1218,7 +1228,7 @@ function test_classify_linearity_logic()
     y = MOI.VariableIndex(2)
     model = MOI.Nonlinear.Model()
     MOI.Nonlinear.set_objective(model, :($x && $y))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(evaluator, [1.0, 1.0]) == 1.0
     @test MOI.eval_objective(evaluator, [0.0, 1.0]) == 0.0
@@ -1235,7 +1245,7 @@ function test_hessian_sparsity_with_subexpressions()
     expr = MOI.Nonlinear.add_expression(model, :($x * $y))
     expr2 = MOI.Nonlinear.add_expression(model, :($expr))
     MOI.Nonlinear.set_objective(model, :(sin($expr2)))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x, y])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     MOI.hessian_lagrangian_structure(evaluator)
     return
@@ -1250,7 +1260,7 @@ function test_toposort_subexpressions()
     d = MOI.Nonlinear.add_expression(model, :($c + $b))
     MOI.Nonlinear.add_constraint(model, :($d), MOI.LessThan(1.0))
     MOI.Nonlinear.add_constraint(model, :($c), MOI.LessThan(1.0))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), [x])
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     g = [NaN, NaN]
     MOI.eval_constraint(evaluator, g, [2.0])
@@ -1272,7 +1282,7 @@ function test_eval_user_defined_operator_ForwardDiff_gradient!()
         :($ex^3 + sin($(x[2])) / ψ($(x[2])) + t($(x[3]), $(x[4]))),
         MOI.LessThan(0.0),
     )
-    d = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    d = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(d, [:Jac])
     X = [1.1, 1.2, 1.3, 1.4]
     g = [NaN]
@@ -1304,7 +1314,7 @@ function test_eval_user_defined_operator_type_mismatch()
         :($ex^3 + sin($(x[2])) / ψ($(x[2])) + t($(x[3]), $(x[4]))),
         MOI.LessThan(0.0),
     )
-    d = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    d = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(d, [:Jac])
     X = [1.1, 1.2, 1.3, 1.4]
     g = [NaN]
@@ -1320,12 +1330,12 @@ end
 function test_generate_hessian_slice_inner()
     # Test that it evaluates without error. The code contents are tested
     # elsewhere.
-    ArrayAD._generate_hessian_slice_inner()
+    ArrayDiff._generate_hessian_slice_inner()
     d = ex = nothing  # These arguments are untyped and not needed for this test
-    for id in [0, ArrayAD.MAX_CHUNK + 1]
+    for id in [0, ArrayDiff.MAX_CHUNK + 1]
         @test_throws(
             ErrorException("Invalid chunk size: $id"),
-            ArrayAD._hessian_slice_inner(d, ex, id),
+            ArrayDiff._hessian_slice_inner(d, ex, id),
         )
     end
     return
@@ -1340,7 +1350,7 @@ function test_hessian_reinterpret_unsafe()
         MOI.EqualTo(0.0),
     )
     MOI.Nonlinear.add_constraint(model, :($(x[4]) * $(x[5])), MOI.EqualTo(1.0))
-    evaluator = MOI.Nonlinear.Evaluator(model, ArrayAD.Mode(), x)
+    evaluator = MOI.Nonlinear.Evaluator(model, ArrayDiff.Mode(), x)
     MOI.initialize(evaluator, [:Hess])
     H_s = MOI.hessian_lagrangian_structure(evaluator)
     H = zeros(length(H_s))