How to deal with an overflow in numpy's root method?

Question

I'm writing an implementation of Coppersmith's method for root finding in modular polynomials and I need to find the roots of a polynomial with really big coefficients (these roots are over the whole range of the integers). When I try running the code I get an overflow in numpy root's method. Is there any possible solution to it? I have been using Decimal all over my code to try to avoid this problem, but as far as I know there is no way to force numpy to use Decimal internally.

Any help is very much appreciated!

This is the code in question. The line that causes the problem is in the coppersmith function, line 83.

import numpy as np
from copy import deepcopy
from math import log, ceil, floor, sqrt, gcd, prod
from decimal import *
from fpylll import IntegerMatrix, LLL


def b2f (b, X):
  return [b[i]//pow(X,i) for i in range(len(b))]

def f2b (f, X):
  return [f[i]*pow(X,i) for i in range(len(f))]

def multPol (p1, p2):
  res = [0]*(len(p1)+len(p2)-1)
  for i in range(len(p1)):
    for j in range(len(p2)):
      res[i+j] += p1[i]*p2[j]
  return res

def raisePol (f, e):
  if e == 0:
    return [1]
  if e == 1:
    return f
  if e%2 == 1:
    return multPol(raisePol(f,e-1), f)
  half = raisePol(f,e//2)
  return multPol(half,half)

# retrurns h_i,j(x) = x^j*N^i*f(x)^(m-i), for 0 <= i < m, 0 <= j < d
def geth (N, f, i, j, m, fix_len):
  # mod = N**(m-i)
  # h = [(N**i)*(x%mod) for x in raisePol(f,m-i)]
  h = [(N**i)*x for x in raisePol(f,m-i)]
  remainder = fix_len - len(h) - j
  if remainder < 0:
    raise Exception(f"fix_len too short when generating h_i={i},j={j}")
  return ([0]*j) + h + ([0]*remainder)

def buildB (f, X, N, m, t):
  d = len(f)-1
  B = []
  fix_len = d*m+t
  for i in range(m, 0, -1):
    for j in range(d):
      B.append(f2b(geth(N,f,i,j,m, fix_len), X))
  for j in range(t):
    B.append(f2b(geth(N,f,0,j,m,fix_len), X))
  return B

# using Horner's method
def evalP (p, x):
  y = p[-1]
  for a in p[0:-1][::-1]:
    y = a + x*y
  return y

def printPol (g):
  l = []
  for i in range(len(g)):
    if g[i] != 0:
      if i == 0:
        l.append(f"{hex(g[i])}")
      elif g[i] == 1:
        l.append(f"x^{i}")
      else:
        l.append(f"{hex(g[i])}*x^{i}")
  return " + ".join(l[::-1])

# based on https://www.math.auckland.ac.nz/~sgal018/crypto-book/ch19.pdf
# and https://www.cits.ruhr-uni-bochum.de/imperia/md/content/may/paper/intro_to_coppersmiths_method.pdf
def coppersmith (f, N, X, m, t):
  if f[-1] != 1:
    raise Exception("Polynomial f(x) is not monic")
  d = len(f)-1
  B = buildB(f,X,N,m,t)
  # B_red = LLLReduction(B)
  BIM = IntegerMatrix.from_matrix(B)
  B_red = LLL.reduction(BIM)
  b0 = list(B_red[0])
  g = b2f(b0,X)
  roots = np.roots(g[::-1]) # OVERFLOW HERE!!!
  # only consider roots in the integers
  return [round(x.real) for x in roots if evalP(g,round(x.real)) == 0]

def linspace (p_min, p_max, dist):
  points = [p_min+dist]
  last = p_min+dist
  yield last
  while points[-1]+dist < p_max:
    last = last + 2*dist
    yield last

# based on https://www.crypto.ruhr-uni-bochum.de/imperia/md/content/may/paper/lll.pdf
def computeParams (beta, delta, epsilon):
  m = ceil(pow(beta,2)/(delta*epsilon))
  t = floor(delta*m*(1/beta-1))
  X = Decimal(N)**Decimal(pow(beta,2)/delta - epsilon)
  X = int(X.to_integral_exact(rounding=ROUND_CEILING))//2
  return m, t, X

# X_max is the maximum value for the solution
# X_lim is the maximum value of X in the algorithm
def getRanges (X_max, X_lim):
  if 2*X_lim >= X_max:
    return X_max
  return X_lim

def getPLowerBits (N, p_min, X_max):
  if X_max >= pow(N,0.25):
    print("WARNING: X_max is too high")
  if p_min < sqrt(N):
    print("INFO: switching p and q")
    p_min = N//(p_min+X_max)
  beta = 0.5
  d = 1
  epsilon = 1/log(N)
  m, t, X_lim = computeParams(beta, d, epsilon)
  X = getRanges(X_max, X_lim)
  c = int(ceil(X_max/(2*X)))
  print(f"INFO: c = {hex(c)}, m = {m}, B of size {m*d+t}")
  i = 0
  for p_center in linspace(0, X_max, X):
    i += 1
    print(f"* Solving center {i}/{c}")
    f = [p_center+p_min,1]
    sols = coppersmith(f, N, X, m, t)
    print("  - Integer roots:")
    for s in sols:
      p_cand = f[0]+s
      print(f"       > {hex(s)}")
      if p_cand != 0 and N%(p_cand) == 0:
        print(f"INFO: found solution p = {hex(p_cand)}")
        return p_cand
  print("INFO: No solutions found!")
  return -1

# N = 0x00b8cb1cca99b6ac41876c18845732a5cbfc875df346ee9002ce608508b5fcf6b60a5ac7722a2d64ef74e1443a338e70a73e63a303f3ac9adf198595699f6e9f30c009d219c7d98c4ec84203610834029c79567efc08f66b4bc3f564bfb571546a06b7e48fb35bb9ccea9a2cd44349f829242078dfa64d525927bfd55d099c024f
# p_aprox = 0x00e700568ff506bd5892af92592125e06cbe9bd45dfeafe931a333c13463023d4f0000000000000000000000000000000000000000000000000000000000000000
# X_max   =                                                                  0x10000000000000000000000000000000000000000000000000000000000000000

N = 0x1533decf500f357917c4b2c359ad0ddc5d3271
p_aprox = 0x55cc5fbccd000000000
X_max   =          0x1000000000

p = getPLowerBits(N, p_aprox, X_max)

print(hex(p))


q = N//p
e = 0x10001
d = pow(e, -1, (p-1)*(q-1))

print(f"d = {hex(d)}")