# factorization and multiplication in polynomial ring

I'm looking for an efficient (i.e. fast enough for $k=4$ and higher) way in Sage to build large rational functions (ratio of polynomials with integer coefficients), and then take residues at simple poles (which I know by construction), by multiplication and taking limit (via subs).

I can do that in the symbolic ring, which I believe uses maxima or sympy, but at some point it becomes slow in factorization (say compared to mathematica), so I tried to work with polynomial rings, hoping that singular would be faster.

However, at present I have issues:

1. since by default it expands whatever factorized polynomial I define, then it has a hard time factorizing it back (see how long it takes even just to return $\chi$).

2. Moreover, I cannot operate with factorization objects, e.g. I need to multiply them, but it seems this is not allowed.

This is my current attempt:

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

def measNew():
k = 3
R = PolynomialRing(ZZ, ['x%d' %p for p in range(k)]+['q1', 'q2', 'q3', 'q4', 'm'], order='invlex')
R.inject_variables()
K = 1 + q2 + q2^2 + q3 + q4
X = R.gens()[:k]
rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
rho.reverse()
chix = prod([ br(X[j]) for j in range(k)])
chim = prod([ br(X[j]/m) 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 # just to check the intermediate step
for xi,rhoi in zip(X,rho):
chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})
return chi.factor()


Any help/suggestions?

Thanks.

edit retag close merge delete

Shouldn't it be X = R.gens()[:k]? At the moment we have, e.g. X[3] = q1 - is this intended?

( 2022-02-17 19:04:19 +0200 )edit

You're right. I never use X[p] with p \geq k, I just didn't know how to code it correctly. Thanks.

( 2022-02-17 21:27:35 +0200 )edit

If you have X = R.gens(), then zip(X,rho) has more terms than needed, including substitution of q1 with something.

( 2022-02-17 22:00:33 +0200 )edit

You're right, that's not correct.

( 2022-02-18 11:11:57 +0200 )edit

I've also reversed rho, to get correct ordering.

( 2022-02-18 11:13:35 +0200 )edit

Sort by » oldest newest most voted

You can try to perform substitution and cancellation of simple poles directly in your br function, which significantly speed up things:

mysubs = []

def br(x):
r = 1-1/x
for t,s in mysubs:
m = 0
while r.subs({t:s})==0:
r = r.derivative(t)
m += 1
r = r.subs({t:s}) / factorial(m)
return r

def measNew():
k = 3
R = PolynomialRing(ZZ, ['x%d' %p for p in range(k)]+['q1', 'q2', 'q3', 'q4', 'm'], order='invlex')
F = R.fraction_field()
R.inject_variables()
K = 1 + q2 + q2^2 + q3 + q4
X = R.gens()[:k]
rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
rho.reverse()

global mysubs
mysubs = list( zip(X,rho) )

chix = prod( br(X[j]) for j in range(k) )
chim = prod( br(X[j]/m) 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(j+1,k) ) 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(j+1,k) ) 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(j+1,k) ) 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(j+1,k) ) for j in range (k) )
dx = prod( X[j] for j in range(k) ).subs(dict(mysubs))
chinum = chim * chiup * chiups
chiden = chix * chido * chidos * dx
chi = F(chinum) / F(chiden)

return chi.factor()

more

Thanks, let me study your suggestions. Is there a reason why you changed the order from invlex? the actual ordering should be exactly the opposite of invlex btw.

( 2022-02-17 21:35:45 +0200 )edit

in case I did not make it clear: the order in which residues are taken is crucial

( 2022-02-17 21:36:29 +0200 )edit

The result should not depend on the term order (I've changed it back to invlex anyway). As for the substitution order, I've changed mysubs into a list to enforce the fixed order.

Btw, please notice that not all poles are simple. E.g., consider the term br(q2*X[i]/X[j]) of chido, which for j=0 and i=1 contributes (x1*q2 - x0)/(x1*q2) to the denominator of chi. This term becomes zero under simultaneous substitution {x0: q2^2, x1: q2}. Multiplying numerator by x0-q2^2 followed by subs({x0: q2^2}) and then multiplying by x1-q2 followed by subs({x1: q2}) may not treat it correctly.

( 2022-02-17 21:54:19 +0200 )edit

Simultaneous substitution is not allowed here: we always compute an iterated residue, for which order matters.

( 2022-02-18 11:22:44 +0200 )edit

Sorry, but this does not look right: cancellations between numerator and denominator of chi will happen along the way (i.e. after some substitutions in my for cycle).

What is really needed is a way to 1. not expand a polynomial that was given in factorized form, and 2. multiply factorized() objects: I'm really surprised this is not doable, and I was hoping it is just due to my ignorance of sage/singular.

( 2022-02-18 11:29:55 +0200 )edit

It might be better to do this computation in the symbolic ring instead in order to avoid expansions, i.e.

# from the question (for completeness)
def br(x):
return 1-1/x

k = 3
gens = var(['x%d' %p for p in range(k)]+['q1', 'q2', 'q3', 'q4', 'm'])
K = 1 + q2 + q2^2 + q3 + q4
X = gens[:k]
rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.operands()]
rho.reverse()
chix = prod([ br(X[j]) for j in range(k)])
chim = prod([ br(X[j]/m) 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)


Then carefully handle the products without calling factor on the entire fraction, but only on the individual factors that we already know.

for xi,rhoi in zip(X,rho):
D = {xi: rhoi}
chi = (chi*(xi-rhoi))
g = prod([f.factor() for f in chi.operands()])  # this seems to have cancelled all fractions, but maybe add assertions for this
dens = [(1/gi).subs(D) for gi in g.operands() if gi.numerator().is_constant()]
nums = [gi.subs(D) for gi in g.operands() if not gi.numerator().is_constant()]
chi = prod(nums) / prod(dens)

chi = prod([f.factor() for f in chi.operands()])


Here, g as well as the last line seem to end up being fractions of factors of actual polynomials (but I am not sure if this will work in every case). Result:

-(q1^2*q2*q3^3 - 1)*(q1*q2^2*q3^3 - 1)*(q1^2*q2^2*q3 - 1)*(m*q1*q2*q3 - 1)*(q1*q2*q3 - 1)*(q1*q2 - q3)*(m - q3)*(m - 1)*(q1 - 1)^2*(q2 - 1)^2*(q3 - 1)/((q1^2*q2^2*q3^3 - 1)*(q1*q2*q3^3 - 1)*(q1^2*q2*q3 - 1)*(q1*q2^2*q3 - 1)*(q1*q2 - 1)^2*(q1*q3 - 1)*(q2*q3 - 1)*(q1 - q3)*(q2 - q3))

more

Thanks, although I thought symbolic computations would be slower here, right?

( 2022-02-25 16:09:08 +0200 )edit