--- a/scripts/ChangeLog	Fri Oct 05 22:40:59 2007 +0000
+++ b/scripts/ChangeLog	Sat Oct 06 04:31:18 2007 +0000
@@ -1,3 +1,13 @@
+2007-10-05  Ben Abbott <bpabbott@mac.com>
+
+	* polynomial/mpoles.m: New function.
+	* polynomial/residue.m: Modified to behave in reciprocal
+	manner.  No longer compute 4th output, "e".  No longer accept
+	tolerance input.  Explicitly set tolerance parameter to 0.001.
+	Respect maximum relative difference in poles when determining
+	their multiplicity.  Use mpoles to determine the multiplicity of
+	poles.
+
 2007-10-05  Peter A. Gustafson  <petegus@umich.edu>
 
 	* plot/__go_draw_axes__.m: Add cbrange to the plot stream 

--- a/scripts/polynomial/Makefile.in	Fri Oct 05 22:40:59 2007 +0000
+++ b/scripts/polynomial/Makefile.in	Sat Oct 06 04:31:18 2007 +0000
@@ -20,7 +20,7 @@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_DATA = @INSTALL_DATA@
 
-SOURCES = compan.m conv.m deconv.m mkpp.m pchip.m poly.m \
+SOURCES = compan.m conv.m deconv.m mkpp.m mpoles.m pchip.m poly.m \
   polyder.m polyderiv.m polyfit.m polygcd.m polyint.m \
   polyout.m polyreduce.m polyval.m polyvalm.m ppval.m residue.m \
   roots.m spline.m unmkpp.m

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/mpoles.m	Sat Oct 06 04:31:18 2007 +0000
@@ -0,0 +1,112 @@
+## Copyright (C) 2007 Ben Abbott
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p})
+## @deftypefnx {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p}, @var{tol})
+## @deftypefnx {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p}, @var{tol}, @var{reorder})
+## Identifiy unique poles in @var{p} and associates their multiplicity,
+## ordering them from largest to smallest.
+## 
+## If the relative difference of the poles is less than @var{tol}, then
+## they are considered to be multiples.  The default value for @var{tol}
+## is 0.001.
+##
+## If the optional parameter @var{reorder} is zero, poles are not sorted.
+##
+## The value @var{multp} is a vector specifying the multiplicity of the
+## poles.  @var{multp}(:) refers to mulitplicity of @var{p}(@var{indx}(:)).
+##
+## For example,
+##
+## @example
+## @group
+## p = [2 3 1 1 2];
+## [m, n] = mpoles(p);
+##   @result{} m = [1; 1; 2; 1; 2]
+##   @result{} n = [2; 5; 1; 4; 3]
+##   @result{} p(n) = [3, 2, 2, 1, 1]
+## @end group
+## @end example
+##
+## @seealso{poly, roots, conv, deconv, polyval, polyderiv, polyinteg, residue}
+## @end deftypefn
+
+## Author: Ben Abbott <bpabbott@mac.com>
+## Created: Sept 30, 2007
+
+function [multp, indx] = mpoles (p, tol, reorder)
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+  endif
+
+   if (nargin < 2 || isempty (tol))
+     tol = 0.001;
+   endif
+
+   if (nargin < 3 || isempty (reorder))
+     reorder = true;
+   endif
+
+  Np = numel (p);
+
+  ## Force the poles to be a column vector.
+
+  p = p(:);
+
+  ## Sort the poles according to their magnitidues, largest first.
+
+  if (reorder)
+    ## Sort with smallest magnitude first.
+    [p, ordr] = sort (p);
+    ## Reverse order, largest maginitude first.
+    n = Np:-1:1;
+    p = p(n);
+    ordr = ordr(n);
+  else
+    ordr = 1:Np;
+  endif
+
+  ## Find pole multiplicty by comparing the relative differnce in the
+  ## poles.
+
+  multp = zeros (Np, 1);
+  indx = [];
+  n = find (multp == 0, 1);
+  while (n)
+    dp = abs (p-p(n));
+    if (p(n) == 0.0)
+      p0 = mean (abs (p(find (abs (p) > 0))));
+      if (isempty (p0))
+        p0 = 1;
+      end
+    else
+      p0 = abs (p(n));
+    endif
+    k = find (dp < tol * p0);
+    m = 1:numel (k);
+    multp(k) = m;
+    indx = [indx; k];
+    n = find (multp == 0, 1);
+  endwhile
+  multp = multp(indx);
+  indx = indx(ordr);
+
+endfunction

