Mercurial > octave-nkf

comparison scripts/polynomial/polygcd.m @ 5216:5ed60b8b1ac4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-03-16 19:51:39 by jwe]
author jwe
date Wed, 16 Mar 2005 19:51:46 +0000
parents
children e88886a6934d
comparison
equal deleted inserted replaced
5215:32c569794216 5216:5ed60b8b1ac4
1 ## Copyright (C) 2000 Paul Kienzle
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 ## -*- texinfo -*-
18 ## @deftypefn {Function File} {[@var{q}]} polygcd (@var{b}, @var{a}, @var{tol})
19 ##
20 ## Find greatest common divisor of two polynomials. This is equivalent
21 ## to the polynomial found by multiplying together all the common roots.
22 ## Together with deconv, you can reduce a ratio of two polynomials.
23 ## Tolerance defaults to
24 ## @example
25 ## sqrt(eps).
26 ## @end example
27 ## Note that this is an unstable
28 ## algorithm, so don't try it on large polynomials.
29 ##
30 ## Example
31 ## @example
32 ## polygcd(poly(1:8),poly(3:12)) - poly(3:8)
33 ## deconv(poly(1:8),polygcd(poly(1:8),poly(3:12))) - poly(1:2)
34 ## @end example
35 ## @end deftypefn
36 ##
37 ## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv,
38 ## residue, filter, polyval, and polyvalm}
39
40 function x = polygcd(b,a,tol)
41 if (nargin<2 || nargin>3)
42 usage("x=polygcd(b,a [,tol])");
43 endif
44 if (nargin<3), tol=sqrt(eps); endif
45 if (length(a)==1 || length(b)==1)
46 if a==0, x=b;
47 elseif b==0, x=a;
48 else x=1;
49 endif
50 return;
51 endif
52 a = a./a(1);
53 while (1)
54 [d, r] = deconv(b, a);
55 nz = find(abs(r)>tol);
56 if isempty(nz)
57 x = a;
58 return;
59 else
60 r = r(nz(1):length(r));
61 endif
62 b = a;
63 a = r./r(1);
64 endwhile
65 endfunction