Mercurial > octave-nkf

changeset 5217:e88886a6934d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-03-16 20:03:01 by jwe]
author jwe
date Wed, 16 Mar 2005 20:03:01 +0000
parents 5ed60b8b1ac4
children 4d036412ccca
files scripts/ChangeLog scripts/polynomial/polyder.m scripts/polynomial/polyderiv.m scripts/polynomial/polygcd.m
diffstat 4 files changed, 87 insertions(+), 77 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Wed Mar 16 19:51:46 2005 +0000
+++ b/scripts/ChangeLog	Wed Mar 16 20:03:01 2005 +0000
@@ -1,3 +1,9 @@
+2005-03-16  Paul Kienzle  <pkienzle@users.sf.net>
+
+	* polynomial/polyderiv.m : Accept a*b, a/b.  Auto-reduce common terms.
+	* polynomial/polyder.m: Ditto.
+        * polynomial/polygcd.m: New function.
+
 2005-03-16  John W. Eaton  <jwe@octave.org>
 
 	* control/base/__stepimp__.m, control/base/bode.m,
--- a/scripts/polynomial/polyder.m	Wed Mar 16 19:51:46 2005 +0000
+++ b/scripts/polynomial/polyder.m	Wed Mar 16 20:03:01 2005 +0000
@@ -25,21 +25,19 @@
 ## @end deftypefn
 
 ## Author: John W. Eaton
-## Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
-##    handle b/a and b*a
 
 function [q, r] = polyder (p, a)
 
   if (nargin == 1)
     q = polyderiv (p);
-  elseif (nargin==2)
-    if (nargout==2)
-      [q, r] = polyderiv (p,a);
+  elseif (nargin == 2)
+    if (nargout == 2)
+      [q, r] = polyderiv (p, a);
     else
-      q = polyderiv (p,a);
+      q = polyderiv (p, a);
     endif
   else
-    usage ("q=polyder(p) or q=polyder(b,a) or [q, r]=polyder(b,a)");
+    usage ("q = polyder (p) or q = polyder (b, a) or [q, r] = polyder (b, a)");
   endif
 
 endfunction
--- a/scripts/polynomial/polyderiv.m	Wed Mar 16 19:51:46 2005 +0000
+++ b/scripts/polynomial/polyderiv.m	Wed Mar 16 20:03:01 2005 +0000
@@ -33,63 +33,61 @@
 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
 ## Created: June 1994
 ## Adapted-By: jwe
-## Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
-##    handle b/a and b*a
 
 function [q, r] = polyderiv (p, a)
 
-  if (nargin < 1 || nargin > 3)
-    usage ("q=polyderiv(p) or q=polyderiv(b,a) or [q, r]=polyderiv(b,a)");
-  endif
-
-  if (! isvector (p))
-    error ("polyderiv: argument must be a vector");
-  endif
-
-  if (nargin == 2)
-    if (! isvector (a))
+  if (nargin == 1 || nargin == 2)
+    if (! isvector (p))
       error ("polyderiv: argument must be a vector");
     endif
-    if (nargout == 1) 
-      ## derivative of p*a returns a single polynomial
-      q = polyderiv(conv(p,a));
+    if (nargin == 2)
+      if (! isvector (a))
+	error ("polyderiv: argument must be a vector");
+      endif
+      if (nargout == 1) 
+	## derivative of p*a returns a single polynomial
+	q = polyderiv (conv (p, a));
+      else
+	## derivative of p/a returns numerator and denominator
+	r = conv (a, a);
+	if (numel (p) == 1)
+	  q = -p * polyderiv (a);
+	elseif (numel (a) == 1)
+	  q = a * polyderiv (p);
+	else
+	  q = conv (polyderiv (p), a) - conv (p, polyderiv (a));
+	  q = polyreduce (q);
+	endif
+
+	## remove common factors from numerator and denominator
+	x = polygcd (q, r);
+	if (length(x) != 1)
+	  q = deconv (q, x);
+	  r = deconv (r, x);
+	endif
+
+	## move all the gain into the numerator
+	q = q/r(1);
+	r = r/r(1);
+      endif
     else
-      ## derivative of p/a returns numerator and denominator
-      r = conv(a, a);
-      if numel(p) == 1
-	q = -p * polyderiv(a);
-      elseif numel(a) == 1
-	q = a * polyderiv(p);
-      else
-      	q = conv(polyderiv(p),a) - conv(p,polyderiv(a));
-      	q = polyreduce(q);
+      lp = numel (p);
+      if (lp == 1)
+	q = 0;
+	return;
+      elseif (lp == 0)
+	q = [];
+	return;
       endif
 
-      ## remove common factors from numerator and denominator
-      x = polygcd(q,r);
-      if length(x)!=1
-      	q=deconv(q,x);
-      	r=deconv(r,x);
-      endif
+      ## Force P to be a row vector.
+      p = p(:).';
 
-      ## move all the gain into the numerator
-      q=q/r(1);
-      r=r/r(1);
+      q = p(1:(lp-1)) .* [(lp-1):-1:1];
     endif
   else
-    lp = numel (p);
-    if (lp == 1)
-      q = 0;
-      return;
-    elseif (lp == 0)
-      q = [];
-      return;
-    end
-
-    ## Force P to be a row vector.
-    p = p(:).';
-
-    q = p (1:(lp-1)) .* [(lp-1):-1:1];
+    usage ("q = polyderiv (p) or q = polyderiv (b, a) or [q, r] = polyderiv (b, a)");
   endif
 
+
 endfunction
--- a/scripts/polynomial/polygcd.m	Wed Mar 16 19:51:46 2005 +0000
+++ b/scripts/polynomial/polygcd.m	Wed Mar 16 20:03:01 2005 +0000
@@ -37,29 +37,37 @@
 ## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv,
 ## residue, filter, polyval, and polyvalm}
 
-function x = polygcd(b,a,tol)
-  if (nargin<2 || nargin>3)
-    usage("x=polygcd(b,a [,tol])");
-  endif
-  if (nargin<3), tol=sqrt(eps); endif
-  if (length(a)==1 || length(b)==1)
-    if a==0, x=b;
-    elseif b==0, x=a;
-    else x=1;
+function x = polygcd (b, a, tol)
+
+  if (nargin == 2 || nargin == 3)
+    if (nargin == 2)
+      tol = sqrt (eps);
     endif
-    return;
+    if (length (a) == 1 || length (b) == 1)
+      if (a == 0)
+	x = b;
+      elseif (b == 0)
+	x = a;
+      else
+	x = 1;
+      endif
+    else
+      a /= a(1);
+      while (1)
+	[d, r] = deconv (b, a);
+	nz = find (abs (r) > tol);
+	if (isempty (nz))
+	  x = a;
+	  break;
+	else
+	  r = r(nz(1):length(r));
+	endif
+	b = a;
+	a /= r(1);
+      endwhile
+    endif
+  else
+    usage ("x = polygcd (b, a [,tol])");
   endif
-  a = a./a(1);
-  while (1)
-    [d, r] = deconv(b, a);
-    nz = find(abs(r)>tol);
-    if isempty(nz)
-      x = a; 
-      return;
-    else
-      r = r(nz(1):length(r));
-    endif
-    b = a;
-    a = r./r(1);
-  endwhile
-endfunction
\ No newline at end of file
+
+endfunction