Improve numerical stability of cov (#85)

`cov(::AbstractVector, ::AbstractVector)` was less stable than other `cov`
methods and than `var` on arrays as it called `sum` on a generator, which
does not use pairwise summation. Use  `dot` instead, which is also faster,
and was used before 7deac81, where it was removed to be able to move
LinearAlgebra to a separate package (when statistics.jl was in Base).

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) =

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