fast factorization of ratios of polynomials

I consider polynomials of the form $P=\prod_i (1-1/y_i)$, where each $y_i$ is a Laurent monomial with unit coefficient in variables $x_1,\ldots,x_k, q_1,q_2,q_3,m$.

I'm interested in taking ratios of the form $P_1 / (P_2 \prod_{n=1}^k x_n)$ of such polynomials, and would like to know what the fastest way is to cancel common factors using sagemath.

The final goal is to take residues at simple poles, so factorizations will happen.

Let me write a full code example:

def br(x):
return 1-1/x

q1,q2,q3,m = var('q1,q2,q3,m')

def mex(q1,q2,q3,m):
q4 = var('q4')
N.<q1,q2,q3,q4> = PolynomialRing(ZZ, 4, order='neglex')
k = 3
K = 1 + q1 + q2
rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
X = [var("x%d" % i) for i in range(k)]
chim = prod([ br(X[j]/m) for j in range(k)])
chix = prod([ br(X[j]) for j in range(k)])
chiup = prod([ prod([ br(q1*q2*X[i]/X[j])*br(q1*q3*X[i]/X[j])*br(q2*q3*X[i]/X[j])*(br(X[i]/X[j]))^2 for i in range(k) if i > j]) for j in range(k)])
chiups = prod([ prod([ br(X[i]/(X[j]*q1*q2))*br(X[i]/(X[j]*q1*q3))*br(X[i]/(X[j]*q2*q3)) for i in range(k) if i > j]) for j in range(k)])
chido = prod([ prod([ br(q1*X[i]/X[j])*br(q2*X[i]/X[j])*br(q3*X[i]/X[j])*br(q1*q2*q3*X[i]/X[j]) for i in range(k) if i > j]) for j in range (k)])
chidos = prod([ prod([ br(X[i]/(X[j]*q1))*br(X[i]/(X[j]*q2))*br(X[i]/(X[j]*q3))*br(X[i]/(X[j]*q1*q2*q3)) for i in range(k) if i > j]) for j in range (k)])
dx = prod([ X[j] for j in range(k)])
chinum = (chim*chiup*chiups)
chiden = (chix*chido*chidos*dx)
chi = chinum/chiden
# return chi.factor()
for xi,rhoi in zip(X,rho):
chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})
return chi


This currently works in the symbolic ring, but gets slow as soon as k,K grow a bit. I was hoping that using singular (Laurent polyn ring, e.g.) would make it faster, but cannot make it such that everything happens there yet.

edit retag close merge delete

It's unclear what is your problem. Your polynomials are already given in the factored form - why do you want to factor them?

( 2022-02-05 21:51:40 +0200 )edit

Well, numerator and denominator are factored separately, but I want to cancel common factors when taking the ratio, namely I want to factor the ratio. Is that more clear?

( 2022-02-06 10:07:52 +0200 )edit

Could you please provide your full code ? In particular, how are q1,q2,q3,x0,x1,x2,x3 defined ?

( 2022-02-06 12:20:12 +0200 )edit

@tmonteil I added a code example, does this clarify my question?

( 2022-02-07 10:19:34 +0200 )edit

Couple of issues in your code: 1 - q1,q2,q3,m = var('q1,q2,q3,m') and q4 = var('q4') are redundant and can be safely removed 2 - defining X as symbolic variables slow down things, you'd better define them as generators of Laurent polynomial ring.

( 2022-02-07 13:41:02 +0200 )edit

Sort by » oldest newest most voted

Represent $y_i$ as a vector of size $k+4$ with components being the degrees of $x_1, \dots, x_k, q_1, q_2, q_3, m$. Then create a dictionary D with such vectors as keys and integer values (counters with zero initial values). With each $y_i$ coming from $P_1$ increase the corresponding value in D, while with every $y_i$ coming from $P_2$ decrease the corresponding value in D. The ratio $P_1/P_2$ can be easily reconstructed from the resulting D and it won't contain any common factors between the numerator and denominator.

