factorization and multiplication in polynomial ring

asked 2022-02-17 15:07:07 +0100

rue82
125 ●2 ●4 ●13

updated 2022-02-18 11:12:56 +0100

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.

edit retag flag offensive close merge delete

Comments

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

Max Alekseyev ( 2022-02-17 19:04:19 +0100 )edit

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

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

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

Max Alekseyev ( 2022-02-17 22:00:33 +0100 )edit

You're right, that's not correct.

rue82 ( 2022-02-18 11:11:57 +0100 )edit

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

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

see more comments

2 Answers

Sort by » oldest newest most voted

answered 2022-02-17 20:05:51 +0100

Max Alekseyev

7303 ●10 ●63 ●178

updated 2022-02-22 15:39:43 +0100

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()

edit flag offensive delete link

Comments

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.

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

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

rue82 ( 2022-02-17 21:36:29 +0100 )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.

Max Alekseyev ( 2022-02-17 21:54:19 +0100 )edit

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

rue82 ( 2022-02-18 11:22:44 +0100 )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.

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

What exactly is not right? I believe my code is equivalent to yours, but just produces the result much faster.

Just in case, I've added an additional tracking for poles being cancelled to make sure that there is one pole in each X[i].

Max Alekseyev ( 2022-02-18 12:48:58 +0100 )edit

Thanks for helping me. Two comments: it seems that R = PolynomialRing(QQ, is better here, e.g. k=1 does not complain. When the code works, it is fast indeed. However, I don't think you can assume # assume that this happens only in the denominator of chi: try e.g. with K = 1+ q1 + q2 + q1*q2, k=4.

rue82 ( 2022-02-18 17:16:08 +0100 )edit

Well, the actual assumption needed here is that all zero factors in the numerator and denominator cancel each other. The result is correct as soon as this assumption is satisfied. My comment and the added poles tracking were not general enough though, I remove them. The code should work fine for K = 1+ q1 + q2 + q1*q2, k=4 (under the above assumption).

Max Alekseyev ( 2022-02-18 19:51:00 +0100 )edit

Yes, this assumption is correct. Let me run a few more checks. Note that k=1 is still broken as I pointed out.

rue82 ( 2022-02-19 21:54:30 +0100 )edit

I've corrected the issue for $k=1$.

Max Alekseyev ( 2022-02-22 02:30:37 +0100 )edit

Thanks a lot! it seems to work well.

rue82 ( 2022-02-22 10:39:29 +0100 )edit

see more comments

answered 2022-02-22 21:38:20 +0100

mwageringel
691 ●10 ●23

updated 2022-02-22 21:39:19 +0100

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))

edit flag offensive delete link

Comments

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

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

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

factorization and multiplication in polynomial ring

Comments

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

factorization and multiplication in polynomial ring edit

Comments

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

factorization and multiplication in polynomial ring