Re: polynomial partial fractions

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* Kevin Ryde [2014-12-11 09:38]:
> I had the urge to break some polynomial ratios into partial fractions.
> I thought maybe a vector result of terms which sum to the original.
> 
>     p = x^4 / ((1-x)*(1-2*x)*(1 - x - 2*x^3))
>     v = polynomial_crack_into_partial_fractions(p)
>     v == [ (1/2) / (1-x),
>            (1/2) / (1-2*x),
>            - (1 + (1/2)*x + x^2) / (1-x-2*x^3) ]
>     vecsum(v) == p
> 
> What would I look at for such a thing?  I know how to build a matrix for
> matsolve() to give the numerators, but perhaps this exists already.
> I saw Henri Cohen's cohen.gp ratdec() but it seems to go polroots()
> where I had in mind only going as far as can be factorized exactly over
> complex or quadratics (so leave the cubic above unchanged).  Could
> supply the desired denominators if necessary.

Here's a script providing the polynomial part and the numerators, sorted
according to the irreducible factors of the denominator

Untested :-)

padic(z,p)=
{ my(n = floor(poldegree(z)/poldegree(p))+1);
  my(i = n, v = vector(n));
  while(z, [z,v[i]] = divrem(z,p); i--);
  v;
}

partial_fractions(z)=
{
  my(A, B, C, Q);
  A = numerator(z);
  B = denominator(z);
  [Q,A] = divrem(A,B);
  my(F = factor(B), P = F[,1], E = F[,2]);
  my(PE = vector(#P, i, P[i]^E[i]));
  C = vector(#P, i, Mod(B \/ PE[i], PE[i]));
  [Q, vector(#P, i, padic(lift(A/C[i]), P[i])), PE];
}


(10:39) gp > partial_fractions(p)
%1 = [0, [[-1/2], [-1/2], [x^2 + 1/2*x + 1]], [x - 1, 2*x - 1, 2*x^3 + x - 1]]

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux1.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux1.fr/  [PARI/GP]
`