Karim Belabas on Thu, 11 Dec 2014 10:44:32 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: polynomial partial fractions |
* 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] `