##
##    Title:	fullparfrac
##    Created:	Wed Jul 31, 1992
##    Author:	Bruno Salvy
##	<salvy@rully.inria.fr>
##
## Description: full partial fraction decomposition.
## Based on an idea in M. Bronstein's "Formulas for Series Computations",
## a complete description of the algorithm is in
##  "Full Partial Fraction Decomposition of Rational Functions"
##           Manuel Bronstein and Bruno Salvy.
##               (see proceedings ISSAC'93)

macro(contribution	= `fullparfrac/contribution`);

fullparfrac:=proc (ratfun, x)
local f, p, q, d, i, j, res, g;
    if not type(ratfun,ratpoly(rational,x)) then
	ERROR(`Not a rational function`,ratfun) fi;
    f:=normal(ratfun); d:=0;
    res[0]:=quo(numer(f)/content(denom(f),x,'q'),q,x,'p'); f:=p/q;
    for i while degree(q,x)>0 do
	q:=gcd(q,diff(q,x),'g');
	res[i]:=contribution(f,quo(d,g,x),x,i-1);
	d:=g;
    od;
    RETURN(convert([seq(res[j],j=0..i-1),contribution(f,d,x,i-1)],`+`))
end: # fullparfrac

contribution:=proc(f,d,x,n)
local Dd, i, u, h, tosubs, ptilde, qtilde, hstar, res, A, al;
    if not has(d,x) then RETURN(0) fi;
    # Precompute a list of derivatives of d
    Dd[1]:=diff(d,x);
    for i to n-1 do Dd[i+1]:=diff(Dd[i],x) od;
    tosubs:=[seq(diff(u(x),x$(n-i))=Dd[n-i+1]/(n-i+1),i=1..n-1),u(x)=Dd[1]];
    h:=normal(f*d^n)/u(x)^n;
    for i from 0 to n-1 do
	gcd(op(subs(tosubs,[numer(h),denom(h)])),'ptilde','qtilde');
	if i<>n-1 then h:=normal(diff(h,x)/(i+1)) fi;
	gcd(d,ptilde,'hstar');
	if not has(hstar,x) then res[i]:=NULL
	else
	    gcdex(qtilde,hstar,x,'A');
	    al:=RootOf(hstar,x);
	    if type(al,RootOf) then # this test should not be needed
		res[i]:=Sum(subs(x='alpha',rem(ptilde*A,hstar,x))/
			(x-'alpha')^(n-i),'alpha'=al)
	    else res[i]:=subs(x=al,rem(ptilde*A,hstar,x))/(x-al)^(n-i) fi
	fi
    od;
    RETURN(convert([seq(res[i],i=0..n-1)],`+`))
end: # contribution

# Test:
# f:=(2*x^7-7*x^5+26*x^3+8*x)/(x^8-5*x^6+6*x^4+4*x^2-8);
# fullparfrac(f,x);
# g:=x^3/(x^21+2*x^20+4*x^19+7*x^18+10*x^17+17*x^16+22*x^15+30*x^14+36*x^13+40
# *x^12+47*x^11+46*x^10+49*x^9+43*x^8+38*x^7+32*x^6+23*x^5+19*x^4+10*x^
# 3+7*x^2+2*x+1);
# fullparfrac(g,x);
# h:=36/(x^5-2*x^4-2*x^3+4*x^2+x-2);
# fullparfrac(h,x);

`help/text/fullparfrac` := TEXT(
`FUNCTION: fullparfrac - full partial fraction decomposition`,
`CALLING SEQUENCE:`,
`   fullparfrac(expr,var)`,
`PARAMETERS:`,
`   expr - a rational function over the rationals`,
`   var  - its variable`,
`SYNOPSIS:   `,
`- The fullparfrac function computes a partial fraction decomposition of expr`,
`over C, i.e. the denominators are linear.`,
`EXAMPLES:   `,
`> f:=expand(x^5*(x-1))/expand((x^2+x+1)^2*(x-2)^3);`,
`   `,
`                                        6    5`,
`                                       x  - x`,
`                 f := ----------------------------------------`,
`                       7      6      5      3      2`,
`                      x  - 4 x  + 3 x  + 9 x  - 6 x  - 4 x - 8`,
`   `,
`> fullparfrac(f,x);`,
`       /              179          135\\   /            37           20 \\`,
`       |  -----    - ---- alpha + ----|   |  -----    ---- alpha + ----|`,
`       |   \\         2401         2401|   |   \\       1029         1029|`,
`       |    )      -------------------| + |    )      -----------------|`,
`       |   /            x - alpha     |   |   /                     2  |`,
`       |  -----                       |   |  -----       (x - alpha)   |`,
`       \\alpha = %1                    /   \\alpha = %1                  /`,
`   `,
`                  1952            464            32`,
`            + ------------ + ------------ + -----------`,
`              2401 (x - 2)              2             3`,
`                             343 (x - 2)    49 (x - 2)`,
`   `,
`                                     2`,
`%1 :=                       RootOf(_Z  + _Z + 1)`,
`   `
):