Re: Partial Fraction Decomposition with imaginary coeff.

*To*: mathgroup at yoda.physics.unc.edu*Subject*: Re: Partial Fraction Decomposition with imaginary coeff.*From*: Brian Evans <evans at gauss.eedsp.gatech.edu>*Date*: Sun, 8 Mar 1992 22:54:48 GMT

You just don't know how difficult Apart has made my life. It is crucial to the implementation of inverse linear transforms. Under 1.2, it would not handle polynomials with rational coefficients expressed in decimal form. Under 1.2 and 2.0, Apart does not break up terms like 1 / (x^2 + 1). About two years ago, I had to write a general purpose routine to work around Apart's drawbacks. I call it MyApart. It is embedded in the signal processing packages for Mathematica and is used by the inverse z- and Laplace transforms when Apart doesn't complete the decomposition. MyApart is, of course, darn slow. Brian L. Evans Digital Signal Processing Laboratory School of Electrical Engineering Georgia Institute of Technology Atlanta, GA 30332-0250 e-mail: evans at eedsp.gatech.edu ----------------------------- cut here ---------------------------------- (* :Title: Miscellaneous Routines *) (* :Authors: Brian Evans, James McClellan *) (* :Summary: To provide a robust partial fractions routine. *) (* :Context: Misc` *) (* :PackageVersion: 2.6 *) (* :Copyright: Copyright 1992 by Brian L. Evans Georgia Tech Research Corporation Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee is hereby granted, provided that the above copyright notice appear in all copies and that both that copyright notice and this permission notice appear in supporting documentation, and that the name of the Georgia Tech Research Corporation, Georgia Tech, or Georgia Institute of Technology not be used in advertising or publicity pertaining to distribution of the software without specific, written prior permission. Georgia Tech makes no representations about the suitability of this software for any purpose. It is provided "as is" without express or implied warranty. *) (* :History: *) (* :Keywords: partial fractions decomposition *) (* :Source: *) (* :Warning: *) (* :Mathematica Version: 1.2 or 2.0 *) (* :Limitation: *) (* :Functions: GetRoot GetRootList MyApart *) (* B E G I N P A C K A G E *) BeginPackage[ "Misc`" ] (* U S A G E I N F O R M A T I O N *) GetRoot::usage = "GetRoot[rule] extracts the value from an expression like \ {z -> 0.}, which is 0. in this case." GetRootList::usage = "GetRootList[p, x] returns a list of the approximate numerical roots \ of expression p, a function of x, with infinite roots removed. \ GetRootList[p, x, filter] applies filter to the list of roots \ returned by the Solve function (defaults to N)." MyApart::usage = "MyApart[ rational_polynomial, x ] decomposes the rational \ polynomial into a sum of fractions whose numerators are of \ the form (x + b)^n where b is a constant and n is an integer. \ MyApart[rational_polynomial, x, filter] specifies a filter to be \ placed on the output of the Solve command used to root the \ denominator: Identity for rational and N for real-valued roots. \ The default is Identity. \ In Mathematica 1.2, MyApart is about 25 times slower than Apart." (* E N D U S A G E I N F O R M A T I O N *) Begin[ "Private`" ] (* GetRoot *) GetRoot[{}] := {} (* no roots *) GetRoot[rule_] := rule[[1]][[2]] (* GetRootList *) GetRootList[p_, x_, filter_:N] := Select[ Map[ GetRoot, filter[ Solve[ p == 0, x ] ] ], (! MatchQ[#1, DirectedInfinity[___]])& ] (* MyApart -- kludge around the way Apart does partial fractions *) (* Root denominator and replace roots with symbols *) MyApart[ratpoly_, x_, filter_:Identity] := Block [ {apart, denom, denomfactored, normfact, numer, partfrac, rootlist, rootmult, rules}, numer = Numerator[ratpoly]; denom = Denominator[ratpoly]; normfact = Last[ CoefficientList[denom, x] ]; numer /= normfact; denom /= normfact; rootlist = Sort[ GetRootList[denom, x, filter] ]; { denomfactored, rules } = multiplicityform[rootlist, x]; apart = Apart[numer / denomfactored, x]; partfrac = apart /. rules; partfrac /. ( a_. / (b_ c_) :> a / ( Together[b] c ) /; FreeQ[b, x] && ! FreeQ[c, x] ) ] multiplicityform[ roots_, x_ ] := Block [ {count = 1, cur, denom = 1, i, last, length, sublist = {}, sym = 1}, Clear[localvar]; (* localvar is global to package *) length = Length[roots]; last = First[roots]; For [ i = 2, i <= length, i++, cur = roots[[i]]; If [ SameQ[ cur, last ], count++, denom *= (x - localvar[sym])^count; PrependTo[ sublist, localvar[sym] -> last ]; sym++; count = 1 ]; last = cur ]; denom *= (x - localvar[sym])^count; PrependTo[ sublist, localvar[sym] -> last ]; { denom, sublist } ] End[] EndPackage[] Protect[ GetRoot, GetRootList, MyApart ] Null