Mercurial > hg > octave-lyh
annotate scripts/polynomial/polygcd.m @ 8286:6f2d95255911
fix @seealso references to point to existing anchors
| author | Thorsten Meyer <thorsten.meyier@gmx.de> |
|---|---|
| date | Wed, 29 Oct 2008 17:52:54 -0400 |
| parents | df9519e9990c |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 2000, 2005, 2006, 2007 Paul Kienzle |
| 5216 | 2 ## |
| 7016 | 3 ## This file is part of Octave. |
| 5216 | 4 ## |
| 7016 | 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. | |
| 5216 | 14 ## |
| 15 ## You should have received a copy of the GNU General Public License | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5216 | 18 |
| 19 ## -*- texinfo -*- | |
| 5382 | 20 ## @deftypefn {Function File} {@var{q} =} polygcd (@var{b}, @var{a}, @var{tol}) |
| 5216 | 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 | |
| 5382 | 34 ## polygcd (poly(1:8), poly(3:12)) - poly(3:8) |
| 6850 | 35 ## @result{} [ 0, 0, 0, 0, 0, 0, 0 ] |
| 7031 | 36 ## deconv (poly(1:8), polygcd (poly(1:8), poly(3:12))) ... |
| 37 ## - poly(1:2) | |
| 6850 | 38 ## @result{} [ 0, 0, 0 ] |
| 5216 | 39 ## @end example |
| 40 ## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv, | |
|
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
7795
diff
changeset
|
41 ## residue, filter, polyval, polyvalm} |
| 5642 | 42 ## @end deftypefn |
| 5216 | 43 |
| 5217 | 44 function x = polygcd (b, a, tol) |
| 45 | |
| 46 if (nargin == 2 || nargin == 3) | |
| 47 if (nargin == 2) | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
48 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
49 tol = sqrt (eps ("single")); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
50 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
51 tol = sqrt (eps); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
52 endif |
| 5216 | 53 endif |
| 5217 | 54 if (length (a) == 1 || length (b) == 1) |
| 55 if (a == 0) | |
| 56 x = b; | |
| 57 elseif (b == 0) | |
| 58 x = a; | |
| 59 else | |
| 60 x = 1; | |
| 61 endif | |
| 62 else | |
| 63 a /= a(1); | |
| 64 while (1) | |
| 65 [d, r] = deconv (b, a); | |
| 66 nz = find (abs (r) > tol); | |
| 67 if (isempty (nz)) | |
| 68 x = a; | |
| 69 break; | |
| 70 else | |
| 71 r = r(nz(1):length(r)); | |
| 72 endif | |
| 73 b = a; | |
| 5382 | 74 a = r / r(1); |
| 5217 | 75 endwhile |
| 76 endif | |
| 77 else | |
| 6046 | 78 print_usage (); |
| 5216 | 79 endif |
| 5217 | 80 |
| 81 endfunction |