Next: Help File for Up: Stability of Polynomials Previous: Mathematical Background.

## Source Code for Hurwitz

```##########################################################################
# (C) Copyright 1989 by Robert M. Corless
#
# Maintenance history
#   --- Originally written in 1989 for AM 511 class demonstration.
#   --- Rewritten Nov 1991 taking into account the share library
#       version by Dominik Gruntz.  In particular, existing code
#       was re-named and the new routine "Hurwitz" was added to
#       decide (when it could) whether a polynomial was Hurwitz or not.
#
#   --- Nov 20, 1991 final changes before putting into share library:
#       saved into file Hurwitz.m, userinfo added, and polynomial is
#       on input converted to rational to fix a multivariate float -
#       gcd weakness.
#
#   --- Future improvements:  need to allow user-specified normalizer.
#
##########################################################################

`Hurwitz/paraconjugate` := proc(f) local e, z;
e := evalc @ conjugate;
unapply(e(f(e(-z))),z);
end:

`Hurwitz/Stieltjes_Fraction` := proc(num,den,z) local q, r;
if den = 0 then NULL
else
q := quo(num,den,z,'r');
evalc(q),`Hurwitz/Stieltjes_Fraction`(den,r,z)
fi
end:

`Hurwitz/data` := proc(expr,z) local g,f,fs;
option remember;  # This proc has option remember because it
# is reasonably likely to be called twice
# by Hurwitz, once to get FAIL and once for 's'.

f := unapply(expr,z):
fs:= `Hurwitz/paraconjugate`(f);
g := evala(Gcd(convert(fs(z),RootOf),convert(f(z),RootOf)));
if nargs > 2 then assign(args[3],g) fi;
if degree(g,z) <>0 then
[[],g]
else
[[`Hurwitz/Stieltjes_Fraction`(f(z)-fs(z),f(z)+fs(z),z)],g];
fi
end:

`Hurwitz/term_ok` := proc(term,z) local a,b,test;
if degree(term,z)<1 then
false
elif degree(term,z) > 1 then
if type(lcoeff(term,z),numeric) then
false
else
FAIL
fi
else
a := coeff(term,z,0); b := coeff(term,z,1);
test := evalb(evalc(Re(a))=0);
if not(test=true or test=false) then test := FAIL fi;
if test then
test := evalb(b>0);
if not(test=true or test=false) then test := FAIL fi
fi;
test
fi
end:

Hurwitz := proc(poly,z) local p,g,a,b,test,i;
if not type(z,name) then ERROR(`2nd argument must be a name`) fi;
if not type(poly,polynom(anything,z)) then
ERROR(`1st argument must be univariate polynomial`) fi;

`Hurwitz/data`(convert(poly,'rational','exact'),z);
p := "[1]; g := ""[2];

userinfo(1,Hurwitz,`Continued fraction computed`);
userinfo(1,Hurwitz,`Starting stability testing...`);

if nargs > 2 then assign(args[3],p); fi;
if nargs > 3 then assign(args[4],g); fi;

if degree(g,z) > 1 then
if type(lcoeff(g,z),numeric) then
false
else
FAIL
fi

elif nops(p) < 1 then
userinfo(1,Hurwitz,`Conditions vacuously satisfied --- check manually`);
true

elif degree(p[1],z) > 1 then
if type(lcoeff(p[1],z),numeric) then
false
else
FAIL
fi

else # p[1] = a + b*z, possibly a = 0 or b=0
a := coeff(p[1],z,0);  b := coeff(p[1],z,1);
test := evalb(evalc(Re(a))=0);
if not(test=true or test=false) then test := FAIL fi;

if test then
test := evalb(b >= 0);
if not(test=true or test=false) then test := FAIL fi
fi;

for i from 2 to nops(p) while test do
test := `Hurwitz/term_ok`(p[i],z)
od;
test
fi
end:
```

Robert Corless
Wed Feb 28 18:25:56 EST 1996