How to Handle Overflow in Double Scalars When Working with Very Large Numbers in Python?

I'm working on a data that involves computations with very large numbers, and I'm encountering an issue with numerical overflow and plotting. The problem arises when I increase the value of a parameter 𝑁 (at the beginning of the code).

For example, with 𝑁=5, my code runs perfectly without any errors and the curve is generated. However, as I gradually increase 𝑁, parts of the plotted curve start to disappear. By the time 𝑁 reaches 20, the curve disappears entirely. the problem is that I should use values for N=1000 or 2000.

the following error occurs RuntimeWarning: overflow encountered in double_scalars,

Here’s the python code:

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()

Given that you need to handle numbers larger than the largest possible float, I would suggest using mpmath for more of this calculation.

Example:

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 = 2000 # 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= mpmath.power(L1, (N*ms))/ mpmath.gamma(N*ms)
    cons2= mpmath.power(L2, (f))/ mpmath.gamma(f)

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

    sup4_num= mpmath.power(L1, (N*ms))*mpmath.gamma(f+1   +  N*ms)
    sup4_den= (f+1)*mpmath.power(( 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 / mpmath.gamma(f)

    Z2= cons2* sup4_num / sup4_den * sup4 / mpmath.gamma(N*ms)   
  
    y = mpmath.log1p(Z1 + Z2) / mpmath.log(2)

    y_axis_points.append(float(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()