diff --git a/src/Statistics.jl b/src/Statistics.jl
index 29aae101..91cf56fb 100644
--- a/src/Statistics.jl
+++ b/src/Statistics.jl
@@ -519,7 +519,7 @@ unscaled_covzm(x::AbstractVector{<:Number})    = sum(abs2, x)
 unscaled_covzm(x::AbstractVector)              = sum(t -> t*t', x)
 unscaled_covzm(x::AbstractMatrix, vardim::Int) = (vardim == 1 ? _conj(x'x) : x * x')
 
-unscaled_covzm(x::AbstractVector, y::AbstractVector) = sum(conj(y[i])*x[i] for i in eachindex(y, x))
+unscaled_covzm(x::AbstractVector, y::AbstractVector) = dot(y, x)
 unscaled_covzm(x::AbstractVector, y::AbstractMatrix, vardim::Int) =
     (vardim == 1 ? *(transpose(x), _conj(y)) : *(transpose(x), transpose(_conj(y))))
 unscaled_covzm(x::AbstractMatrix, y::AbstractVector, vardim::Int) =
diff --git a/test/runtests.jl b/test/runtests.jl
index 00cdad10..2e061ee4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -434,10 +434,19 @@ Y = [6.0  2.0;
         @inferred cov(X, Y, dims=vd, corrected=cr)
     end
 
-    @testset "floating point accuracy for `cov` of large numbers" begin
+    @testset "floating point accuracy for `cov`" begin
+        # Large numbers
         A = [4.0, 7.0, 13.0, 16.0]
         C = A .+ 1.0e10
-        @test cov(A, A) ≈ cov(C, C)
+        @test cov(A) ≈ cov(A, A) ≈
+            cov(reshape(A, :, 1))[1] ≈ cov(reshape(A, :, 1))[1] ≈
+            cov(C, C) ≈ var(C)
+
+        # Large Vector{Float32}
+        A = 20*randn(Float32, 10_000_000) .+ 100
+        @test cov(A) ≈ cov(A, A) ≈
+            cov(reshape(A, :, 1))[1] ≈ cov(reshape(A, :, 1))[1] ≈
+            var(A)
     end
 end