[FIXED] How to get Python ufuncs to return floats larger than the normal float size limit?

Issue

I have a NumPy array that includes a number that is a little over 260, and I am running the gamma ufunc from scipy.special on it and it currently returns inf for this entry of the array because the result is over 1e514 and the largest float Python can ordinarily handle is about 1.8e308. I have tried changing my array to the float128 type in the hope this would fix things, but no, this causes me to receive the error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: ufunc 'gamma' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Is there a way around this error, or are Python ufuncs incapable of handling all floats over 1.8e308?

Solution

ufuncs use highly optimized C code to perform the operations required. C uses 64-bit arithmetic by default, so ufuncs cannot perform more precise arithmetic than 64-bit.

I found a workaround for this problem, however, by translating my Python code to Julia and using the BigFloat type for arbitrary precision arithmetic. Python has mpmath for arbitrary precision arithmetic, but mpmath's functions cannot be run on NumPy arrays and my code heavily utilizes arrays. The BigFloat type in Julia is compatible with most functions and most functions in Julia can be run on arrays without a problem, so my Julia code ran without any major issues.

Answered By - Josh Pinto