Enforce int types in math.factorial

[jakelishman]

Using integer-like floats in math.factorial is deprecated as of Python 3.9.

Glancing over the rest of the code, I'm fairly sure math.factorial is only called on floats formed by (e.g.) N / 2 + 0.5, which is guaranteed be an integer for all odd N integers, but to be safe I inserted the same test that math.factorial will do as well.

By the way: depending on how accurate you actually want/need to be with your degeneracy calculations, a common way to deal with these binomial quantities is to work in the logarithmic space --

def factln(x):
    return scipy.special.gammaln(x + 1)

def energy_degeneracy(N, m):
    return int(np.exp(Decimal(factln(N) - factln(N/2 + m) - factln(N/2 - m))))

This is pretty much guaranteed to be faster, but a little less precise; double-precision floats have ~15 decimal digits of precision compared to the Decimal default of 28. You have to be careful that the np.exp call doesn't overflow (unless you can use the number in logarithmic space as well), but you can just use a single Decimal instance like I did if it really matters to you to have huuuuge numbers output. I suspect it doesn't, since you multiply it by a float right after, which will overflow to inf.

Timings:

from math import factorial
from decimal import Decimal
import numpy as np
from scipy.special import gammaln

def degeneracy_all_decimal(N, m):
    return int(Decimal(factorial(N)) / (Decimal(factorial(int(N/2+m))) * Decimal(factorial(int(N/2-m)))))

def degeneracy_log_then_decimal(N, m):
    return int(np.exp(Decimal(gammaln(N+1) - gammaln(N/2+m+1) - gammaln(N/2-m+1))))

def degeneracy_log(N, m):
    return int(np.exp(gammaln(N+1) - gammaln(N/2+m+1) - gammaln(N/2-m+1)))

In [3]: %timeit degeneracy_all_decimal(1000, 0)
   ...: %timeit degeneracy_log_then_decimal(1000, 0)
   ...: %timeit degeneracy_log(1000, 0)
831 µs ± 1.91 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
31.9 µs ± 190 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
4.63 µs ± 86.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

[coveralls]

Coverage increased (+0.003%) to 64.332% when pulling 8aa0cb5 on jakelishman:fix-piqs-factorial into c5be0ee on qutip:master.

[quantshah]

I think we do not have very large numbers as output and the log degeneracy formula should work for us. Thanks for the neat suggestion and the benchmark. Let me see how much improvement we get in the construction with this faster degeneracy computation. I am going to merge this and discuss the change in a new PR.

[jakelishman]

I think you won't get a huge speed up unless you're dealing with pretty large numbers - if I remember correctly from looking through the Cython backing as part of the dev.major swap-over, you're most likely dominated by interaction with the Python runtime in your tight Cython loops. That said, for huge numbers it may be the limit - arbitrary-precision integer factorials have exponential time complexity.