Hi,
I implemented an extension to the powerseries command, designed to handle cases
that fall under the following theorem:
Rational Expansion Theorem for Distinct Roots
Graham, Knuth, Patashnik, "Concrete Math" 2nd ed, p 340
If R(z) = P(z)/Q(z), where Q(z) = q0*(1-r_1*z)*...*(1-r_l*z) and the
numbers (r_1 ... r_l) are distinct, and if P(z) is a polynomial of degree less
than l, then
[z^n]R(z) = a_1*r_1^n + ... + a_l*r_l^n, where a_k = -r_k*P(1/r_k)/Q'(1/r_k)
The code is here:
https://sourceforge.net/tracker/?func=detail&atid=304933&aid=1722156&group_id=4933
Here are some test cases:
/* two simple roots */
powerseries((1)/((1-2*x)*(1-3*x)), x, 0);
/* fibonacci */
powerseries((1)/(1-x-x^2), x, 0);
/* 1 1 2 2 4 4 8 8 */
powerseries((1+x)/(1-2*x^2), x, 0);
/* multiple root */
powerseries((1+x)/(1-x)^2, x, 0);
/* numerator higher order poly than denom */
powerseries((1+x^3)/(1-x-x^2), x, 0);
/* zero root in denom */
powerseries((1)/((1-2*x)*(x)), x, 0);
/* one simple and one repeated root in denom */
powerseries((1+x+x^2)/((1-2*x)*(1+x)^2), x, 0);
/* gcd of exps is two */
powerseries((1-x^2)/(1-4*x^2+x^4), x, 0);
Here is the behavior of the new version:
1
(%i2) powerseries(-------------------, x, 0)
(1 - 2 x) (1 - 3 x)
inf
====
\ i1 + 1 i1 i1
(%o2) > (3 - 2 2 ) x
/
====
i1 = 0
1
(%i3) powerseries(------------, x, 0)
2
- x - x + 1
inf
==== - i2 - 1 i2 - i2 - 1 i2
\ 2 (sqrt(5) - 1) 2 2 (sqrt(5) + 1) (- 2)
(%o3) - > (- --------------------------- - -------------------------------)
/ sqrt(5) sqrt(5)
====
i2 = 0
i2
x
x + 1
(%i4) powerseries(--------, x, 0)
2
1 - 2 x
inf 1 i3/2 i3 1 i3/2
==== (------- - 1) 2 (- 1) (- ------- - 1) 2
\ sqrt(2) sqrt(2) i3
(%o4) > (- --------------------------- - ---------------------) x
/ 2 2
====
i3 = 0
x + 1
(%i5) powerseries(--------, x, 0)
2
(1 - x)
inf
====
\ i4 i4
(%o5) > (2 (i4 + 1) x - x )
/
====
i4 = 0
3
x + 1
(%i6) powerseries(------------, x, 0)
2
- x - x + 1
inf
==== - i5 - 1 i5
\ 2 (sqrt(5) - 1) 2
(%o6) - 2 x > (- ---------------------------
/ sqrt(5)
====
i5 = 0
- i5 - 1 i5
2 (sqrt(5) + 1) (- 2) i5
- -------------------------------) x - x + 1
sqrt(5)
1
(%i7) powerseries(-----------, x, 0)
(1 - 2 x) x
inf
====
\ i6 i6
> 2 x
/
====
i6 = 0
(%o7) --------------
x
2
x + x + 1
(%i8) powerseries(------------------, x, 0)
2
(1 - 2 x) (x + 1)
inf
==== i7 i7 i7 i7 i7 i7
\ 7 2 x (i7 + 1) (- 1) x (- 1) x
(%o8) > (--------- + -------------------- - -----------)
/ 9 3 9
====
i7 = 0
2
1 - x
(%i9) powerseries(-------------, x, 0)
4 2
x - 4 x + 1
inf
==== - i8 - 1
\ (- sqrt(3) - 1) (sqrt(3) + 2)
(%o9) > (- -------------------------------------
/ 2 (sqrt(3) + 2) - 4
====
i8 = 0
- i8 - 1
(2 - sqrt(3)) (sqrt(3) - 1) 2 i8
- -----------------------------------) x
2 (2 - sqrt(3)) - 4
Here is the behavior of the old version:
1
(%i2) powerseries(-------------------, x, 0)
(1 - 2 x) (1 - 3 x)
inf inf
==== ====
\ i1 i1 + 1 i1 \ i1 i1 + 1 i1
(%o2) 2 > (- 2) (- 1) x - 3 > (- 3) (- 1) x
/ /
==== ====
i1 = 0 i1 = 0
1
(%i3) powerseries(------------, x, 0)
2
- x - x + 1
(%o3) Unable to expand
x + 1
(%i4) powerseries(--------, x, 0)
2
1 - 2 x
(%o4) Unable to expand
x + 1
(%i5) powerseries(--------, x, 0)
2
(1 - x)
inf inf
==== ====
\ i4 \ i4
(%o5) 2 > (i4 + 1) x - > x
/ /
==== ====
i4 = 0 i4 = 0
3
x + 1
(%i6) powerseries(------------, x, 0)
2
- x - x + 1
(%o6) Unable to expand
1
(%i7) powerseries(-----------, x, 0)
(1 - 2 x) x
inf
====
\ i6 i6 + 1 i6
> (- 2) (- 1) x
/
====
i6 = 0
(%o7) - ------------------------------
x
2
x + x + 1
(%i8) powerseries(------------------, x, 0)
2
(1 - 2 x) (x + 1)
inf inf
==== ====
\ i7 i7 + 1 i7 \ i7 i7
7 > (- 2) (- 1) x > (i7 + 1) (- 1) x
/ /
==== ====
i7 = 0 i7 = 0
(%o8) - -------------------------------- + ---------------------------
9 3
inf
====
\ i7 i7
> (- 1) x
/
====
i7 = 0
- ------------------
9
2
1 - x
(%i9) powerseries(-------------, x, 0)
4 2
x - 4 x + 1
(%o9) Unable to expand
Could a developer please take a look at this proposed change?