fast factorization of ratios of polynomials

asked 2022-02-05 19:02:55 +0200

rue82
113 ●4 ●13

updated 2022-02-10 14:24:55 +0200

I consider polynomials of the form $P=\prod_i (1-1/y_i)$, where each $y_i$ is a Laurent monomial with unit coefficient in variables $x_1,\ldots,x_k, q_1,q_2,q_3,m$.

I'm interested in taking ratios of the form $P_1 / (P_2 \prod_{n=1}^k x_n)$ of such polynomials, and would like to know what the fastest way is to cancel common factors using sagemath.

The final goal is to take residues at simple poles, so factorizations will happen.

Let me write a full code example:

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

q1,q2,q3,m = var('q1,q2,q3,m')

def mex(q1,q2,q3,m):
    q4 = var('q4')
    N.<q1,q2,q3,q4> = PolynomialRing(ZZ, 4, order='neglex')
    k = 3
    K = 1 + q1 + q2
    rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
    X = [var("x%d" % i) for i in range(k)]
    chim = prod([ br(X[j]/m) for j in range(k)])
    chix = prod([ br(X[j]) 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.factor()
    for xi,rhoi in zip(X,rho):
         chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})
    return chi

This currently works in the symbolic ring, but gets slow as soon as k,K grow a bit. I was hoping that using singular (Laurent polyn ring, e.g.) would make it faster, but cannot make it such that everything happens there yet.

Comments

It's unclear what is your problem. Your polynomials are already given in the factored form - why do you want to factor them?

Max Alekseyev ( 2022-02-05 21:51:40 +0200 )edit

Well, numerator and denominator are factored separately, but I want to cancel common factors when taking the ratio, namely I want to factor the ratio. Is that more clear?

rue82 ( 2022-02-06 10:07:52 +0200 )edit

Could you please provide your full code ? In particular, how are q1,q2,q3,x0,x1,x2,x3 defined ?

tmonteil ( 2022-02-06 12:20:12 +0200 )edit

@tmonteil I added a code example, does this clarify my question?

rue82 ( 2022-02-07 10:19:34 +0200 )edit

Couple of issues in your code: 1 - q1,q2,q3,m = var('q1,q2,q3,m') and q4 = var('q4') are redundant and can be safely removed 2 - defining X as symbolic variables slow down things, you'd better define them as generators of Laurent polynomial ring.

Max Alekseyev ( 2022-02-07 13:41:02 +0200 )edit

see more comments

answered 2022-02-06 15:42:56 +0200

Max Alekseyev

5748 ●5 ●37 ●123

updated 2022-02-07 17:05:42 +0200

Represent $y_i$ as a vector of size $k+4$ with components being the degrees of $x_1, \dots, x_k, q_1, q_2, q_3, m$. Then create a dictionary D with such vectors as keys and integer values (counters with zero initial values). With each $y_i$ coming from $P_1$ increase the corresponding value in D, while with every $y_i$ coming from $P_2$ decrease the corresponding value in D. The ratio $P_1/P_2$ can be easily reconstructed from the resulting D and it won't contain any common factors between the numerator and denominator.

ADDED. Here a rough translation of your code to these settings. In fact, it shows that no cancellations happen in your polynomials as they have no common factors, while the slow down is caused by size of the resulting object.

def mul_br(D,x,inc=1):
    for t in x.dict():
        D.setdefault(t,0)
        D[t] += inc
        if D[t]==0:
            del D[t]

def mex():
    k = 3
    L = LaurentPolynomialRing(ZZ, 5+k, 'mqx')
    m,q1,q2,q3,q4 = L.gens()[:5]
    X = L.gens()[5:]

    D = dict()

    #chim = prod([ br(X[j]/m) for j in range(k)])
    for j in range(k):
        mul_br(D,X[j]/m)
    #chix = prod([ br(X[j]) for j in range(k)])
    for j in range(k):
        mul_br(D,X[j],-1)
    #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)])
    for j in range(k):
        for i in range(k):
            if i > j:
                mul_br(D,q1*q2*X[i]/X[j])
                mul_br(D,q1*q3*X[i]/X[j])
                mul_br(D,q2*q3*X[i]/X[j])
                mul_br(D,X[i]/X[j],2)
    #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)])
    for j in range(k):
        for i in range(k):
            if i > j:
                mul_br(D,X[i]/(X[j]*q1*q2))
                mul_br(D,X[i]/(X[j]*q1*q3))
                mul_br(D,X[i]/(X[j]*q2*q3))
    #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)])
    for j in range (k):
        for i in range(k):
            if i > j:
                mul_br(D,q1*X[i]/X[j],-1)
                mul_br(D,q2*X[i]/X[j],-1)
                mul_br(D,q3*X[i]/X[j],-1)
                mul_br(D,q1*q2*q3*X[i]/X[j],-1)
    #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)])
    for j in range (k):
        for i in range(k):
            if i > j:
                mul_br(D,X[i]/(X[j]*q1),-1)
                mul_br(D,X[i]/(X[j]*q2),-1)
                mul_br(D,X[i]/(X[j]*q3),-1)
                mul_br(D,X[i]/(X[j]*q1*q2*q3),-1)
    dx = prod(X[j] for j in range(k))
    #chinum = (chim*chiup*chiups)
    #chiden = (chix*chido*chidos*dx)
    #chi = chinum/chiden

    return D

Running mex() returns an answer in form of a dictionary with 51 elements. If you want an actual polynomial from it, you can get it by running

prod( (1-1/prod(g^e for g,e in zip(L.gens(),d)))^p for d,p in D.items() ) / dx

although it may be time consuming.

edit flag offensive delete link

Comments

I added a full explanation, sorry for being not complete. I'm not sure using a dictionary is the fastest way here?

rue82 ( 2022-02-10 14:25:51 +0200 )edit

I'm not sure if it's the fastest but it does the job quite efficiently. Also, it shows that the issue is not in canceling common factors (as there are none) and/or factorization, but in the size of the resulting object.

Max Alekseyev ( 2022-02-10 15:02:22 +0200 )edit

But see the edited question, as there will be cancellations.

rue82 ( 2022-02-10 16:23:39 +0200 )edit

do you think using the dictionary could work for the above?

rue82 ( 2022-02-10 16:24:13 +0200 )edit

How exactly did you change the question? I have already demonstrated how the dictionary works. So, I don't quite follow your concern.

Max Alekseyev ( 2022-02-10 20:33:59 +0200 )edit

see more comments

Your Answer

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

fast factorization of ratios of polynomials

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

fast factorization of ratios of polynomials edit

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

fast factorization of ratios of polynomials