← Back to branch summary

~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to scripts/polynomial/polygcd.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
Import upstream version 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/polynomial/polygcd.m
 
1
## Copyright (C) 2000, 2005, 2006, 2007 Paul Kienzle
 
2
##
 
3
## This file is part of Octave.
 
4
##
 
5
## Octave is free software; you can redistribute it and/or modify it
 
6
## under the terms of the GNU General Public License as published by
 
7
## the Free Software Foundation; either version 3 of the License, or (at
 
8
## your option) any later version.
 
9
##
 
10
## Octave is distributed in the hope that it will be useful, but
 
11
## WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
13
## General Public License for more details.
 
14
##
 
15
## You should have received a copy of the GNU General Public License
 
16
## along with Octave; see the file COPYING.  If not, see
 
17
## <http://www.gnu.org/licenses/>.
 
18
 
 
19
## -*- texinfo -*-
 
20
## @deftypefn {Function File} {@var{q} =} polygcd (@var{b}, @var{a}, @var{tol})
 
21
##
 
22
## Find greatest common divisor of two polynomials.  This is equivalent
 
23
## to the polynomial found by multiplying together all the common roots.
 
24
## Together with deconv, you can reduce a ratio of two polynomials.
 
25
## Tolerance defaults to 
 
26
## @example 
 
27
## sqrt(eps).
 
28
## @end example
 
29
##  Note that this is an unstable
 
30
## algorithm, so don't try it on large polynomials.
 
31
##
 
32
## Example
 
33
## @example
 
34
## polygcd (poly(1:8), poly(3:12)) - poly(3:8)
 
35
## @result{} [ 0, 0, 0, 0, 0, 0, 0 ]
 
36
## deconv (poly(1:8), polygcd (poly(1:8), poly(3:12))) ...
 
37
##   - poly(1:2)
 
38
## @result{} [ 0, 0, 0 ]
 
39
## @end example
 
40
## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv,
 
41
## residue, filter, polyval, and polyvalm}
 
42
## @end deftypefn
 
43
 
 
44
function x = polygcd (b, a, tol)
 
45
 
 
46
  if (nargin == 2 || nargin == 3)
 
47
    if (nargin == 2)
 
48
      tol = sqrt (eps);
 
49
    endif
 
50
    if (length (a) == 1 || length (b) == 1)
 
51
      if (a == 0)
 
52
        x = b;
 
53
      elseif (b == 0)
 
54
        x = a;
 
55
      else
 
56
        x = 1;
 
57
      endif
 
58
    else
 
59
      a /= a(1);
 
60
      while (1)
 
61
        [d, r] = deconv (b, a);
 
62
        nz = find (abs (r) > tol);
 
63
        if (isempty (nz))
 
64
          x = a;
 
65
          break;
 
66
        else
 
67
          r = r(nz(1):length(r));
 
68
        endif
 
69
        b = a;
 
70
        a = r / r(1);
 
71
      endwhile
 
72
    endif
 
73
  else
 
74
    print_usage ();
 
75
  endif
 
76
 
 
77
endfunction