Ask Your Question

Revision history [back]

If you don't mind I will just focus on the example of desc_prod_recip as in your previous question. It should be possible to generalize this approach. You want to work with a symbolic expression like

sage: n,t = var('n t')
sage: E1 = desc_prod_recip(n,t).expand().simplify(); E1
1/384*n^(-t - 4)*t^8 + 1/96*n^(-t - 4)*t^7 + 1/48*n^(-t - 3)*t^6 - 1/576*n^(-t - 4)*t^6 + 1/48*n^(-t - 3)*t^5 - 1/30*n^(-t - 4)*t^5 + 1/8*n^(-t - 2)*t^4 - 1/16*n^(-t - 3)*t^4 - 5/1152*n^(-t - 4)*t^4 - 1/12*n^(-t - 2)*t^3 - 1/48*n^(-t - 3)*t^3 + 1/32*n^(-t - 4)*t^3 + 1/2*n^(-t - 1)*t^2 - 1/8*n^(-t - 2)*t^2 + 1/24*n^(-t - 3)*t^2 + 1/288*n^(-t - 4)*t^2 - 1/2*n^(-t - 1)*t + 1/12*n^(-t - 2)*t - 1/120*n^(-t - 4)*t + 1/n^t

This is almost a rational function in n (which would be easy to manipulate), except that you have this pesky symbolic exponent t. Suggestion: introduce a new variable z which represents n^t, so that the expression becomes a rational function in n and z.

def desc_prod_recip(x,y,z):
    t0 = 1
    t1 = y*(y-1)/(2*x)
    t2 = (y*(y+1)*(y-1)*(3*y-2))/(24*x^2)
    t3 = (y^2*(y-1)^2*(y+1)*(y+2))/(48*x^3)
    t4 = y*(y-1)*(y+1)*(y+2)*(y+3)*(15*y^3-15*y^2-10*y+8)/(5760*x^4)
    return (z^(-1))*(t0 + t1 + t2 + t3 + t4)

To simplify our job, use the fraction field of a polynomial ring rather than the symbolic ring:

sage: R.<n,t,z> = PolynomialRing(QQ)
sage: E1 = desc_prod_recip(n,t,z); E1
(34560*t^8 + 276480*n*t^6 + 138240*t^7 + 1658880*n^2*t^4 + 276480*n*t^5 - 23040*t^6 + 6635520*n^3*t^2 - 1105920*n^2*t^3 - 829440*n*t^4 - 442368*t^5 + 13271040*n^4 - 6635520*n^3*t - 1658880*n^2*t^2 - 276480*n*t^3 - 57600*t^4 + 1105920*n^2*t + 552960*n*t^2 + 414720*t^3 + 46080*t^2 - 110592*t)/(13271040*n^4*z)

To get the (actual) exponents of n, get the exponents of monomials in n and z and weigh them appropriately:

denom = E1.denominator()
denom_coeff = denom.coefficients()[0]
denom_monom = denom.monomials()[0]
for (c,m) in E1.numerator():
    n_degree = m.degree(n) - denom_monom.degree(n) + t*(m.degree(z) - denom_monom.degree(z))
    t_degree = m.degree(t) - denom_monom.degree(t)
    print(c/denom_coeff, SR(t)^t_degree * SR(n)^SR(n_degree))

Output:

1/384 n^(-t - 4)*t^8
1/48 n^(-t - 3)*t^6
1/96 n^(-t - 4)*t^7
1/8 n^(-t - 2)*t^4
1/48 n^(-t - 3)*t^5
-1/576 n^(-t - 4)*t^6
1/2 n^(-t - 1)*t^2
-1/12 n^(-t - 2)*t^3
-1/16 n^(-t - 3)*t^4
-1/30 n^(-t - 4)*t^5
1 1/(n^t)
-1/2 n^(-t - 1)*t
-1/8 n^(-t - 2)*t^2
-1/48 n^(-t - 3)*t^3
-5/1152 n^(-t - 4)*t^4
1/12 n^(-t - 2)*t
1/24 n^(-t - 3)*t^2
1/32 n^(-t - 4)*t^3
1/288 n^(-t - 4)*t^2
-1/120 n^(-t - 4)*t

You can also collect terms according to n_degree, etc.

