diff --git a/Project.toml b/Project.toml
index bc45575..d0a3a53 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 [compat]
 ExplicitImports = "1.14.0"
 JET = "0.9, 0.10, 0.11"
-PrecompileTools = "1"
+PrecompileTools = "1.2"
 julia = "1.10"
 
 [extras]
diff --git a/src/LightweightStats.jl b/src/LightweightStats.jl
index f6a065e..d53a95f 100644
--- a/src/LightweightStats.jl
+++ b/src/LightweightStats.jl
@@ -426,15 +426,17 @@ function cov(X::AbstractMatrix; dims::Int = 1, corrected::Bool = true)
         n == 0 && return fill(oftype(real(zero(eltype(X))) / 1, NaN), p, p)
 
         means = vec(mean(X; dims = 1))
-        C = zeros(float(real(eltype(X))), p, p)
-
-        # Center the data once using broadcasting
-        X_centered = X .- means'
+        Tc = float(eltype(X))
+        C = zeros(Tc, p, p)
 
         for i in 1:p
+            mi = means[i]
             for j in i:p
-                # Use views and broadcasting for column operations
-                s = sum(view(X_centered, :, i) .* _conj(view(X_centered, :, j)))
+                mj = means[j]
+                s = zero(Tc)
+                for k in 1:n
+                    s += (X[k, i] - mi) * _conj(X[k, j] - mj)
+                end
                 C[i, j] = corrected ? s / (n - 1) : s / n
                 if i != j
                     C[j, i] = _conj(C[i, j])
@@ -443,19 +445,21 @@ function cov(X::AbstractMatrix; dims::Int = 1, corrected::Bool = true)
         end
         return C
     elseif dims == 2
-        n, p = size(X')
+        p, n = size(X)
         n == 0 && return fill(oftype(real(zero(eltype(X))) / 1, NaN), p, p)
 
         means = vec(mean(X; dims = 2))
-        C = zeros(float(real(eltype(X))), p, p)
-
-        # Center the data once using broadcasting
-        X_centered = X .- means
+        Tc = float(eltype(X))
+        C = zeros(Tc, p, p)
 
         for i in 1:p
+            mi = means[i]
             for j in i:p
-                # Use views and broadcasting for row operations
-                s = sum(view(X_centered, i, :) .* _conj(view(X_centered, j, :)))
+                mj = means[j]
+                s = zero(Tc)
+                for k in 1:n
+                    s += (X[i, k] - mi) * _conj(X[j, k] - mj)
+                end
                 C[i, j] = corrected ? s / (n - 1) : s / n
                 if i != j
                     C[j, i] = _conj(C[i, j])
diff --git a/test/regression_tests.jl b/test/regression_tests.jl
index 897718b..6e6051a 100644
--- a/test/regression_tests.jl
+++ b/test/regression_tests.jl
@@ -208,6 +208,19 @@ using Random
         end
     end
 
+    @testset "cov - complex matrices" begin
+        # Off-diagonal covariances are genuinely complex here, so the result
+        # matrix must have a complex eltype (regression for an InexactError when
+        # the result was allocated as a real matrix).
+        Xc = ComplexF64[1 + 1im 2 + 3im; 3 - 1im 4 + 2im; 5 + 0im 6 - 1im]
+        for dims in (1, 2)
+            C_lw = LightweightStats.cov(Xc; dims = dims)
+            C_st = Statistics.cov(Xc; dims = dims)
+            @test eltype(C_lw) <: Complex
+            @test C_lw ≈ C_st
+        end
+    end
+
     @testset "cor - correlation vectors" begin
         x = randn(30)
         y = randn(30)