--- a/scripts/polynomial/residue.m	Fri Oct 05 22:40:59 2007 +0000
+++ b/scripts/polynomial/residue.m	Sat Oct 06 04:31:18 2007 +0000
@@ -1,4 +1,5 @@
 ## Copyright (C) 1996, 1997 John W. Eaton
+## Copyright (C) 2007 Ben Abbott
 ##
 ## This file is part of Octave.
 ##
@@ -18,17 +19,22 @@
 ## 02110-1301, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} residue (@var{b}, @var{a}, @var{tol})
-## If @var{b} and @var{a} are vectors of polynomial coefficients, then
-## residue calculates the partial fraction expansion corresponding to the
-## ratio of the two polynomials.
+## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}] =} residue (@var{b}, @var{a})
+## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k})
+## In the instance of two inputs, they are assumed to correspond to vectors
+## of polynomial coefficients, @var{b} and @var{a}. From these polynomial
+## coefficients, the function residue calculates the partial fraction 
+## expansion corresponding to the ratio of the two polynomials.
 ## @cindex partial fraction expansion
 ##
-## The function @code{residue} returns @var{r}, @var{p}, @var{k}, and
-## @var{e}, where the vector @var{r} contains the residue terms, @var{p}
-## contains the pole values, @var{k} contains the coefficients of a direct
-## polynomial term (if it exists) and @var{e} is a vector containing the
-## powers of the denominators in the partial fraction terms.
+## In this instance, the function @code{residue} returns @var{r}, @var{p},
+## and @var{k}, where the vector @var{r} contains the residue terms,
+## @var{p} contains the pole values, @var{k} contains the coefficients of a 
+## direct polynomial term (if it exists).
+## 
+## In the instance of three inputs, the function 'residue' performs the
+## reciprocal task. Meaning the partial fraction expansion is reconstituted
+## into the corresponding pair of polynomial coefficients.
 ##
 ## Assuming @var{b} and @var{a} represent polynomials
 ## @iftex
@@ -60,25 +66,19 @@
 ## @noindent
 ## where @math{M} is the number of poles (the length of the @var{r},
 ## @var{p}, and @var{e} vectors) and @math{N} is the length of the
-## @var{k} vector.
+## @var{k} vector. The @var{e} vector specifies the multiplicity of the
+## mth residue's pole.
 ##
-## The argument @var{tol} is optional, and if not specified, a default
-## value of 0.001 is assumed.  The tolerance value is used to determine
-## whether poles with small imaginary components are declared real.  It is
-## also used to determine if two poles are distinct.  If the ratio of the
-## imaginary part of a pole to the real part is less than @var{tol}, the
-## imaginary part is discarded.  If two poles are farther apart than
-## @var{tol} they are distinct.  For example,
+## For example,
 ##
 ## @example
 ## @group
-##  b = [1, 1, 1];
-##  a = [1, -5, 8, -4];
-##  [r, p, k, e] = residue (b, a);
+## b = [1, 1, 1];
+## a = [1, -5, 8, -4];
+## [r, p, k] = residue (b, a);
 ##      @result{} r = [-2, 7, 3]
 ##      @result{} p = [2, 2, 1]
 ##      @result{} k = [](0x0)
-##      @result{} e = [1, 2, 1]
 ## @end group
 ## @end example
 ##
@@ -98,84 +98,61 @@
 ##    ------------------- = ----- + ------- + -----
 ##    s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
 ## @end example
+##
 ## @end ifinfo