If you don't mind I will just focus on the example of desc_prod_recip as in your previous question. It should be possible to generalize this approach. You want to work with a symbolic expression like

sage: n,t = var('n t')
sage: E1 = desc_prod_recip(n,t).expand().simplify(); E1
1/384*n^(-t - 4)*t^8 + 1/96*n^(-t - 4)*t^7 + 1/48*n^(-t - 3)*t^6 - 1/576*n^(-t - 4)*t^6 + 1/48*n^(-t - 3)*t^5 - 1/30*n^(-t - 4)*t^5 + 1/8*n^(-t - 2)*t^4 - 1/16*n^(-t - 3)*t^4 - 5/1152*n^(-t - 4)*t^4 - 1/12*n^(-t - 2)*t^3 - 1/48*n^(-t - 3)*t^3 + 1/32*n^(-t - 4)*t^3 + 1/2*n^(-t - 1)*t^2 - 1/8*n^(-t - 2)*t^2 + 1/24*n^(-t - 3)*t^2 + 1/288*n^(-t - 4)*t^2 - 1/2*n^(-t - 1)*t + 1/12*n^(-t - 2)*t - 1/120*n^(-t - 4)*t + 1/n^t

This is almost a rational function in n (which would be easy to manipulate), except that you have this pesky symbolic exponent t. Suggestion: introduce a new variable z which represents n^t, so that the expression becomes a rational function in n and z.

def desc_prod_recip(x,y,z):
desc_prod_recip(x,y,z=None):
    if z is None: z = x^y
    t0 = 1
    t1 = y*(y-1)/(2*x)
    t2 = (y*(y+1)*(y-1)*(3*y-2))/(24*x^2)
    t3 = (y^2*(y-1)^2*(y+1)*(y+2))/(48*x^3)
    t4 = y*(y-1)*(y+1)*(y+2)*(y+3)*(15*y^3-15*y^2-10*y+8)/(5760*x^4)
    return (z^(-1))*(t0 + t1 + t2 + t3 + t4)

To simplify our job, use the fraction field of a polynomial ring rather than the symbolic ring:

sage: R.<n,t,z> = PolynomialRing(QQ)
sage: E1 = desc_prod_recip(n,t,z); E1
(34560*t^8 + 276480*n*t^6 + 138240*t^7 + 1658880*n^2*t^4 + 276480*n*t^5 - 23040*t^6 + 6635520*n^3*t^2 - 1105920*n^2*t^3 - 829440*n*t^4 - 442368*t^5 + 13271040*n^4 - 6635520*n^3*t - 1658880*n^2*t^2 - 276480*n*t^3 - 57600*t^4 + 1105920*n^2*t + 552960*n*t^2 + 414720*t^3 + 46080*t^2 - 110592*t)/(13271040*n^4*z)

To get the (actual) exponents of n, get the exponents of monomials in n and z and weigh them appropriately:

denom = E1.denominator()
denom_coeff = denom.coefficients()[0]
denom_monom = denom.monomials()[0]
for (c,m) in E1.numerator():
    n_degree = m.degree(n) - denom_monom.degree(n) + t*(m.degree(z) - denom_monom.degree(z))
    t_degree = m.degree(t) - denom_monom.degree(t)
    print(c/denom_coeff, SR(t)^t_degree * SR(n)^SR(n_degree))

Output:

1/384 n^(-t - 4)*t^8
1/48 n^(-t - 3)*t^6
1/96 n^(-t - 4)*t^7
1/8 n^(-t - 2)*t^4
1/48 n^(-t - 3)*t^5
-1/576 n^(-t - 4)*t^6
1/2 n^(-t - 1)*t^2
-1/12 n^(-t - 2)*t^3
-1/16 n^(-t - 3)*t^4
-1/30 n^(-t - 4)*t^5
1 1/(n^t)
-1/2 n^(-t - 1)*t
-1/8 n^(-t - 2)*t^2
-1/48 n^(-t - 3)*t^3
-5/1152 n^(-t - 4)*t^4
1/12 n^(-t - 2)*t
1/24 n^(-t - 3)*t^2
1/32 n^(-t - 4)*t^3
1/288 n^(-t - 4)*t^2
-1/120 n^(-t - 4)*t

You can also collect terms according to n_degree, etc.