Mercurial > octave-nkf

comparison scripts/polynomial/polygcd.m @ 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 4c8a2e4e0717
comparison
equal deleted inserted replaced
5216:5ed60b8b1ac4 5217:e88886a6934d
35 ## @end deftypefn 35 ## @end deftypefn
36 ## 36 ##
37 ## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv, 37 ## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv,
38 ## residue, filter, polyval, and polyvalm} 38 ## residue, filter, polyval, and polyvalm}
39 39
40 function x = polygcd(b,a,tol) 40 function x = polygcd (b, a, tol)
41 if (nargin<2 || nargin>3) 41
42 usage("x=polygcd(b,a [,tol])"); 42 if (nargin == 2 || nargin == 3)
43 if (nargin == 2)
44 tol = sqrt (eps);
45 endif
46 if (length (a) == 1 || length (b) == 1)
47 if (a == 0)
48 x = b;
49 elseif (b == 0)
50 x = a;
51 else
52 x = 1;
53 endif
54 else
55 a /= a(1);
56 while (1)
57 [d, r] = deconv (b, a);
58 nz = find (abs (r) > tol);
59 if (isempty (nz))
60 x = a;
61 break;
62 else
63 r = r(nz(1):length(r));
64 endif
65 b = a;
66 a /= r(1);
67 endwhile
68 endif
69 else
70 usage ("x = polygcd (b, a [,tol])");
43 endif 71 endif
44 if (nargin<3), tol=sqrt(eps); endif 72
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 73 endfunction