-## @seealso{poly, roots, conv, deconv, polyval, polyderiv, and polyinteg}
+## A similar, but reciprocal example, where the fraction's polynomials are 
+## reconstituted from the residues, poles, and direct term is
+##
+## @example
+## @group
+## r = [-2, 7, 3];
+## p = [2, 2, 1];
+## k = [1 0];
+## [b, a] = residue (r, p, k);
+##      @result{} b = [1, -5, 9, -3, 1]
+##      @result{} a = [1, -5, 8, 4]
+## @end group
+## @end example
+##
+## @noindent
+## which implies the following partial fraction expansion
+## @iftex
+## @tex
+## $$
+## {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+##
+## @example
+##    s^4 - 5s^3 + 9s^2 - 3s + 1    -2        7        3
+##    -------------------------- = ----- + ------- + ----- + s
+##        s^3 - 5s^2 + 8s - 4      (s-2)   (s-2)^2   (s-1)
+## @end example
+## @end ifinfo
+## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg}
 ## @end deftypefn
 
 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
+## Author: Ben Abbott <bpabbott@mac.com>
 ## Created: June 1994
 ## Adapted-By: jwe
 
-function [r, p, k, e] = residue (b, a, toler)
-
-  ## Here's the method used to find the residues.
-  ## The partial fraction expansion can be written as:
-  ##
-  ##
-  ##   P(s)    D   M(k)      A(k,m)
-  ##   ---- =  #    #    -------------
-  ##   Q(s)   k=1  m=1   (s - pr(k))^m
-  ##
-  ## (# is used to represent a summation) where D is the number of
-  ## distinct roots, pr(k) is the kth distinct root, M(k) is the
-  ## multiplicity of the root, and A(k,m) is the residue cooresponding
-  ## to the kth distinct root with multiplicity m.  For example,
-  ##
-  ##        s^2            A(1,1)   A(2,1)    A(2,2)
-  ## ------------------- = ------ + ------ + -------
-  ## s^3 + 4s^2 + 5s + 2    (s+2)    (s+1)   (s+1)^2
-  ##
-  ## In this case there are two distinct roots (D=2 and pr = [-2 -1]),
-  ## the first root has multiplicity one and the second multiplicity
-  ## two (M = [1 2]) The residues are actually stored in vector format as
-  ## r = [ A(1,1) A(2,1) A(2,2) ].
-  ##
-  ## We then multiply both sides by Q(s).  Continuing the example:
-  ##
-  ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
-  ##
-  ## or
-  ##
-  ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
-  ##
-  ## The coefficients of the polynomials on the right are stored in a row
-  ## vector called rhs, while the coefficients of the polynomial on the
-  ## left is stored in a row vector called lhs.  If the multiplicity of
-  ## any root is greater than one we'll also need derivatives of this
-  ## equation of order up to the maximum multiplicity minus one.  The
-  ## derivative coefficients are stored in successive rows of lhs and
-  ## rhs.
-  ##
-  ## For our example lhs and rhs would be:
-  ##
-  ##       | 1 0 0 |
-  ## lhs = |       |
-  ##       | 0 2 0 |
-  ##
-  ##       | 1 2 1 1 3 2 0 1 2 |
-  ## rhs = |                   |
-  ##       | 0 2 2 0 2 3 0 0 1 |
-  ##
-  ## We then form a vector B and a matrix A obtained by evaluating the
-  ## polynomials in lhs and rhs at the pole values. If a pole has a
-  ## multiplicity greater than one we also evaluate the derivative
-  ## polynomials (successive rows) at the pole value.
-  ##
-  ## For our example we would have
-  ##
-  ## | 4|   | 1 0 0 |   | r(1) |
-  ## | 1| = | 0 0 1 | * | r(2) |
-  ## |-2|   | 0 1 1 |   | r(3) |
-  ##
-  ## We then solve for the residues using matrix division.
+function [r, p, k] = residue (b, a, varargin)
 
   if (nargin < 2 || nargin > 3)
     print_usage ();
   endif
 
