diff --git a/Project.toml b/Project.toml
index 36a8c43..1d17f04 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,6 +38,10 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
+QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
+QuanticsGrids = "634c7f73-3e90-4749-a1bd-001b8efc642d"
+SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d"
 
 [targets]
-test = ["Test", "Random", "ReTestItems", "Aqua", "PartitionedMPSs"]
+test = ["Test", "Random", "ReTestItems", "Aqua", "PartitionedMPSs", "TensorCrossInterpolation", "QuanticsTCI", "QuanticsGrids", "SparseIR"]
diff --git a/src/fouriertransform.jl b/src/fouriertransform.jl
index f9cf2c7..40540f1 100644
--- a/src/fouriertransform.jl
+++ b/src/fouriertransform.jl
@@ -19,7 +19,7 @@ We denote the input and output MPS's by ``X`` and ``Y``, respectively.
 * ``Y(y_N, ..., y_1) = Y_1(y_N) ... Y_N (y_1)``.
 
 """
-function _qft(sites; cutoff::Float64=1e-14, sign::Int=1)
+function _qft(sites; cutoff::Float64=1e-25, sign::Int=1)
     if !all(dim.(sites) .== 2)
         error("All siteinds for qft must has Qubit tag")
     end
@@ -108,8 +108,8 @@ function fouriertransform(M::MPS;
         tag::String="",
         sitessrc=nothing,
         sitesdst=nothing,
-        originsrc::Float64=0.0,
-        origindst::Float64=0.0,
+        originsrc::Real=0.0,
+        origindst::Real=0.0,
         cutoff_MPO=1e-25, kwargs...)
     sites = siteinds(M)
     sitepos, target_sites = _find_target_sites(M; sitessrc=sitessrc, tag=tag)
@@ -122,12 +122,15 @@ function fouriertransform(M::MPS;
         error("Invalid target_sites")
     end
 
+    p_back = precision(BigFloat)
+    setprecision(BigFloat, 256)
+
     # Prepare MPO for QFT
     MQ_ = _qft(target_sites; sign=sign, cutoff=cutoff_MPO)
     MQ = matchsiteinds(MQ_, sites)
 
     # Phase shift from origindst
-    M_result = phase_rotation(M, sign * 2π * origindst / (2.0^length(sitepos));
+    M_result = phase_rotation(M, sign * 2π * BigFloat(origindst) / (BigFloat(2)^length(sitepos));
         targetsites=target_sites, kwargs...)
 
     # Apply QFT
@@ -139,10 +142,13 @@ function fouriertransform(M::MPS;
     end
 
     # Phase shift from originsrc
-    M_result = phase_rotation(M_result, sign * 2π * originsrc / (2.0^length(sitepos));
+    M_result = phase_rotation(M_result, sign * 2π * BigFloat(originsrc) / (BigFloat(2)^length(sitepos));
         targetsites=sitesdst, kwargs...)
 
-    M_result *= exp(sign * im * 2π * originsrc * origindst / 2.0^length(sitepos))
+    tmp = Float64(mod(sign * 2π * BigFloat(originsrc) * BigFloat(origindst) / BigFloat(2)^length(sitepos), 2 * π))
+    M_result *= exp(im * tmp)
+
+    setprecision(BigFloat, p_back)
 
     return M_result
 end
diff --git a/src/imaginarytime.jl b/src/imaginarytime.jl
index 39827d5..5a2257c 100644
--- a/src/imaginarytime.jl
+++ b/src/imaginarytime.jl
@@ -20,21 +20,27 @@ function to_wn(stat::Statistics, gtau::MPS, beta::Float64; sitessrc=nothing, tag
         sitesdst=nothing, kwargs...)::MPS
     sitepos, _ = _find_target_sites(gtau; sitessrc=sitessrc, tag=tag)
     nqbit_t = length(sitepos)
-    originwn = 0.5 * (-2.0^nqbit_t + _stat_shift(stat))
+    p_back = precision(BigFloat)
+    setprecision(BigFloat, 256)
+    originwn = 0.5 * (-BigFloat(2)^BigFloat(nqbit_t) + _stat_shift(stat))
     giv = fouriertransform(gtau; tag=tag, sitessrc=sitessrc, sitesdst=sitesdst,
         origindst=originwn, kwargs...)
     giv *= (beta * 2^(-nqbit_t / 2))
+    setprecision(BigFloat, p_back)
     return giv
 end
 
 function to_tau(stat::Statistics, giv::MPS, beta::Float64; sitessrc=nothing, tag="",
         sitesdst=nothing, kwargs...)::MPS
     sitepos, _ = _find_target_sites(giv; sitessrc=sitessrc, tag=tag)
+    p_back = precision(BigFloat)
+    setprecision(BigFloat, 256)
     nqbit_t = length(sitepos)
-    originwn = 0.5 * (-2.0^nqbit_t + _stat_shift(stat))
+    originwn = 0.5 * (-BigFloat(2)^BigFloat(nqbit_t) + _stat_shift(stat))
     gtau = fouriertransform(giv; sign=-1, tag=tag, sitessrc=sitessrc, sitesdst=sitesdst,
         originsrc=originwn, kwargs...)
     gtau *= ((2^(nqbit_t / 2)) / beta)
+    setprecision(BigFloat, p_back)
     return gtau
 end
 
diff --git a/src/transformer.jl b/src/transformer.jl
index 9d6b00b..9d03d82 100644
--- a/src/transformer.jl
+++ b/src/transformer.jl
@@ -116,7 +116,7 @@ end
 """
 Multiply by exp(i θ x), where x = (x_1, ..., x_R)_2.
 """
-function phase_rotation(M::MPS, θ::Float64; targetsites=nothing, tag="", kwargs...)::MPS
+function phase_rotation(M::MPS, θ::Real; targetsites=nothing, tag="", kwargs...)::MPS
     transformer = phase_rotation_mpo(siteinds(M), θ; targetsites=targetsites, tag=tag)
     _apply(transformer, M; kwargs...)
 end
@@ -126,18 +126,22 @@ Create an MPO for multiplication by `exp(i θ x)`, where `x = (x_1, ..., x_R)_2`
 
 `sites`: site indices for `x_1`, `x_2`, ..., `x_R`.
 """
