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:
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$).
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.
Shouldn't it be
X = R.gens()[:k]
? At the moment we have, e.g.X[3] = 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()
, thenzip(X,rho)
has more terms than needed, including substitution ofq1
with something.You're right, that's not correct.
I've also reversed rho, to get correct ordering.