Как справиться с переполнением в двойных скалярах при работе с очень большими числами в Python?

Я работаю над данными, которые включают вычисления с очень большими числами, и столкнулся с проблемой числового переполнения и построения графиков. Проблема возникает, когда я увеличиваю значение параметра 𝑁 (в начале кода).

Например, при 𝑁=5 мой код работает идеально, без каких-либо ошибок, и кривая генерируется. Однако по мере того, как я постепенно увеличиваю 𝑁, части построенной кривой начинают исчезать. К тому времени, когда 𝑁 достигает 20, кривая полностью исчезает. проблема в том, что мне следует использовать значения для N=1000 или 2000.

возникает следующая ошибка RuntimeWarning: в double_scalars обнаружено переполнение,

Вот код Python:

import numpy as np
import cmath
import math
from math import e
import matplotlib.pyplot as plt
import mpmath
from scipy.special import gamma, factorial
import scipy.special as sc

N = 5 # VARIABLE OF INTEREST (SHOULD BE INCREASED)



g1, g2= 1.31928e-06, 4.947323e-07

H, ms = 1/6, 6

x_axis_points = np.arange(0, 40, 1)
E1 = 10.0 ** ((x_axis_points-30) / 10.0)/1e-11
y_axis_points = []
for QS in E1:
 
    A, B = N * 0.979405**2, N * 0.07986762

    w= (B/A)**2*(A**2/B)*(A**2/B+1)
    r= (B/A)**4*(A**2/B)*(A**2/B+1)*(A**2/B+2)*(A**2/B+3)


    f= w**2/ (r-w**2)
    c= g2**2 * QS*(r-w**2)/w
  


    BB=H *QS * g1**2
    L1, L2= 1/BB, 1/c


################### THE PROBLEM APPEARS IN THE COMMING LINES
#(GIVES INFINITY WHEN INCREASING N) ##############

    cons1= L1**(N*ms)/ sc.gamma(N*ms)  
    cons2= L2**(f)/ sc.gamma(f  )    

    sup3_num= L2**(f)*sc.gamma(N*ms+1   +  f)
    sup3_den= (N*ms+1)*( L2 + L1)**( N*ms+1 +f)               
    sup3=mpmath.hyp2f1(1, N*ms+1 + f, N*ms+1+1, L1/(L1+L2))

    sup4_num= L1**(N*ms)*sc.gamma(f+1   +  N*ms)
    sup4_den= (f+1)*( L2 + L1)**( f+1 +N*ms)                   
    sup4=mpmath.hyp2f1(1, f+1 + N*ms, f+1+1, L2/(L1+L2))



    Z1= cons1* sup3_num/sup3_den * sup3/ sc.gamma(f)

    Z2= cons2* sup4_num/sup4_den * sup4/ sc.gamma(N*ms)   
  
    y = math.log2(1+ Z1+Z2 ) 

    y_axis_points.append((y))            
plt.plot(x_axis_points, y_axis_points,    linewidth='2.5', color='blue')

plt.xlabel('x')
plt.ylabel('y')
plt.grid()

plt.xlim([0, 40])
plt.ylim([0 , 3 ])
plt.show()