Mercurial > octave-nkf

comparison scripts/polynomial/polygcd.m @ 10549:95c3e38098bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Untabify .m scripts
author Rik <code@nomad.inbox5.com>
date Fri, 23 Apr 2010 11:28:50 -0700
parents f6e0404421f4
children 693e22af08ae
comparison
equal deleted inserted replaced
10548:479536c5bb10 10549:95c3e38098bf
46 function x = polygcd (b, a, tol) 46 function x = polygcd (b, a, tol)
47 47
48 if (nargin == 2 || nargin == 3) 48 if (nargin == 2 || nargin == 3)
49 if (nargin == 2) 49 if (nargin == 2)
50 if (isa (a, "single") || isa (b, "single")) 50 if (isa (a, "single") || isa (b, "single"))
51 tol = sqrt (eps ("single")); 51 tol = sqrt (eps ("single"));
52 else 52 else
53 tol = sqrt (eps); 53 tol = sqrt (eps);
54 endif 54 endif
55 endif 55 endif
56 if (length (a) == 1 || length (b) == 1) 56 if (length (a) == 1 || length (b) == 1)
57 if (a == 0) 57 if (a == 0)
58 x = b; 58 x = b;
59 elseif (b == 0) 59 elseif (b == 0)
60 x = a; 60 x = a;
61 else 61 else
62 x = 1; 62 x = 1;
63 endif 63 endif
64 else 64 else
65 a /= a(1); 65 a /= a(1);
66 while (1) 66 while (1)
67 [d, r] = deconv (b, a); 67 [d, r] = deconv (b, a);
68 nz = find (abs (r) > tol); 68 nz = find (abs (r) > tol);
69 if (isempty (nz)) 69 if (isempty (nz))
70 x = a; 70 x = a;
71 break; 71 break;
72 else 72 else
73 r = r(nz(1):length(r)); 73 r = r(nz(1):length(r));
74 endif 74 endif
75 b = a; 75 b = a;
76 a = r / r(1); 76 a = r / r(1);
77 endwhile 77 endwhile
78 endif 78 endif
79 else 79 else
80 print_usage (); 80 print_usage ();
81 endif 81 endif