-function phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Float64;
+function phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real;
         targetsites=nothing, tag="")::MPO where {T}
     _, target_sites = _find_target_sites(sites; sitessrc=targetsites, tag=tag)
     transformer = _phase_rotation_mpo(target_sites, θ)
     return matchsiteinds(transformer, sites)
 end
 
-function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Float64)::MPO where {T}
+function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real)::MPO where {T}
     R = length(sites)
     tensors = [ITensor(true) for _ in 1:R]
+    θ_mod = mod(θ, 2π)
+    p_back = precision(BigFloat)
+    setprecision(BigFloat, 256)
     for n in 1:R
-        tensors[n] = ITensors.SiteTypes.op("Phase", sites[n]; ϕ=θ * 2^(R - n))
+        tmp = Float64(mod(θ_mod * BigFloat(2)^BigFloat(R - n), 2 * π))
+        tensors[n] = ITensors.SiteTypes.op("Phase", sites[n]; ϕ=tmp)
     end
     links = [Index(1, "Link,l=$l") for l in 1:(R - 1)]
     tensors[1] = ITensor(
@@ -149,6 +153,7 @@ function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Float64)::MPO
     tensors[end] = ITensor(
         Array(tensors[end], sites[end]', sites[end]), links[end], sites[end], sites[end]')
 
+    setprecision(BigFloat, p_back)
     return MPO(tensors)
 end
 
diff --git a/test/imaginarytime_tests.jl b/test/imaginarytime_tests.jl
index 48d783a..cd785c7 100644
--- a/test/imaginarytime_tests.jl
+++ b/test/imaginarytime_tests.jl
@@ -1,11 +1,18 @@
 @testitem "imaginarytime_tests.jl/imaginarytime" begin
     using Test
     using Quantics
+    import Quantics
     using ITensors: siteinds, Index
     using ITensors.SiteTypes: op
     import ITensors
     import SparseIR: Fermionic, Bosonic, FermionicFreq, valueim
 
+    import ITensorMPS: MPS, onehot
+    import QuanticsGrids as QG
+    import TensorCrossInterpolation as TCI
+    import QuanticsTCI: quanticscrossinterpolate
+    import TCIITensorConversion
+
     function _test_data_imaginarytime(nbit, β)
         ω = 0.5
         N = 2^nbit
@@ -80,6 +87,52 @@
         @test maximum(abs, (gtau - gtau_smpl)[trunc(Int, 0.2 * nτ):trunc(Int, 0.8 * nτ)]) <
               1e-2
     end
+
+
+    @testset "ImaginaryTimeFT.to_tau with large R" begin
+        function fermionic_wn(n, β)
+            return (2 * n + 1) * π / β
+        end
+
+        function inv_iwn(n, β; ϵ=0.0)
+            return 1.0 / (im * fermionic_wn(n, β) - ϵ)
+        end
+
+        _evaluate(Ψ::MPS, sites, index::Vector{Int}) = only(reduce(
+            *, Ψ[n] * onehot(sites[n] => index[n]) for n in 1:length(Ψ)))
+
+        β = 100.0
+        R = 50
+        N = 2^R
+        N_half = 2^(R - 1)
+        tol = 1e-16
+        maxdim_TCI = 100
+        maxdim_contract = 1000
+        cutoff_mpo = 1e-30
+        cutoff_contract = 1e-30
+        τ_check = 0.99 * β
+
+        ngrid = QG.InherentDiscreteGrid{1}(R, -N_half)
+        τgrid = QG.DiscretizedGrid{1}(R, 0, β)
+
+        sitesiω = [Index(2, "Qubit, iω=$n") for n in 1:R]
+        sitesτ = [Index(2, "Qubit, τ=$n") for n in 1:R]
+
+        inv_iwn_tci(n) = inv_iwn(n, β; ϵ= 0.0)
+        qtci2, ranks2, errors2 = quanticscrossinterpolate(
+            ComplexF64, inv_iwn_tci, ngrid; tolerance=tol, maxbonddim=maxdim_TCI)
+
+        inv_iwn_tt = TCI.TensorTrain(qtci2.tci)
+        iwmps = MPS(inv_iwn_tt; sites=sitesiω)
+
+        fourier_inv_iw = Quantics.to_tau(
+            Fermionic(), iwmps, β; tag="iω", sitesdst=sitesτ, cutoff_MPO=cutoff_mpo,
+            cutoff=cutoff_contract, maxdim=maxdim_contract, alg="naive")
+
+        gtau_reconst = _evaluate(fourier_inv_iw, reverse(sitesτ), reverse(QG.origcoord_to_quantics(τgrid, τ_check)))
+
+        @test abs(gtau_reconst - (-1/2)) < 1e-11
+    end
 end
 
 @testitem "imaginarytime_tests.jl/poletomps" begin
@@ -112,4 +165,5 @@ end
         gtauref = gtauf.(LinRange(0, β, 2^nqubit + 1)[1:(end - 1)])
         @test maximum(abs, gtauref .- gtauvec) < 1e-14
     end
+
 end