next up previous
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