ADDED. Here a rough translation of your code to these settings. In fact, it shows that no cancellations happen in your polynomials as they have no common factors, while the slow down is caused by size of the resulting object.

def mul_br(D,x,inc=1):
for t in x.dict():
D.setdefault(t,0)
D[t] += inc
if D[t]==0:
del D[t]

def mex():
k = 3
L = LaurentPolynomialRing(ZZ, 5+k, 'mqx')
m,q1,q2,q3,q4 = L.gens()[:5]
X = L.gens()[5:]

D = dict()

#chim = prod([ br(X[j]/m) for j in range(k)])
for j in range(k):
mul_br(D,X[j]/m)
#chix = prod([ br(X[j]) for j in range(k)])
for j in range(k):
mul_br(D,X[j],-1)
#chiup = prod([ prod([ br(q1*q2*X[i]/X[j])*br(q1*q3*X[i]/X[j])*br(q2*q3*X[i]/X[j])*(br(X[i]/X[j]))^2 for i in range(k) if i > j]) for j in range(k)])
for j in range(k):
for i in range(k):
if i > j:
mul_br(D,q1*q2*X[i]/X[j])
mul_br(D,q1*q3*X[i]/X[j])
mul_br(D,q2*q3*X[i]/X[j])
mul_br(D,X[i]/X[j],2)
#chiups = prod([ prod([ br(X[i]/(X[j]*q1*q2))*br(X[i]/(X[j]*q1*q3))*br(X[i]/(X[j]*q2*q3)) for i in range(k) if i > j]) for j in range(k)])
for j in range(k):
for i in range(k):
if i > j:
mul_br(D,X[i]/(X[j]*q1*q2))
mul_br(D,X[i]/(X[j]*q1*q3))
mul_br(D,X[i]/(X[j]*q2*q3))
#chido = prod([ prod([ br(q1*X[i]/X[j])*br(q2*X[i]/X[j])*br(q3*X[i]/X[j])*br(q1*q2*q3*X[i]/X[j]) for i in range(k) if i > j]) for j in range (k)])
for j in range (k):
for i in range(k):
if i > j:
mul_br(D,q1*X[i]/X[j],-1)
mul_br(D,q2*X[i]/X[j],-1)
mul_br(D,q3*X[i]/X[j],-1)
mul_br(D,q1*q2*q3*X[i]/X[j],-1)
#chidos = prod([ prod([ br(X[i]/(X[j]*q1))*br(X[i]/(X[j]*q2))*br(X[i]/(X[j]*q3))*br(X[i]/(X[j]*q1*q2*q3)) for i in range(k) if i > j]) for j in range (k)])
for j in range (k):
for i in range(k):
if i > j:
mul_br(D,X[i]/(X[j]*q1),-1)
mul_br(D,X[i]/(X[j]*q2),-1)
mul_br(D,X[i]/(X[j]*q3),-1)
mul_br(D,X[i]/(X[j]*q1*q2*q3),-1)
dx = prod(X[j] for j in range(k))
#chinum = (chim*chiup*chiups)
#chiden = (chix*chido*chidos*dx)
#chi = chinum/chiden

return D


Running mex() returns an answer in form of a dictionary with 51 elements. If you want an actual polynomial from it, you can get it by running

prod( (1-1/prod(g^e for g,e in zip(L.gens(),d)))^p for d,p in D.items() ) / dx


although it may be time consuming.

more

I added a full explanation, sorry for being not complete. I'm not sure using a dictionary is the fastest way here?

( 2022-02-10 14:25:51 +0200 )edit

I'm not sure if it's the fastest but it does the job quite efficiently. Also, it shows that the issue is not in canceling common factors (as there are none) and/or factorization, but in the size of the resulting object.

( 2022-02-10 15:02:22 +0200 )edit

But see the edited question, as there will be cancellations.

( 2022-02-10 16:23:39 +0200 )edit

do you think using the dictionary could work for the above?

( 2022-02-10 16:24:13 +0200 )edit

How exactly did you change the question? I have already demonstrated how the dictionary works. So, I don't quite follow your concern.

( 2022-02-10 20:33:59 +0200 )edit