Use muladd in aleph calculation in _quantile to avoid some rounding

errors.

In some cases when aleph is an integer (in exact) arithmetic, the
calculation of aleph via m can introduce some rounding errors that
makes the value non-integer. Using muladd avoids this in some cases.

src/Statistics.jl

@@ -1035,8 +1035,9 @@ end
     @assert n > 0 # this case should never happen here
 
     m = alpha + p * (one(alpha) - alpha - beta)
-    aleph = n*p + oftype(p, m)
-    j = clamp(trunc(Int, aleph), 1, n-1)
+    # Using muladd here avoids some rounding errors when aleph is an integer
+    aleph = muladd(n, p, m)
+    j = clamp(trunc(Int, aleph), 1, n - 1)
     γ = clamp(aleph - j, 0, 1)
 
     if n == 1

test/runtests.jl

@@ -796,6 +796,11 @@ end
     @test quantile(v, 0.8, alpha=1.0, beta=1.0) ≈ 10.6
     @test quantile(v, 1.0, alpha=0.0, beta=0.0) ≈ 21.0
     @test quantile(v, 1.0, alpha=1.0, beta=1.0) ≈ 21.0
+
+    @testset "avoid some rounding" begin
+        @test [quantile(1:10, i/9) for i in 0:9] == 1:10
+        @test [quantile(1:14, i/13) for i in 0:13] == 1:14
+    end
 end
 
 # StatsBase issue 164