Big numbers in lambdify

Question

I have (automatically generated by SymPy) expression, including sqrt function with big numbers under sqrt and small multiplier in front of it. So overall result must be within float range, but value under sqrt - is not. I coonvert this expression to Python function using lambdify command. Call of this function gets exception:

from sympy import *
t = symbols("t")
k = 10
f = 2 * 10 ** (- k) * sqrt(10 ** (2 * k) * t ** 2 + 1)
print(f)
F = lambdify(t, f)
t0 = 10 ** 10
T = np.arange(t0, 2 * t0, t0 / 4)
print(T)
F(T)

Output:

2.0e-10*sqrt(100000000000000000000*t**2 + 1)
[1.00e+10 1.25e+10 1.50e+10 1.75e+10]
AttributeError: 'float' object has no attribute 'sqrt'

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_15740/1035914544.py in <module>
      8 T = np.arange(t0, 2 * t0, t0 / 4)
      9 print(T)
---> 10 F(T)

<lambdifygenerated-11> in _lambdifygenerated(t)
      1 def _lambdifygenerated(t):
----> 2     return 2.0e-10*sqrt(100000000000000000000*t**2 + 1)

TypeError: loop of ufunc does not support argument 0 of type float which has no callable sqrt method

For k = 2 code works properly:

0.02*sqrt(10000*t**2 + 1)
[1.00e+10 1.25e+10 1.50e+10 1.75e+10]
array([2.0e+10, 2.5e+10, 3.0e+10, 3.5e+10])

Is there any way to fix this problem without manually rewriting expression?

UPD: Looks like it is a problem of NumPy:

import numpy as np
k = 10
def F1(t):
    return np.sqrt( (10 ** (- k)) ** 2 * 10 ** (2 * k) * t ** 2 + 1)
def F2(t):
    return 10 ** (- k) * np.sqrt(10 ** (2 * k) * t ** 2 + 1)
print(F1(10 ** 5))
print(F2(10 ** 5))

First call works, second call - not!

Answer 1

Converting argument of np.sqrt to numpy.double solves problem:

def Sqrt(x):
  return np.sqrt(np.double(x))
F = lambdify(t, f, {'sqrt': Sqrt})

Answer 2

Because of the large multiplier, the np.sqrt argument is object dtype array:

In [3]: 100000000000000000000 * T**2
Out[3]: array([1e+40, 1.5625e+40, 2.25e+40, 3.0625e+40], dtype=object)

With object dtype arrays, numpy iterates (at list comprehension speed), applying a 'method' to each element. In effect

 1e+40.sqrt() etc

Hence the no method error.

Your fix:

In [3]: np.double(100000000000000000000 * T**2)
Out[3]: array([1.0000e+40, 1.5625e+40, 2.2500e+40, 3.0625e+40])

In [4]: np.sqrt(_)
Out[4]: array([1.00e+20, 1.25e+20, 1.50e+20, 1.75e+20])

or

In [6]: np.sqrt((100000000000000000000 * T**2).astype(float))
Out[6]: array([1.00e+20, 1.25e+20, 1.50e+20, 1.75e+20])