-  if (nargin == 2)
-    toler = .001;
-  endif
+  toler = .001;
+
+  if (nargin == 3)
+    ## The inputs are the residue, pole, and direct part. Solve for the
+    ## corresponding numerator and denominator polynomials
+    [r, p] = rresidue (b, a, varargin{1}, toler);
+    return
+  end
 
   ## Make sure both polynomials are in reduced form.
 
@@ -206,17 +183,26 @@
 
   ## Determine if the poles are (effectively) zero.
 
-  p(abs (p) < toler) = 0;
+  small = max (abs (p));
+  small = max ([small, 1] ) * 1e-8 * (1 + numel (p))^2;
+  p(abs (p) < small) = 0;
+
+  ## Determine if the poles are (effectively) real, or imaginary.
 
-  ## Determine if the poles are (effectively) real.
+  index = (abs (imag (p)) < small);
+  p(index) = real (p(index));
+  index = (abs (real (p)) < small);
+  p(index) = 1i * imag (p(index));
 
-  index = (abs (p) >= toler & (abs (imag (p)) ./ abs (p) < toler));
-  p(index) = real (p(index));
+  ## Sort poles so that multiplicity loop will work.
+
+  [e, indx] = mpoles (p, toler, 1);
+  p = p (indx);
 
   ## Find the direct term if there is one.
 
   if (lb >= la)
-    ## Also returns the reduced numerator.
+    ## Also return the reduced numerator.
     [k, b] = deconv (b, a);
     lb = length (b);
   else
@@ -225,123 +211,106 @@
 
   if (lp == 1)
     r = polyval (b, p);
-    e = 1;
     return;
   endif
 
-  ## We need to determine the number and multiplicity of the roots.
-  ##
-  ## D  is the number of distinct roots.
-  ## M  is a vector of length D containing the multiplicity of each root.
-  ## pr is a vector of length D containing only the distinct roots.
-  ## e  is a vector of length lp which indicates the power in the partial
-  ##    fraction expansion of each term in p.
-
-  ## Set initial values.  We'll remove elements from pr as we find
-  ## multiplicities.  We'll shorten M afterwards.
-
-  e = ones (lp, 1);
-  M = zeros (lp, 1);
-  pr = p;
-  D = 1;
-  M(1) = 1;
-
-  old_p_index = 1;
-  new_p_index = 2;
-  M_index = 1;
-  pr_index = 2;
+  ## Determine the order of the denominator and remaining numerator.
+  ## With the direct term removed the potential order of the numerator
+  ## is one less than the order of the denominator.
 
-  while (new_p_index <= lp)
-    if (abs (p (new_p_index) - p (old_p_index)) < toler)
-      ## We've found a multiple pole.
-      M (M_index) = M (M_index) + 1;
-      e (new_p_index) = e (new_p_index-1) + 1;
-      ## Remove the pole from pr.
-      pr (pr_index) = [];
-    else
-      ## It's a different pole.
-      D++;
-      M_index++;
-      M (M_index) = 1;
-      old_p_index = new_p_index;
-      pr_index++;
-    endif
-    new_p_index++;
-  endwhile
-
-  ## Shorten M to it's proper length
-
-  M = M (1:D);
-
-  ## Now set up the polynomial matrices.
-
-  MM = max(M);
-
-  ## Left hand side polynomial
+  aorder = numel (a) - 1;
+  border = aorder - 1;
 
