extend factorial(n) beyond n=170
Reported by Samuel GOUGEON (@sgougeon)
BUG DESCRIPTION:
----------------
Presently, any input argument n>170 returns %inf. This is definitely not informative.
This means that factorial() can actually work only for the 170 first integers.
This is poor.
I propose to extend factorial() for n up to 1e13.
To do that, 2 new output arguments are proposed:
[f, p] = factorial()
[f, p, m] = factorial()
where
* p are the powers-of-10 of the results
* m are the mantissae of the results, in ]10, 1].
The general Stirling series can be used to approximate the results.
The log10 of its factors is considered to compute p and then m.
Then,
171! is computed within a relative error < 2e-14
200! to 900! are computed within a relative error < 7e-13
1000! to 10000! are computed within a relative error < 1e-11
20000! to 100000! are computed within a relative error < 3e-11
1e6! to 1e13! are computed within a relative error increasing from 3e-9 to 0.07
Beyond, the value of mantissa becomes more uncertain.
Beyond 2e14!, even the power-of-10 becomes uncertain, since it reaches values > 1/%eps
HOW TO REPRODUCE THE BUG:
-------------------------
factorial(171)