# extract coefficients from products and sums of descending products

I need to extract coefficients of $n$ from some powers, products and sums of descending products and their reciprocals. The first 5 terms (in $x$) of the descending product $(x)_y$ are $$x^y \left (1 - \frac{y(y-1)}{2x} + \frac{y(y-1)(y-2)(3y-1)}{24x^2} - \frac{y^2(y - 1)^2(y - 2)(y - 3)}{48x^3} + \frac{y(y-1)(y-2)(y-3)(y-4)(15y^3-30y^2+5y+2)}{5760x^4}\right)$$ and the first 5 terms of its reciprocal are $$x^{-y} \left(1 + \frac{y(y-1)}{2x} + \frac{y(y-1)(y+1)(3y-2)}{24x^2} + \frac{y^2(y-1)^2(y+1)(y+2)}{48x^3} + \frac{y(y-1)(y+1)(y+2)(y+3)(15y^3-15y^2-10y+8)}{5760x^4}\right)$$

I need to extract the coefficients of various powers of $n$ in expressions like

$$\frac{(l(v-1)-1)(v-1)l(l)_t}{(n-1)_{t+1}} - \left(\frac{(l(v-1)-1)(v-1)l(l)_t}{(n-1)_{t+1}}\right)^2 + \frac{n(v-1)l(l)_{2t}(l(v-1)-1)^2(l(v-1)-2)}{(n)_{2t+3}}$$ where $l = n/v$ and $n$, $t$ and $v$ are symbolic variables.

To encapsulate these descending products and their powers I wrote

def desc_prod(x,y):
t0 = 1
t1 = - y*(y-1)/(2*x)
t2 = (y*(y-1)*(y-2)*(3*y-1))/(24*x^2)
t3 = -(y^2*(y-1)^2*(y-2)*(y-3))/(48*x^3)
t4 = y*(y-1)*(y-2)*(y-3)*(y-4)*(15*y^3-30*y^2 + 5*y + 2)/(5760*x^4)
return (x^y)*(t0 + t1 + t2 + t3 + t4 )

def desc_prod_recip(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 (x^(-y))*(t0 + t1 + t2 + t3 + t4)

def desc_prod_power(x,y,e):
t0 = 1
t1 = - y*(y-1)/(2*x)
t2 = (y*(y-1)*(y-2)*(3*y-1))/(24*x^2)
t3 = -(y^2*(y-1)^2*(y-2)*(y-3))/(48*x^3)
t4 = y*(y-1)*(y-2)*(y-3)*(y-4)*(15*y^3-30*y^2 + 5*y + 2)/(5760*x^4)
return (x^(y*e))*(1 + multinomial(e-1,1)*t1 + multinomial(e-2,2)*t1^2 + multinomial(e-1,1)*t2  + multinomial(e-3,3)*t1^3 + multinomial(e-2,1,1)*t1*t2 + multinomial(e-4,4)*t1^4 + multinomial(e-2,2)*t2^2 + multinomial(e-3,2,1)*t1^2*t2 + multinomial(e-2,1,1)*t1*t3 + multinomial(e-1,1)*t4 )

def desc_prod_recip_power(x,y,e):
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 (x^(-y*e))*(1 + multinomial(e-1,1)*t1 + multinomial(e-2,2)*t1^2 + multinomial(e-1,1)*t2  + multinomial(e-3,3)*t1^3 + multinomial(e-2,1,1)*t1*t2 + multinomial(e-4,4)*t1^4 + multinomial(e-2,2)*t2^2 + multinomial(e-3,2,1)*t1^2*t2 + multinomial(e-2,1,1)*t1*t3 + multinomial(e-1,1)*t4 )


but sagemath does not correctly compute the coefficients in $n$. In an attempt to figure this out I asked a question here about computing coefficients of $n$ correctly and you can see there a very simple instance of sagemath failing to compute the coefficient. In my question I tried to give the simplest situation, where I do not have powers, products nor sums of my descnding products and their reciprocals, but just enough to showcase the behaviour of sage that is getting me stuck. In his answer @dan_fulea suggests that I need to approach my problem differently. Can anyone please give me some advice about the best way to compute the coefficients of any given power of $n$ in expressions like mine.

edit retag close merge delete

Sort by » oldest newest most voted

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=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.

more

1

Thank you @rburing this was very helpful. I have been able to move forward with my computations.

( 2020-06-01 23:43:50 +0100 )edit