Ask Your Question
1

fast factorization of ratios of polynomials

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

rue82 gravatar image

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

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.

edit retag flag offensive close merge delete

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 gravatar imageMax Alekseyev ( 2022-02-05 21:51:40 +0100 )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 gravatar imagerue82 ( 2022-02-06 10:07:52 +0100 )edit

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

tmonteil gravatar imagetmonteil ( 2022-02-06 12:20:12 +0100 )edit

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

rue82 gravatar imagerue82 ( 2022-02-07 10:19:34 +0100 )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 gravatar imageMax Alekseyev ( 2022-02-07 13:41:02 +0100 )edit

1 Answer

Sort by ยป oldest newest most voted
2

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

Max Alekseyev gravatar image

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

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 more

Comments

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

rue82 gravatar imagerue82 ( 2022-02-10 14:25:51 +0100 )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 gravatar imageMax Alekseyev ( 2022-02-10 15:02:22 +0100 )edit

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

rue82 gravatar imagerue82 ( 2022-02-10 16:23:39 +0100 )edit

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

rue82 gravatar imagerue82 ( 2022-02-10 16:24:13 +0100 )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 gravatar imageMax Alekseyev ( 2022-02-10 20:33:59 +0100 )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-05 19:02:55 +0100

Seen: 388 times

Last updated: Feb 10 '22