diff --git a/src/integral.jl b/src/integral.jl
index 1a1c7878..0fee71fb 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -100,7 +100,7 @@ function _integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, N), ones(FP, N); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, N), _ones(FP, N); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index fc5905b3..7b0f44ab 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -82,7 +82,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index 4d4c4c89..cb87717e 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -68,11 +68,11 @@ function integral(
 
     # Integrate f along the Line
     differential(line, x) = t′(x) * _units(line(0))
-    integrand(x::AbstractVector) = f(line(t(x[1]))) * differential(line, x[1])
+    integrand(xs) = f(line(t(xs[1]))) * differential(line, xs[1])
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = FP[0.5]
+    testpoint_parametriccoord = (FP(0.5),)
     integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index 1b34405e..b6de0e5d 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -91,7 +91,9 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, -ones(FP, 2), ones(FP, 2); rule.kwargs...)[1]
+    a = .-_ones(FP, 2)
+    b = _ones(FP, 2)
+    value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index 9008b4a8..84d6faf2 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -82,7 +82,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP, 1), ones(FP, 1); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, _zeros(FP, 1), _ones(FP, 1); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index c63e715b..ea0c7c11 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -87,7 +87,7 @@ function integral(
         v = R * (1 - b / (a + b))
         return f(triangle(u, v)) * R / (a + b)^2
     end
-    ∫ = HCubature.hcubature(integrand, zeros(FP, 2), FP[1, π / 2], rule.kwargs...)[1]
+    ∫ = HCubature.hcubature(integrand, _zeros(FP, 2), (FP(1), FP(π / 2)), rule.kwargs...)[1]
 
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
     return 2 * Meshes.area(triangle) .* ∫
diff --git a/src/utils.jl b/src/utils.jl
index e98ade74..e69a38ef 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -14,6 +14,14 @@ function _error_unsupported_combination(geometry, rule)
     throw(ArgumentError(msg))
 end
 
+# Return an NTuple{N, T} of zeros; same interface as Base.zeros() but faster
+_zeros(T::DataType, N::Int64) = ntuple(_ -> zero(T), N)
+_zeros(N::Int) = _zeros(Float64, N)
+
+# Return an NTuple{N, T} of ones; same interface as Base.ones() but faster
+_ones(T::DataType, N::Int64) = ntuple(_ -> one(T), N)
+_ones(N::Int) = _ones(Float64, N)
+
 ################################################################################
 #                           DifferentiationMethod
 ################################################################################
diff --git a/test/utils.jl b/test/utils.jl
index f5c26308..a42ae24f 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1,5 +1,6 @@
 @testitem "Utilities" setup=[Setup] begin
     using LinearAlgebra: norm
+    using MeshIntegrals: _units, _zeros, _ones
 
     # _KVector
     v = Meshes.Vec(3, 4)
@@ -7,7 +8,15 @@
 
     # _units
     p = Point(1.0u"cm", 2.0u"mm", 3.0u"m")
-    @test MeshIntegrals._units(p) == u"m"
+    @test _units(p) == u"m"
+
+    # _zeros
+    @test _zeros(2) == (0.0, 0.0)
+    @test _zeros(Float32, 2) == (0.0f0, 0.0f0)
+
+    # _ones
+    @test _ones(2) == (1.0, 1.0)
+    @test _ones(Float32, 2) == (1.0f0, 1.0f0)
 end
 
 @testitem "DifferentiationMethod" setup=[Setup] begin