# 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 = q1 - is this intended?

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

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

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.

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

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.

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

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.

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?