-  lhs = zeros (MM, lb);
-  rhs = zeros (MM, lp*lp);
-  lhs (1, :) = b;
-  rhi = 1;
-  dpi = 1;
-  mpi = 1;
-  while (dpi <= D)
-    for ind = 1:M(dpi)
-      if (mpi > 1 && (mpi+ind) <= lp)
-        cp = [p(1:mpi-1); p(mpi+ind:lp)];
-      elseif (mpi == 1)
-        cp = p (mpi+ind:lp);
-      else
-        cp = p (1:mpi-1);
-      endif
-      rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp, 0, 2);
-      rhi = rhi + lp;
-    endfor
-    mpi = mpi + M (dpi);
-    dpi++;
-  endwhile
-  if (MM > 1)
-    for index = 2:MM
-      lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb, 0, 2);
-      ind = 1;
-      for rhi = 1:lp
-        cp = rhs (index-1, ind:ind+lp-1);
-        rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp, 0, 2);
-        ind = ind + lp;
-      endfor
-    endfor
-  endif
+  ## Construct a system of equations relating the individual
+  ## contributions from each residue to the complete numerator.
 
-  ## Now lhs contains the numerator polynomial and as many derivatives as
-  ## are required.  rhs is a matrix of polynomials, the first row
-  ## contains the corresponding polynomial for each residue and
-  ## successive rows are derivatives.
-
-  ## Now we need to evaluate the first row of lhs and rhs at each
-  ## distinct pole value.  If there are multiple poles we will also need
-  ## to evaluate the derivatives at the pole value also.
-
-  B = zeros (lp, 1);
-  A = zeros (lp, lp);
-
-  dpi = 1;
-  row = 1;
-  while (dpi <= D)
-    for mi = 1:M(dpi)
-      B (row) = polyval (lhs (mi, :), pr (dpi));
-      ci = 1;
-      for col = 1:lp
-        cp = rhs (mi, ci:ci+lp-1);
-        A (row, col) = polyval (cp, pr(dpi));
-        ci = ci + lp;
-      endfor
-      row++;
-    endfor
-    dpi++;
-  endwhile
+  A = zeros (border+1, border+1);
+  B = prepad (reshape (b, [numel(b), 1]), border+1, 0);
+  for ip = 1:numel(p)
+    ri = zeros (size (p));
+    ri(ip) = 1;
+    A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).';
+  endfor
 
   ## Solve for the residues.
 
   r = A \ B;
 
 endfunction
+
+function [pnum, pden] = rresidue (r, p, k, toler)
+
+  ## Reconstitute the numerator and denominator polynomials from the
+  ## residues, poles, and direct term.
+
+  if (nargin < 2 || nargin > 4)
+    print_usage ();
+  endif
+
+  if (nargin < 4)
+    toler = [];
+  endif
+
+  if (nargin < 3)
+    k = [];
+  endif
+  
+  [multp, indx] = mpoles (p, toler, 0);
+
+  p = p (indx);
+  r = r (indx);
+
+  indx = 1:numel(p);
+
+  for n = indx
+    pn = [1, -p(n)];
+    if n == 1
+      pden = pn;
+    else
+      pden = conv (pden, pn);
+    endif
+  endfor
+
+  ## D is the order of the denominator
+  ## K is the order of the direct polynomial
+  ## N is the order of the resulting numerator
+  ## pnum(1:(N+1)) is the numerator's polynomial
+  ## pden(1:(D+1)) is the denominator's polynomial
+  ## pm is the multible pole for the nth residue
+  ## pn is the numerator contribution for the nth residue
+
+  D = numel (pden) - 1;
+  K = numel (k) - 1;
+  N = K + D;
+  pnum = zeros (1, N+1);
+  for n = indx(abs(r)>0)
+    p1 = [1, -p(n)];
+    for m = 1:multp(n)
+      if m == 1
+        pm = p1;
+      else
+        pm = conv (pm, p1);
+      endif
+    endfor
+    pn = deconv (pden, pm);
+    pn = r(n) * pn;
+    pnum = pnum + prepad ( pn, N+1, 0);
+  endfor
+
+  ## Add the direct term.
+
+  if (numel (k))
+    pnum = pnum + conv (pden, k);
+  endif
+
+  ## Check for leading zeros and trim the polynomial coefficients.
+
+  small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
+
+  pnum (abs (pnum) < small) = 0;
+  pden (abs (pden) < small) = 0;
+
+  pnum = polyreduce (pnum);
+  pden = polyreduce (pden);
+
+endfunction