Binomial overflow when constructing spin coherent state with large S

[saulpila]

The code

b=SpinBasis(40)
coherentspinstate(b,0.1,0.2)

throws OverflowError: binomial(80, 21) overflows

There is an easy solution to this. When dealing with potentially large numbers, it is better to operate their logarithms and exponentiate at the end. This should fix it:

using SpecialFunctions
logbinomial(N::Integer,n::Integer) = logfactorial(N) - logfactorial(n) - logfactorial(N-n)

function coherentspinstate(b::SpinBasis, theta::Real, phi::Real)
    result = Ket(b)
    data = result.data
    N = length(b)-1
    expphi = exp(0.5im*phi)
    expphi_con = conj(expphi)
    if theta==0 #prevents a log(0) below
        data[1]=expphi_con^N
    else
        sinth = sin(0.5theta)
        costh = cos(0.5theta)
        logcosthexpphi_con=log(costh*expphi_con)
        logsinthexpphi=log(sinth*expphi)
        @inbounds for n=0:N
            data[n+1] =  exp(logbinomial(N, n)/2.0  + (N-n)*logcosthexpphi_con + logsinthexpphi*n)
        end
    end
    return result
 end

[david-pl]

Thanks for the report, the current implementation really shouldn't use binomial like that. The logbinomial you use is a pretty cool solution. Alternatively, we can also recursively multiply to immediately apply the sqrt by using binomial(N,n) == (N-n+1)/n *binomial(N,n-1) This should also avoid overflowing.

function coherentspinstate(b::SpinBasis, theta::Real, phi::Real)
    result = Ket(b)
    data = result.data
    N = length(b)-1
    sinth = sin(0.5theta)
    costh = cos(0.5theta)
    expphi = exp(0.5im*phi)
    expphi_con = conj(expphi)
    # avoid 1/n if n == 0
    data[1] = (costh * expphi_con)^N
    binom_sqrt = 1.0
    @inbounds for n=1:N
        binom_sqrt *= sqrt(N - n + 1) / sqrt(n)
        data[n+1] = binom_sqrt * (sinth*expphi)^n * (costh*expphi_con)^(N-n)
    end
    return result
end

Eventually the factor binom_sqrt could still get pretty large. We could even go a step further though and build the entire data recursively like so

# special handling for theta==pi
if iszero(costh)
    data[end] = 1
    return data
end
# avoid 1/n if n == 0
data[1] = (costh * expphi_con)^N
factor = sinth*expphi^2 / costh
@inbounds for n=1:N
    data[n+1] = sqrt(N - n + 1) / sqrt(n) * factor * data[n]
end

I guess that one only overflows if the exp you use also overflows. Not sure if I got that special handling for costh == 0 right here and it could be rewritten a bit more nicely.

The solutions here might be a bit more simple than yours, but I'm not sure which one is best to use.

[saulpila]

That works too! The only argument I might have in favor of logarithms is that in principle they allow to keep better numerical precision when the factors like sinth^n are very small, but to be honest I could not see any difference even for S=1000.

It seems that your second option would be the most efficient and it solves the issue quite nicely :)

[david-pl]

@saulpila looks like @goerz beat me to it: qojulia/QuantumOptics.jl#327 ;)

[goerz]

Wow, I totally missed this issue!