--- a/scripts/polynomial/residue.m
+++ b/scripts/polynomial/residue.m
@@ -1,130 +1,131 @@
-# Copyright (C) 1996 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.
+### Copyright (C) 1996 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
+  ## 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.
+  ## 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.
+  ## 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])");
@@ -134,7 +135,7 @@
     toler = .001;
   endif
 
-# Make sure both polynomials are in reduced form.
+  ## Make sure both polynomials are in reduced form.
 
   a = polyreduce (a);
   b = polyreduce (b);
@@ -145,7 +146,7 @@
   la = length (a);
   lb = length (b);
 
-# Handle special cases here.
+  ## Handle special cases here.
 
   if (la == 0 || lb == 0)
     k = r = p = e = [];
@@ -156,22 +157,22 @@
     return;
   endif
 
-# Find the poles.
+  ## Find the poles.
 
   p = roots (a);
   lp = length (p);
 
-# Determine if the poles are (effectively) real.
+  ## 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.
+  ## Find the direct term if there is one.
 
   if (lb >= la)
-# Also returns the reduced numerator.
+    ## Also returns the reduced numerator.
     [k, b] = deconv (b, a);
     lb = length (b);
   else
@@ -185,16 +186,16 @@
   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.
+  ## 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.
+  ## 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);
@@ -209,13 +210,13 @@
 
   while (new_p_index <= lp)
     if (abs (p (new_p_index) - p (old_p_index)) < toler)
-# We've found a multiple pole.
+      ## 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.
+      ## Remove the pole from pr.
       pr (pr_index) = [];
     else
-# It's a different pole.
+      ## It's a different pole.
       D++;
       M_index++;
       M (M_index) = 1;
@@ -225,15 +226,15 @@
     new_p_index++;
   endwhile
 
-# Shorten M to it's proper length
+  ## Shorten M to it's proper length
 
   M = M (1:D);
 
-# Now set up the polynomial matrices.
+  ## Now set up the polynomial matrices.
 
   MM = max(M);
 
-# Left hand side polynomial
+  ## Left hand side polynomial
 
   lhs = zeros (MM, lb);
   rhs = zeros (MM, lp*lp);
@@ -268,14 +269,14 @@
     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 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.
+  ## 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);
@@ -296,7 +297,7 @@
     dpi++;
   endwhile
 
-# Solve for the residues.
+  ## Solve for the residues.
 
   r = A \ B;