同一Julia程序在不同硬件上数值结果存在差异的技术咨询

我使用同一版本的Julia在两台不同PC上运行同一程序，却得到略有差异的数值结果。以下是两台机器的运行结果：

1  1.0000000000582077  15.124999999941792
2  1.9999999999417923  27.536972172616515
3  3.000000000523869  45.722282028989866

1  1.0  15.125
2  2.0  27.53697217302397
3  3.0  45.722282028466

问题在于这些舍入误差，我不清楚误差产生的环节。这并非格式化错误（程序完全相同），也不属于机器epsilon误差（两台机器的机器epsilon均>2.220446049250313e-16）。目前我怀疑是使用的SparseArrays、LinearAlgebra或FFTW包引入了这类误差，但无法确定。

编辑补充 #

我整理了一个最小示例（可能无需使用由稀疏数组构成的[a,a]形式数组），修改程序以定位问题，以下结果与上述略有不同，但均显示了不可忽略的误差：

using LinearAlgebra
using SparseArrays
using FFTW

const Nsites = 10000

function Opbuild1(::Val{true})
    Up = spzeros(Complex{Float64}, Nsites, Nsites)
    N = div(Nsites,2)
    ii = 1
    for p in 0:(N-1)
        Up[ii,ii] = exp(-1.0im*p^2*0.25)
        Up[ii + N, ii + N] = exp(-1.0im*(p - N)^2*0.25)
        ii += 1
    end
    return kron(sparse(Matrix(1.0I, 2, 2)), Up)
end

function Opbuild2(::Val{true})
    pop = spzeros(Float64, Nsites, Nsites)
    N = div(Nsites,2)
    ii = 1
    for p in 0:(N-1)
        pop[ii,ii] = p
        pop[ii+N,ii+N] = (p-N)
        ii += 1
    end
    return kron(sparse(Matrix(1.0I, 2, 2)), pop)
end

function initst1(BoolN, irealiz)
    psi = spzeros(Complex{Float64}, Nsites)
    a = range(680,stop=783,step=1)
    x = a[irealiz]
    psi[x] = 1.0
    psi = kron([1.0,1.0im]/sqrt(2), psi)
    return psi, x
end

function ifftcustom(array)
    array = Array(array)
    array1 = copy(array)
    array1[1:Nsites] = ifft(array[1:Nsites])
    array1[Nsites+1:end] = ifft(array[Nsites+1:end])
    return array1/norm(array1)
end

function fftcustom(array)
    array = Array(array)
    array1 = copy(array)
    array1[1:Nsites] = fft(array[1:Nsites])
    array1[Nsites+1:end] = fft(array[Nsites+1:end])
    return array1/norm(array1)
end

function main()
    boolN = rem(Nsites,2) == 0
    Operator = Opbuild1(Val(boolN))
    Operator_p = Opbuild2(Val(boolN))
    #                                                                                                                                                                                                              
    psi0_inp, pinit = initst1(boolN, 1) #!                                                                                                                                                                         
    psi0_inp = round.(ifftcustom(Operator*psi0_inp), digits=15) # (*)See below
    psi0_inp = round.(fftcustom(psi0_inp), digits=15)
    psi0_inp = round.(psi0_inp, digits=15)
    @show psi0_inp'*Operator_p*psi0_inp
end

main()

差异情况如下：若将标记为(*)的代码行替换为psi0_inp = round.(ifftcustom(psi0_inp), digits=15)，两台机器的@show部分均得到679.000000000001；而使用原代码时，一台机器得到679.0000000000001 + 0.0im，另一台得到679.0 + 1.2036944776504516e-16im。我不在意虚部的1e-16量级数值（双精度下可视为0），但实部的1e-13量级差异并非双精度下的可忽略零值。

内容的提问来源于stack exchange，提问作者user2820579