use scipy.special.factorial for non-integer values

[ChristopherChudzicki]

This essential uses factorial(n)=gamma(n+1), which gives us the correct algebraic properties but let's us be a little sloppier with sampling. Similar to the approach we've taken with square roots and complex numbers.

[coveralls]

Coverage remained the same at 100.0% when pulling 303cf1fa23d98fb3c8e319ccbec77e99fb4903d5 on factorial-as-gamma+1 into 710adaa1881859d5b550e8699e46830bd3186e3a on master.

[jolyonb]

Important question: How does special.factorial go with complex arguments? Do we need to forbid that? Should we change to an actual complex gamma function? When you cast to a float, that precludes any complex result.

[jolyonb]

You've got failing builds.

Also, you need to update calc.py, around line 333.

[ChristopherChudzicki]

Re:

Also, you need to update calc.py, around line 333.

I already updated that :) The failing builds were rounding differences in a doctest.

Complex numbers ... good point.

So I've switched to gamma(z+1) because plugging in a few numbers leads me to believe scipy's factorial has a bug:

# Expected:
factorial(3.1 + 2j)     # (-3.5657350914240893+1.8577259455190758j)
gamma(3.1+2j + 1)       # (-3.5657350914240893+1.8577259455190758j)

factorial(1e-6 - 1j)    #(0.49801588208594105+0.15494930676741031j)
gamma(1e-6 - 1j + 1)    #(0.49801588208594105+0.15494930676741031j)

# unexpected:
factorial(0 - 1j)    # 0j ???
gamma(0 - 1j + 1)    #(0.4980156681183554+0.15494982830181081j)

This persists in recent versions of scipy on Python 3. Maybe if I get bored some weekend I'l try to fix it =).

Anyway, should be ready now.

[jolyonb]

Ok, sounds good to me!