--- a/scripts/polynomial/residue.m
+++ b/scripts/polynomial/residue.m
@@ -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