Tag: multiplicities

Complex Factors of Real and Complex Polynomials – part 2

Last week I posted a Maxima utility to factor polynomials using complex roots.  The idea was to use the built-in solver solve to find each root a with multiplicity m, and then construct a factor of the form (z-a)^m.  Typically, solve can fail to identify some roots that are returned by the more robust to_poly_solve.  My problem at that time was that I didn’t know how to access the multiplicities for the roots returned by to_poly_solve.

I’ve written a utility to determine multiplicities for the roots returned by to_poly_solve. With that, here’s an updated version of factorC that takes advantage of the more robust solver.

/* factor a polynomial using all its complex roots */
factorC(_f,_z):=block(
[s,n,m,fp,j],
fp:1,
ss:to_poly_solve(_f,_z),
s:create_list(args(ss)[k][1],k,1,length(ss)),
m:tpsmult1(_f,_z,ss),
/*These lines were needed before I figured out how
 how to handle multiplicities with to_poly_solve
s:solve(_f,_z),
m:multiplicities,*/
n:length(s),
for j:1 thru n do 
 if lhs(s[j])#0 
 then fp:fp*(_z-(rhs(s[j])))^m[j],
 fp:fp*divide(_f,fp)[1],
fp
);
Advertisement

Roots, Multiplicities, and to_poly_solve

A few days ago I posted some code to factor and perform partial fraction decomposition using all available complex roots.  The approach is to find roots a — taking note of the multiplicity m of each —  and then construct factors of the form (z-a)^m.

As I’ve noted before, the built-in command solve can fail to find solutions.  The more robust Maxima solver is to_poly_solve.    That said, after a call to to_poly_solve, I don’t know how to get my hands on the multiplicities.  So, I’ve written a little utility to compute and return a list of multiplicities for the roots returned by to_poly_solve.  Here’s an example, followed by the code.  Notice that because to_poly_solve returns solutions using the %union construction, extracting the jth root requires something like root:rhs(part(s,j,1))

tpsmult1

/* tpsmult1 assumes results of to_poly_solve applied to _f, 
 are passed in the argument s */
tpsmult1(_f,_z,s):=block(
[d,root,nroots,mult,m1],
nroots:length(s),
mult:[],
for j:1 thru nroots do (
root:rhs(part(s,j,1)),
m1:rootmult(_f,_z,root),
mult: endcons(m1,mult)
),
mult
);

/* find the multiplicity of a root _a for the polynomial _f */
rootmult(_f,_z,_a):=block(
[d,val,m,mm],
m:1,
d:divide(_f,(_z-_a)),
val:subst(_z=_a,part(d,1)),
if (val=0) then (
 mm:rootmult(part(d,1),_z,_a),
 m:m+mm
 ),
m
);