diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index 7dd655b6..6b474bef 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -35,6 +35,7 @@ export
     # extrapolation/extrapolation.jl
     # monotonic/monotonic.jl
     # scaling/scaling.jl
+    # hermite/cubic.jl
 
 using LinearAlgebra, SparseArrays
 using StaticArrays, WoodburyMatrices, Ratios, AxisAlgorithms, OffsetArrays
diff --git a/src/hermite/cubic.jl b/src/hermite/cubic.jl
index 127380c4..dc67ce74 100644
--- a/src/hermite/cubic.jl
+++ b/src/hermite/cubic.jl
@@ -1,3 +1,5 @@
+export CubicHermite
+
 """
     CubicHermite
 
@@ -148,8 +150,3 @@ function hessian(ch::CubicHermite, x::Float64)::Float64
     d2 = 2 * (3x - 2x1 - x2) / (h * h)
     c1 * (y1 - y2) + d1 * dydx1 + d2 * dydx2
 end
-
-
-export CubicHermite
-export gradient
-export hessian
diff --git a/test/cubic_hermite.jl b/test/cubic_hermite.jl
index 9921dcf7..e351e30b 100644
--- a/test/cubic_hermite.jl
+++ b/test/cubic_hermite.jl
@@ -25,25 +25,25 @@ function cubic_hermite_test()
     for i = 1:length(xs)-1
         x = Float64(i)
         @test ch(x) == 1.0
-        @test gradient(ch, x) == 0.0
-        @test hessian(ch, x) == 0.0
+        @test Interpolations.gradient(ch, x) == 0.0
+        @test Interpolations.hessian(ch, x) == 0.0
         x = x + 0.25
         @test ch(x) == 1.0
-        @test gradient(ch, x) == 0.0
-        @test hessian(ch, x) == 0.0
+        @test Interpolations.gradient(ch, x) == 0.0
+        @test Interpolations.hessian(ch, x) == 0.0
         x = x + 0.25
         @test ch(x) == 1.0
-        @test gradient(ch, x) == 0.0
-        @test hessian(ch, x) == 0.0
+        @test Interpolations.gradient(ch, x) == 0.0
+        @test Interpolations.hessian(ch, x) == 0.0
         x = x + 0.25
         @test ch(x) == 1.0
-        @test gradient(ch, x) == 0.0
-        @test hessian(ch, x) == 0.0
+        @test Interpolations.gradient(ch, x) == 0.0
+        @test Interpolations.hessian(ch, x) == 0.0
     end
     # Ensure that the ight endpoint query doesn't read past the end of the array:
     @test ch(5.0) == 1.0
-    @test gradient(ch, 5.0) == 0.0
-    @test hessian(ch, 5.0) == 0.0
+    @test Interpolations.gradient(ch, 5.0) == 0.0
+    @test Interpolations.hessian(ch, 5.0) == 0.0
 
     # Now linear functions:
     a = 7.2
@@ -56,24 +56,24 @@ function cubic_hermite_test()
     for i = 1:length(xs)-1
         x = Float64(i)
         @test ch(x) ≈ a * x + b
-        @test gradient(ch, x) ≈ a
-        @test abs(hessian(ch, x)) < 3e-14
+        @test Interpolations.gradient(ch, x) ≈ a
+        @test abs(Interpolations.hessian(ch, x)) < 3e-14
         x = x + 0.25
         @test ch(x) ≈ a * x + b
-        @test gradient(ch, x) ≈ a
-        @test abs(hessian(ch, x)) < 3e-14
+        @test Interpolations.gradient(ch, x) ≈ a
+        @test abs(Interpolations.hessian(ch, x)) < 3e-14
         x = x + 0.25
         @test ch(x) ≈ a * x + b
-        @test gradient(ch, x) ≈ a
-        @test abs(hessian(ch, x)) < 3e-14
+        @test Interpolations.gradient(ch, x) ≈ a
+        @test abs(Interpolations.hessian(ch, x)) < 3e-14
         x = x + 0.25
         @test ch(x) ≈ a * x + b
-        @test gradient(ch, x) ≈ a
-        @test abs(hessian(ch, x)) < 3e-14
+        @test Interpolations.gradient(ch, x) ≈ a
+        @test abs(Interpolations.hessian(ch, x)) < 3e-14
     end
     @test ch(last(xs)) ≈ a * last(xs) + b
-    @test gradient(ch, last(xs)) ≈ a
-    @test abs(hessian(ch, last(xs))) < 3e-14
+    @test Interpolations.gradient(ch, last(xs)) ≈ a
+    @test abs(Interpolations.hessian(ch, last(xs))) < 3e-14
 
     # Now the interpolation condition:
     xs = zeros(50)
@@ -92,7 +92,7 @@ function cubic_hermite_test()
 
     for i = 1:50
         @test ch(xs[i]) ≈ ys[i]
-        @test gradient(ch, xs[i]) ≈ dydxs[i]
+        @test Interpolations.gradient(ch, xs[i]) ≈ dydxs[i]
     end
 
     # Now quadratics:
@@ -107,13 +107,13 @@ function cubic_hermite_test()
     for i = 1:200
         x = rand() * last(xs)
         @test ch(x) ≈ a * x * x + b * x + c
-        @test gradient(ch, x) ≈ 2 * a * x + b
-        @test hessian(ch, x) ≈ 2 * a
+        @test Interpolations.gradient(ch, x) ≈ 2 * a * x + b
+        @test Interpolations.hessian(ch, x) ≈ 2 * a
     end
     x = last(xs)
     @test ch(x) ≈ a * x * x + b * x + c
-    @test gradient(ch, x) ≈ 2 * x * a + b
-    @test hessian(ch, x) ≈ 2 * a
+    @test Interpolations.gradient(ch, x) ≈ 2 * x * a + b
+    @test Interpolations.hessian(ch, x) ≈ 2 * a
     # Cannot extrapolate:
     @test_throws DomainError ch(x + 0.1)
     @test_throws DomainError ch(xs[1] - 0.1)
diff --git a/test/gpu_support.jl b/test/gpu_support.jl
index 54d9d26c..811c5723 100644
--- a/test/gpu_support.jl
+++ b/test/gpu_support.jl
@@ -9,18 +9,18 @@ JLArrays.allowscalar(false)
     idx = 2.0:0.17:19.0
     jlidx = jl(collect(idx))
     @test itp.(idx) == collect(jlitp.(idx)) == collect(jlitp.(jlidx))
-    @test gradient.(Ref(itp), idx) ==
-        collect(gradient.(Ref(jlitp), idx)) ==
-        collect(gradient.(Ref(jlitp), jlidx))
+    @test Interpolations.gradient.(Ref(itp), idx) ==
+        collect(Interpolations.gradient.(Ref(jlitp), idx)) ==
+        collect(Interpolations.gradient.(Ref(jlitp), jlidx))
 
     sitp = scale(itp, A_x)
     jlsitp = jl(sitp)
     idx = 1.0:0.4:39.0
     jlidx = jl(collect(idx))
     @test sitp.(idx) == collect(jlsitp.(idx)) == collect(jlsitp.(jlidx))
-    @test gradient.(Ref(sitp), idx) ==
-        collect(gradient.(Ref(jlsitp), idx)) ==
-        collect(gradient.(Ref(jlsitp), jlidx))
+    @test Interpolations.gradient.(Ref(sitp), idx) ==
+        collect(Interpolations.gradient.(Ref(jlsitp), idx)) ==
+        collect(Interpolations.gradient.(Ref(jlsitp), jlidx))
 
 
     esitp = extrapolate(sitp, Flat())
@@ -28,18 +28,18 @@ JLArrays.allowscalar(false)
     idx = -1.0:0.84:41.0
     jlidx = jl(collect(idx))
     @test esitp.(idx) == collect(jlesitp.(idx)) == collect(jlesitp.(jlidx))
-    @test gradient.(Ref(esitp), idx) ==
-        collect(gradient.(Ref(jlesitp), idx)) ==
-        collect(gradient.(Ref(jlesitp), jlidx))
+    @test Interpolations.gradient.(Ref(esitp), idx) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), idx)) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), jlidx))
 
     esitp = extrapolate(sitp, 0.0)
     jlesitp = jl(esitp)
     idx = -1.0:0.84:41.0
     jlidx = jl(collect(idx))
     @test esitp.(idx) == collect(jlesitp.(idx)) == collect(jlesitp.(jlidx))
-    @test gradient.(Ref(esitp), idx) ==
-        collect(gradient.(Ref(jlesitp), idx)) ==
-        collect(gradient.(Ref(jlesitp), jlidx))
+    @test Interpolations.gradient.(Ref(esitp), idx) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), idx)) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), jlidx))
 end
 
 @testset "2d GPU Interpolation" begin
@@ -50,44 +50,44 @@ end
     idx = 2.0:0.17:19.0
     jlidx = jl(collect(idx))
     @test itp.(idx, idx') == collect(jlitp.(idx, idx')) == collect(jlitp.(jlidx, jlidx'))
-    @test gradient.(Ref(itp), idx, idx') ==
-        collect(gradient.(Ref(jlitp), idx, idx')) ==
-        collect(gradient.(Ref(jlitp), jlidx, jlidx'))
-    @test hessian.(Ref(itp), idx, idx') ==
-        collect(hessian.(Ref(jlitp), idx, idx')) ==
-        collect(hessian.(Ref(jlitp), jlidx, jlidx'))
+    @test Interpolations.gradient.(Ref(itp), idx, idx') ==
+        collect(Interpolations.gradient.(Ref(jlitp), idx, idx')) ==
+        collect(Interpolations.gradient.(Ref(jlitp), jlidx, jlidx'))
+    @test Interpolations.hessian.(Ref(itp), idx, idx') ==
+        collect(Interpolations.hessian.(Ref(jlitp), idx, idx')) ==
+        collect(Interpolations.hessian.(Ref(jlitp), jlidx, jlidx'))
 
     sitp = scale(itp, A_x, A_x)
     jlsitp = jl(sitp)
     idx = 1.0:0.4:39.0
     jlidx = jl(collect(idx))
     @test sitp.(idx, idx') == collect(jlsitp.(idx, idx')) == collect(jlsitp.(jlidx, jlidx'))
-    @test gradient.(Ref(sitp), idx, idx') ==
-        collect(gradient.(Ref(jlsitp), idx, idx')) ==
-        collect(gradient.(Ref(jlsitp), jlidx, jlidx'))
-    @test hessian.(Ref(sitp), idx, idx') ==
-        collect(hessian.(Ref(jlsitp), idx, idx')) ==
-        collect(hessian.(Ref(jlsitp), jlidx, jlidx'))
+    @test Interpolations.gradient.(Ref(sitp), idx, idx') ==
+        collect(Interpolations.gradient.(Ref(jlsitp), idx, idx')) ==
+        collect(Interpolations.gradient.(Ref(jlsitp), jlidx, jlidx'))
+    @test Interpolations.hessian.(Ref(sitp), idx, idx') ==
+        collect(Interpolations.hessian.(Ref(jlsitp), idx, idx')) ==
+        collect(Interpolations.hessian.(Ref(jlsitp), jlidx, jlidx'))
 
     esitp = extrapolate(sitp, Flat())
     jlesitp = jl(esitp)
     idx = -1.0:0.84:41.0
     jlidx = jl(collect(idx))
     @test esitp.(idx, idx') == collect(jlesitp.(idx, idx')) == collect(jlesitp.(jlidx, jlidx'))
-    # gradient for `extrapolation` is currently broken under CUDA
-    @test gradient.(Ref(esitp), idx, idx') ==
-        collect(gradient.(Ref(jlesitp), idx, idx')) ==
-        collect(gradient.(Ref(jlesitp), jlidx, jlidx'))
+    # Interpolations.gradient for `extrapolation` is currently broken under CUDA
+    @test Interpolations.gradient.(Ref(esitp), idx, idx') ==
+        collect(Interpolations.gradient.(Ref(jlesitp), idx, idx')) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), jlidx, jlidx'))
 
     esitp = extrapolate(sitp, 0.0)
     jlesitp = jl(esitp)
     idx = -1.0:0.84:41.0
     jlidx = jl(collect(idx))
     @test esitp.(idx, idx') == collect(jlesitp.(idx, idx')) == collect(jlesitp.(jlidx, jlidx'))
-    # gradient for `extrapolation` is currently broken under CUDA
-    @test gradient.(Ref(esitp), idx, idx') ==
-        collect(gradient.(Ref(jlesitp), idx, idx')) ==
-        collect(gradient.(Ref(jlesitp), jlidx, jlidx'))
+    # Interpolations.gradient for `extrapolation` is currently broken under CUDA
+    @test Interpolations.gradient.(Ref(esitp), idx, idx') ==
+        collect(Interpolations.gradient.(Ref(jlesitp), idx, idx')) ==
+        collect(Interpolations.gradient.(Ref(jlesitp), jlidx, jlidx'))
 end
 
 @testset "Lanczos on gpu" begin