# Copyright (C) 1995 John W. Eaton
# 
# 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

function [r, p, k, e] = residue (b, a, toler)

# usage: [r, p, k, e] = residue (b, a)
#
# If b and a are vectors of polynomial coefficients, then residue
# calculates the partial fraction expansion corresponding to the
# ratio of the two polynomials. The vector r contains the residue
# terms, p contains the pole values, k contains the coefficients of
# a direct polynomial term (if it exists) and e is a vector containing
# the powers of the denominators in the partial fraction terms.
# Assuming b and a represent polynomials P(s) and Q(s) we have:
#
#  P(s)    M       r(m)         N
#  ---- =  #  -------------  +  # k(n)*s^(N-n)
#  Q(s)   m=1 (s-p(m))^e(m)    n=1
#
# (# represents summation) where M is the number of poles (the length of
# the r, p, and e vectors) and N is the length of the k vector.
#
# [r p k e] = residue(b,a,tol)
#
# This form of the function call may be used to set a tolerance value.
# The default value is 0.001. 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 tol, the
# imaginary part is discarded. If two poles are farther apart than tol
# they are distinct.
#
# Example:
#  b = [1,  1, 1];
#  a = [1, -5, 8, -4];
#
#  [r, p, k, e] = residue (b, a)
#
#  returns
#
#  r = [-2, 7, 3]; p = [2, 2, 1]; k = []; e = [1, 2, 1];
#
#  which implies the following partial fraction expansion
#
#        s^2 + s + 1         -2       7        3
#    ------------------- = ----- + ------- + -----
#    s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
#
# SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg

# Written by Tony Richardson (amr@mpl.ucsd.edu) June 1994.

# 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.

  if (nargin < 2 || nargin > 3)
    usage ("residue (b, a [, toler])");
  endif

  if (nargin == 2)
    toler = .001;
  endif

# Make sure both polynomials are in reduced form.

  a = polyreduce (a);
  b = polyreduce (b);

  b = b / a(1);
  a = a / a(1);

  la = length (a);
  lb = length (b);

# Handle special cases here.

  if (la == 0 || lb == 0)
    k = r = p = e = [];
    return;
  elseif (la == 1)
    k = b / a;
    r = p = e = [];
    return;
  endif

# Find the poles.

  p = roots (a);
  lp = length (p);

# Determine if the poles are (effectively) real.

  index = find (abs (imag (p) ./ real (p)) < toler);
  if (length (index) != 0)
    p (index) = real (p (index));
  endif

# Find the direct term if there is one.

  if (lb >= la)
# Also returns the reduced numerator.
    [k, b] = deconv (b, a);
    lb = length (b);
  else
    k = [];
  endif

  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;

  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

  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);
      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);
      ind = 1;
      for rhi = 1:lp
        cp = rhs (index-1, ind:ind+lp-1);
        rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp);
        ind = ind + lp;
      endfor
    endfor
  endif

# 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

# Solve for the residues.

  r = A \ B;

endfunction