Revision history [back]

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()
    rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
    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
    for xi,rhoi in zip(X,rho):
        chi = (chi*(xi-rhoi)).factor().subs({xi: rhoi})
    return chi.factor()

measNew()

Any help/suggestions?

Thanks.

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()
    rho = [p.substitute({q4:(q1*q2*q3)^-1}) for p in K.monomials()]
    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
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()

measNew()

Any help/suggestions?

Thanks.

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