Ask Your Question
1

factorization and multiplication in polynomial ring

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

rue82 gravatar image

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

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 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 gravatar imageMax Alekseyev ( 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.

rue82 gravatar imagerue82 ( 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.

Max Alekseyev gravatar imageMax Alekseyev ( 2022-02-17 22:00:33 +0200 )edit

You're right, that's not correct.

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

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

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

2 Answers

Sort by ยป oldest newest most voted
1

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

Max Alekseyev gravatar image

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

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 more

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 gravatar imagerue82 ( 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

rue82 gravatar imagerue82 ( 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.

Max Alekseyev gravatar imageMax Alekseyev ( 2022-02-17 21:54:19 +0200 )edit

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

rue82 gravatar imagerue82 ( 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.

rue82 gravatar imagerue82 ( 2022-02-18 11:29:55 +0200 )edit
1

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

mwageringel gravatar image

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

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 more

Comments

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

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

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2022-02-17 15:07:07 +0200

Seen: 96 times

Last updated: Feb 22