diff --git a/src/intervals/arithmetic/trigonometric.jl b/src/intervals/arithmetic/trigonometric.jl
index 0005ba24e..69fb4c641 100644
--- a/src/intervals/arithmetic/trigonometric.jl
+++ b/src/intervals/arithmetic/trigonometric.jl
@@ -3,13 +3,23 @@
 # in Section 10.5.3
 
 # helper functions
+# NOTE: According to J-M Muller, "Elementary Functions: Algorithms and
+# implementations", Birkhäuser, 3rd ed (2005), chap. 11, pp 203,
+# "The first solution consists of using multiple-precision arithmetic,
+# but this may make the computation significantly slower. Moreover,
+# it is not that easy to predict the accuracy with which the computation
+# should be performed." Later (pp 210), it estimates that the precision needed
+# to compute correctly the worst case of the range reduction algorithm
+# for Float64 is about 61 "guard digits". Below we use `rem2pi` for
+# BigFloats, with at least 200 bits of precision.
+# (As noted before, `rem2pi(x, RoundNearest)`, with Float64 indeed may have inaccuracies.)
+#
+# For consistency, `_quadrantpi` uses `rem(big(x), 2, RoundNearest)`
 
 function _quadrant(x::AbstractFloat)
-    # NOTE: this algorithm may be flawed as it relies on `rem2pi(x, RoundNearest)`
-    # to yield a very tight result. This is not guaranteed by Julia, see e.g.
-    # https://github.com/JuliaLang/julia/blob/9669eecc99bc4553e28d94d7dd3dc9fd40b3bf3f/base/mpfr.jl#L845-L846
-    PI_LO, PI_HI = bounds(bareinterval(typeof(x), π))
-    r = rem2pi(x, RoundNearest)
+    @assert precision(BigFloat) > 200
+    PI_LO, PI_HI = bounds(bareinterval(BigFloat, π))
+    r = rem2pi(big(x), RoundNearest)
     r2 = 2r # should be exact for floats
     r2 ≤ -PI_HI && return 2 # [-π, -π/2)
     r2 < -PI_LO && return throw(ArgumentError("could not determine the quadrant, the remainder $r of the division of $x by 2π is lesser or greater than -π/2"))
@@ -20,7 +30,8 @@ function _quadrant(x::AbstractFloat)
 end
 
 function _quadrantpi(x::AbstractFloat) # used in `sinpi` and `cospi`
-    r = rem(x, 2) # [-2, 2], should be exact for floats
+    @assert precision(BigFloat) > 200
+    r = rem(big(x), 2, RoundNearest) # [-2, 2], should be exact for floats
     2r < -3 && return 0 # [-2π, -3π/2)
     r  < -1 && return 1 # [-3π/2, -π)
     2r < -1 && return 2 # [-π, -π/2)
diff --git a/test/interval_tests/trigonometric.jl b/test/interval_tests/trigonometric.jl
index 24180ec80..5bdceb4a2 100644
--- a/test/interval_tests/trigonometric.jl
+++ b/test/interval_tests/trigonometric.jl
@@ -307,4 +307,16 @@ end
     @test isequal_interval(sin(x), interval(-1, 1))
     @test isequal_interval(cos(x), interval(-1, 1))
     @test isequal_interval(tan(x), interval(-Inf, Inf))
+
+    setprecision(100) do
+        y = 6381956970095103 * 2.0^797 # worst case for C = pi/2 for Float64
+        yb = 6381956970095103 * big(2.0)^797
+        @test_throws AssertionError IntervalArithmetic._quadrant(y) == IntervalArithmetic._quadrant(yb)
+    end
+    setprecision(1000) do
+        y = 6381956970095103 * 2.0^797 # worst case for C = pi/2 for Float64
+        yb = 6381956970095103 * big(2.0)^797
+        @test IntervalArithmetic._quadrant(y) == IntervalArithmetic._quadrant(yb)
+    end
+
 end