diff --git a/src/SparseIR.jl b/src/SparseIR.jl
index 8b141db..dbab101 100644
--- a/src/SparseIR.jl
+++ b/src/SparseIR.jl
@@ -13,14 +13,19 @@ export AugmentedBasis, TauConst, TauLinear, MatsubaraConst
 export TauSampling, MatsubaraSampling, evaluate, fit, evaluate!, fit!,
     MatsubaraSampling64F, MatsubaraSampling64B, TauSampling64, sampling_points
 
-using MultiFloats: Float64x2
+using MultiFloats: MultiFloats, Float64x2, _call_big
+# If we call MultiFloats.use_bigfloat_transcendentals() like MultiFloats
+# recommends, we get an error during precompilation:
+for name in (:exp, :sinh, :cosh)
+    eval(:(function Base.$name(x::Float64x2)
+        Float64x2(_call_big($name, x, precision(Float64x2) + 20))
+    end))
+end
 using LinearAlgebra: LinearAlgebra, cond, dot, svd, SVD, QRIteration, mul!
 using QuadGK: gauss, quadgk
 using Bessels: sphericalbesselj
 using PrecompileTools
 
-include("_multifloat_funcs.jl")
-
 include("_linalg.jl")
 include("_roots.jl")
 include("_specfuncs.jl")
diff --git a/src/_multifloat_funcs.jl b/src/_multifloat_funcs.jl
deleted file mode 100644
index 33fec50..0000000
--- a/src/_multifloat_funcs.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-# FIXME: These are piracy, but needed to make MultiFloats work for us.
-function Base.sinh(x::Float64x2)
-    if iszero(x) || isnan(x) || isinf(x)
-        return x
-    end
-
-    return if abs(x) > log(2)
-        exp_x = exp(x)
-        exp_minus_x = 1 / exp_x
-        (exp_x - exp_minus_x) / 2
-    else
-        term = x
-        sum = x
-        n = 2
-        x² = x^2
-        while true
-            term *= x² / (n * (n + 1))
-            sum_new = sum + term
-            sum == sum_new && break
-            sum = sum_new
-            n += 2
-        end
-        sum
-    end
-end
-
-function Base.cosh(x::Float64x2)
-    iszero(x) && return one(x)
-    isnan(x) && return x
-    isinf(x) && return abs(x)
-
-    return if abs(x) > log(2)
-        exp_x = exp(x)
-        exp_minus_x = 1 / exp_x
-        (exp_x + exp_minus_x) / 2
-    else
-        term = one(x)
-        sum = one(x)
-        n = 2
-        x² = x^2
-        while true
-            term *= x² / ((n - 1) * n)
-            sum_new = sum + term
-            sum == sum_new && break
-            sum = sum_new
-            n += 2
-        end
-        sum
-    end
-end
diff --git a/test/_multifloat_funcs.jl b/test/_multifloat_funcs.jl
index faae4c3..a9e850d 100644
--- a/test/_multifloat_funcs.jl
+++ b/test/_multifloat_funcs.jl
@@ -18,12 +18,6 @@ logrange(x1, x2, length) = (exp10(y) for y in range(log10(x1), log10(x2); length
         @test cosh(x) isa Float64x2
     end
 
-    @testset "allocations" begin
-        x = rand(xx)
-        @test iszero(@allocated sinh(x))
-        @test iszero(@allocated cosh(x))
-    end
-
     @testset "sinh(x) where x = $x" for x in xx
         @test sinh(x)≈sinh(big(x)) rtol=eps(Float64x2)
     end