Merge pull request #184 from JuliaStats/an/quantilemuladd

JuliaStats · Jan 3, 2025 · bfa5c6b · bfa5c6b
2 parents 793733e + 6bd1531
commit bfa5c6b
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 2 deletions.
diff --git a/src/Statistics.jl b/src/Statistics.jl
@@ -1035,8 +1035,10 @@ 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 fma here avoids some rounding errors when aleph is an integer
+    # The use of oftype supresses the promotion caused by alpha and beta
+    aleph = fma(n, p, oftype(p, m))
+    j = clamp(trunc(Int, aleph), 1, n - 1)
     γ = clamp(aleph - j, 0, 1)
 
     if n == 1

diff --git a/test/runtests.jl b/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