I have added a new py_Fact function to the py_Scipy module, that returns the factorial of an integer as a float if it is up to 15 digits long, or a string for longer numbers.
The revised code, and the example shown below in the py_Special-Dist.xlsb spreadsheet, can be downloaded from:
@xl_func(category = "Scipy", help_topic="https://docs.scipy.org/doc/scipy/reference/special.factorial.html!0")
@xl_arg('n', 'numpy_array<int>', ndim=1)
@xl_arg('exact', 'bool')
@xl_return('numpy_array<var>')
def py_fact(n, exact = False):
res = scipy.special.factorial(n, exact = exact)
if exact:
for i in range(0, res.shape[0]):
if res[i] > 1E15: res[i] = str(res[i])
return res
If the “exact” argument is set to True, or omitted, the Python function returns a long integer of the required length, but if this is returned to Excel as a number it would be converted to a float, so all the results greater than 1E15 are converted to a string so no digits are lost, and the resulting array is returned to Excel as a variant array.
In the example below the last 13 digits of the strings have been compared with the last 13 digits of the values returned as floats. Because the factorial numbers always have trailing zeros the integer and float results are exactly equal up to 22!, after which the difference rapidly starts to increase.

Recently someone told me that they were going to a birthday party of a neighbour who was turning 106! and I can now point out that this is equal to:
114628056373470835453434738414834942870388487424139673389282723476762012382449946252660360871841673476016298287096435143747350528228224302506311680000000000000000000000000.