diff --git a/test/testsuite/reductions.jl b/test/testsuite/reductions.jl
index de2ced00..dc9ebec0 100644
--- a/test/testsuite/reductions.jl
+++ b/test/testsuite/reductions.jl
@@ -25,18 +25,21 @@ test_shapes(ET, all, few) = ET in (Float32, ComplexF32) ? all : few
 @testsuite "reductions/mapreducedim!" (AT, eltypes)->begin
     @testset "$ET" for ET in eltypes
         range = ET <: Real ? (ET(1):ET(10)) : ET
+        # ten ComplexF16 values from [0, 1) multiply to a magnitude `norm` cannot resolve in
+        # Float16, so the array comparison tolerates no error; unit values keep it near one
+        prod_data(sz) = ET == ComplexF16 ? ET.(cis.(2π .* rand(sz...))) : rand(range, sz)
         for (sz,red) in test_shapes(ET, [(10,)=>(1,), (10,10)=>(1,1), (10,10,10)=>(1,1,1), (10,10,10)=>(10,10,10),
                                          (10,10,10)=>(1,10,10), (10,10,10)=>(10,1,10), (10,10,10)=>(10,10,1),
                                          (0,)=>(1,)],
                                     [(10,)=>(1,), (10,10,10)=>(10,1,10), (0,)=>(1,)])
             # mapreducedim!
             @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, sz), zeros(ET, red))
-            @test compare((A,R)->Base.mapreducedim!(identity, *, R, A), AT, rand(range, sz), ones(ET, red))
+            @test compare((A,R)->Base.mapreducedim!(identity, *, R, A), AT, prod_data(sz), ones(ET, red))
             @test compare((A,R)->Base.mapreducedim!(x->x+x, +, R, A), AT, rand(range, sz), zeros(ET, red))
 
             # reducedim!
             @test compare((A,R)->Base.reducedim!(+, R, A), AT, rand(range, sz), zeros(ET, red))
-            @test compare((A,R)->Base.reducedim!(*, R, A), AT, rand(range, sz), ones(ET, red))
+            @test compare((A,R)->Base.reducedim!(*, R, A), AT, prod_data(sz), ones(ET, red))
         end
 
         # implicit singleton dimensions
diff --git a/test/testsuite/statistics.jl b/test/testsuite/statistics.jl
index fb9edfc2..29886163 100644
--- a/test/testsuite/statistics.jl
+++ b/test/testsuite/statistics.jl
@@ -30,9 +30,11 @@ using Statistics
             @test compare(mean, AT, rand(range, 2, 2))
             @test compare(A -> mean(A; dims = 2), AT, rand(range, 2, 2))
             @test compare(A -> mean(A; dims = [1, 3]), AT, rand(range, 2, 2, 2))
-            @test compare(A -> mean(sin, A), AT, rand(range, 2, 2))
-            @test compare(A -> mean(sin, A; dims = 2), AT, rand(range, 2, 2))
-            @test compare(A -> mean(sin, A; dims = [1, 3]), AT, rand(range, 2, 2, 2))
+            # sin of `range` can cancel to near zero, where the summation order decides the
+            # result; values in [0, 1) keep every term positive
+            @test compare(A -> mean(sin, A), AT, rand(ET, 2, 2))
+            @test compare(A -> mean(sin, A; dims = 2), AT, rand(ET, 2, 2))
+            @test compare(A -> mean(sin, A; dims = [1, 3]), AT, rand(ET, 2, 2, 2))
         end
     end
 
@@ -41,8 +43,9 @@ using Statistics
         if !(ET in [Float32, Float64, Float16, ComplexF64])
             continue
         end
-        # Larger range for Float16 to avoid numerical instability
-        range = ET == Float16 ? (ET(1):ET(100)) : ET
+        # Integers for Float16, so that `cor` of two-element rows have no subnormal variances,
+        # but small ones, so that `cov` of 100 samples does not overflow
+        range = ET == Float16 ? (ET(1):ET(10)) : ET
 
         @testset "cov" begin
             s = 100