Use recursive formula for binomial coefficient in CSS construction

[goerz]

This uses a recursive formula for the binomial coefficient in the construction of a coherentspinstate in order to avoid the overflow associated with the built-in binomial for larger values of S

Closes #326

[goerz]

Actually, you could still make this faster by a factor of 2-8 (depending on the size of the Hilbert space) and numerically more stable by using recursion not just of the binomial coefficient, but also for the rest of the equation. I can add another commit for that...

[goerz]

Ok, should be ready now.

The performance benefit might not seem very relevant just for initializing a state, but it really speeds things up when you use the routine for creating the usual visualizations (especially animations) of spin-squeezed states like this:

where every point on the sphere is the overlap with the CSS pointing in that direction.

[codecov[bot]]

Codecov Report

Merging #327 (14d438d) into master (f4e996f) will increase coverage by 0.00%. The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master     #327   +/-   ##
=======================================
  Coverage   98.27%   98.27%           
=======================================
  Files          16       16           
  Lines        1331     1335    +4     
=======================================
+ Hits         1308     1312    +4     
  Misses         23       23           

Impacted Files	Coverage Δ
src/phasespace.jl	98.71% <100.00%> (+0.01%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update f4e996f...14d438d. Read the comment docs.

[david-pl]

@goerz Thanks for the PR! There's actually a duplicate issue here: https://github.com/qojulia/QuantumOpticsBase.jl/issues/43, which is also solved by your changes.

[saulpila]

By the way, now that we are talking about coherentspinstate efficiency, if you additionally change

https://github.com/qojulia/QuantumOptics.jl/blob/871183ffbec52f0a41af255415283a6c7206a24f/src/phasespace.jl#L290-L293

to

function coherentspinstate(b::SpinBasis, theta::Real, phi::Real, result::Ket=Ket(b))
        data = result.data

Then it allows to reuse the same array when doing visualizations of the Husimi function and prevents excessive memory allocation. This should speed things up in e.g. qfuncsu2 by changing

https://github.com/qojulia/QuantumOptics.jl/blob/871183ffbec52f0a41af255415283a6c7206a24f/src/phasespace.jl#L338-L340

to

 coherent = Ket(b) 
 @inbounds  for i = 0:Ntheta-1, j = 0:Nphi-1 
     coherentspinstate(b,pi-i*pi/(Ntheta-1),j*2pi/(Nphi-1)-pi,coherent)
     result[i+1,j+1] = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherent.data) 
 end 

[goerz]

I agree that could be useful, but it should probably just be a coherentspinstate! function, like so:

coherentspinstate!(Psi::Ket, theta::Real, phi::Real)

where the basis is extracted from Psi.

Then, coherentspinstate can delegate to coherentspinstate!.