# HG changeset patch # User jwe # Date 1111002706 0 # Node ID 5ed60b8b1ac4e6dc06bfde949a5707e5884d46f3 # Parent 32c569794216f27ba66c10daf5b5ae98064a2f97 [project @ 2005-03-16 19:51:39 by jwe] diff -r 32c569794216 -r 5ed60b8b1ac4 scripts/polynomial/polyder.m --- a/scripts/polynomial/polyder.m Wed Mar 16 18:54:52 2005 +0000 +++ b/scripts/polynomial/polyder.m Wed Mar 16 19:51:46 2005 +0000 @@ -19,17 +19,27 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} polyder (@var{c}) +## @deftypefnx {Function File} {[@var{q}] =} polyder (@var{b}, @var{a}) +## @deftypefnx {Function File} {[@var{q}, @var{r}] =} polyder (@var{b}, @var{a}) ## See polyderiv. ## @end deftypefn -## Author: jwe +## Author: John W. Eaton +## Paul Kienzle +## handle b/a and b*a -function q = polyder (p) +function [q, r] = polyder (p, a) if (nargin == 1) q = polyderiv (p); + elseif (nargin==2) + if (nargout==2) + [q, r] = polyderiv (p,a); + else + q = polyderiv (p,a); + endif else - usage ("polyder (vector)"); + usage ("q=polyder(p) or q=polyder(b,a) or [q, r]=polyder(b,a)"); endif endfunction diff -r 32c569794216 -r 5ed60b8b1ac4 scripts/polynomial/polyderiv.m --- a/scripts/polynomial/polyderiv.m Wed Mar 16 18:54:52 2005 +0000 +++ b/scripts/polynomial/polyderiv.m Wed Mar 16 19:51:46 2005 +0000 @@ -19,39 +19,77 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} polyderiv (@var{c}) +## @deftypefnx {Function File} {[@var{q}] =} polyder (@var{b}, @var{a}) +## @deftypefnx {Function File} {[@var{q}, @var{r}] =} polyder (@var{b}, @var{a}) ## Return the coefficients of the derivative of the polynomial whose -## coefficients are given by vector @var{c}. +## coefficients are given by vector @var{c}. If a pair of polynomials +## is given @var{b} and @var{a}, the derivative of the product is +## returned in @var{q}, or the quotient numerator in @var{q} and the +## quotient denominator in @var{r}. ## @end deftypefn -## ## @seealso{poly, polyinteg, polyreduce, roots, conv, deconv, residue, -## filter, polyval, and polyvalm} +## filter, polygcd, polyval, and polyvalm} ## Author: Tony Richardson ## Created: June 1994 ## Adapted-By: jwe - -function q = polyderiv (p) +## Paul Kienzle +## handle b/a and b*a - if (nargin != 1) - usage ("polyderiv (vector)"); +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 - lp = numel (p); - if (lp == 1) - q = 0; - return; - elseif (lp == 0) - q = []; - return; - end + 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 - ## Force P to be a row vector. - p = p(:).'; + ## remove common factors from numerator and denominator + x = polygcd(q,r); + if length(x)!=1 + q=deconv(q,x); + r=deconv(r,x); + endif - q = p(1:(lp-1)) .* [(lp-1):-1:1]; + ## move all the gain into the numerator + q=q/r(1); + r=r/r(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]; + endif endfunction diff -r 32c569794216 -r 5ed60b8b1ac4 scripts/polynomial/polygcd.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/polygcd.m Wed Mar 16 19:51:46 2005 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 2000 Paul Kienzle +## +## This program 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 of the License, or +## (at your option) any later version. +## +## This program 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 this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{q}]} polygcd (@var{b}, @var{a}, @var{tol}) +## +## Find greatest common divisor of two polynomials. This is equivalent +## to the polynomial found by multiplying together all the common roots. +## Together with deconv, you can reduce a ratio of two polynomials. +## Tolerance defaults to +## @example +## sqrt(eps). +## @end example +## Note that this is an unstable +## algorithm, so don't try it on large polynomials. +## +## Example +## @example +## polygcd(poly(1:8),poly(3:12)) - poly(3:8) +## deconv(poly(1:8),polygcd(poly(1:8),poly(3:12))) - poly(1:2) +## @end example +## @end deftypefn +## +## @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; + endif + return; + 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