I am new to Julia and have only written in python before which is most likely also reflected by my coding style.
Unfortunately, my applications require very high numerical precision which is why I am resorting to
BigInt in Julia. For some reason, this leads to a major increase in allocations besides the obvious increase in computational time.
I am wondering where my code can be improved.
I have attached the current version of my code below with example parameters that allow for execution under 1s.
using Combinatorics function gamma_mpnq(m::Int, p::Int, n::Int, q::Int)::BigFloat return sqrt(factorial(big(m)) * factorial(big(n))) / (big(2)^(p + q) * factorial(big(m - 2 * p)) * factorial(big(n - 2 * q)) * factorial(big(q)) * factorial(big(p))) end function general_binomial(alpha::Float64, k::Int) """ Implementation of the generalized binomial. """ _prod = BigFloat(1.0) for kk in 1:k _prod *= (alpha - kk + 1.) end return _prod / factorial(big(k)) end function construct_M(nmax::Int, mmax::Int)::Array """ Construct the matrix M containing all moments <C^n * C*^m> up to nmax and mmax. """ _M = zeros(BigFloat, nmax, mmax) for mm in 1:nmax for nn in 1:mmax if nn == 1 _M[nn, mm] = 1. / doublefactorial(BigInt(2*mm - 1)) elseif mm == 1 _M[nn, mm] = 1. / doublefactorial(BigInt(2*nn - 1)) else _M[nn, mm] = (((nn-1) * _M[nn-1, mm] + (mm-1) * _M[nn, mm-1]) / (2. * (nn-mm)^2 + (nn-1) + (mm-1))) end end end return _M end function maclaurin_exp(n::Int, m::Int, p::Int, q::Int, M::Array, cutoff::Int)::BigFloat """ Calculate the expectation value of ((1+C)/(1+C*))^((n-m)/2) * C^p * C*^q using the Maclaurin expansion and the matrix containing all moments M. """ inner::BigFloat = 0.0 @views for nn in 1:cutoff for mm in 1:cutoff inner += (general_binomial((n-m)/2., mm-1) * general_binomial((m-n)/2., nn-1) * M[q+nn, p+mm]) end end return inner end function Hmn(m::Int, n::Int, M::Array, cutoff::Int)::Float64 """ Calculates a single entry of the Hmn matrix. """ ppmax = div(m, 2) qqmax = div(n, 2) outer = 0.0 for pp in 0:ppmax for qq in 0:qqmax outer += (gamma_mpnq(m, pp, n, qq) * maclaurin_exp(n, m, pp, qq, M, cutoff)) # println(pp, " ", qq, " ", outer) end end return outer end dim = 128 cutoff = 64 M = construct_M(dim+cutoff, dim+cutoff) @btime Hmn(12, 2, M, cutoff) ```