Mercurial > octave-nkf
annotate scripts/polynomial/polygcd.m @ 20595:c1a6c31ac29a
eliminate more simple uses of error_state
* ov-classdef.cc: Eliminate simple uses of error_state.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 06 Oct 2015 00:20:02 -0400 |
parents | f1d0f506ee78 |
children |
rev | line source |
---|---|
19731
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19258
diff
changeset
|
1 ## Copyright (C) 2000-2015 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 -*- | |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
20 ## @deftypefn {Function File} {@var{q} =} polygcd (@var{b}, @var{a}) |
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
21 ## @deftypefnx {Function File} {@var{q} =} polygcd (@var{b}, @var{a}, @var{tol}) |
5216 | 22 ## |
20200
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
23 ## Find the greatest common divisor of two polynomials. |
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
24 ## |
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
25 ## This is equivalent to the polynomial found by multiplying together all the |
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
26 ## common roots. Together with deconv, you can reduce a ratio of two |
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
27 ## polynomials. |
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
28 ## |
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
29 ## The tolerance @var{tol} defaults to @code{sqrt (eps)}. |
5216 | 30 ## |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
31 ## @strong{Caution:} This is a numerically unstable algorithm and should not |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
32 ## be used on large polynomials. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
33 ## |
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
34 ## Example code: |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
35 ## |
5216 | 36 ## @example |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
37 ## @group |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
38 ## polygcd (poly (1:8), poly (3:12)) - poly (3:8) |
6850 | 39 ## @result{} [ 0, 0, 0, 0, 0, 0, 0 ] |
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
40 ## deconv (poly (1:8), polygcd (poly (1:8), poly (3:12))) - poly (1:2) |
6850 | 41 ## @result{} [ 0, 0, 0 ] |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
42 ## @end group |
5216 | 43 ## @end example |
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
44 ## @seealso{poly, roots, conv, deconv, residue} |
5642 | 45 ## @end deftypefn |
5216 | 46 |
5217 | 47 function x = polygcd (b, a, tol) |
48 | |
49 if (nargin == 2 || nargin == 3) | |
50 if (nargin == 2) | |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
51 if (isa (a, "single") || isa (b, "single")) |
10549 | 52 tol = sqrt (eps ("single")); |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
53 else |
10549 | 54 tol = sqrt (eps); |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7031
diff
changeset
|
55 endif |
5216 | 56 endif |
5217 | 57 if (length (a) == 1 || length (b) == 1) |
58 if (a == 0) | |
10549 | 59 x = b; |
5217 | 60 elseif (b == 0) |
10549 | 61 x = a; |
5217 | 62 else |
10549 | 63 x = 1; |
5217 | 64 endif |
65 else | |
66 a /= a(1); | |
67 while (1) | |
10549 | 68 [d, r] = deconv (b, a); |
69 nz = find (abs (r) > tol); | |
70 if (isempty (nz)) | |
71 x = a; | |
72 break; | |
73 else | |
74 r = r(nz(1):length(r)); | |
75 endif | |
76 b = a; | |
77 a = r / r(1); | |
5217 | 78 endwhile |
79 endif | |
80 else | |
6046 | 81 print_usage (); |
5216 | 82 endif |
5217 | 83 |
84 endfunction | |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
85 |
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
86 |
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
87 %!test |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
88 %! poly1 = [1 6 11 6]; # (x+1)(x+2)(x+3); |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
89 %! poly2 = [1 3 2]; # (x+1)(x+2); |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
90 %! poly3 = polygcd (poly1, poly2); |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
91 %! assert (poly3, poly2, sqrt (eps)); |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
92 |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
93 %!assert (polygcd (poly (1:8), poly (3:12)), poly (3:8), sqrt (eps)) |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
94 %!assert (deconv (poly (1:8), polygcd (poly (1:8), poly (3:12))), poly (1:2), sqrt (eps)) |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
95 |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13045
diff
changeset
|
96 %!test |
19258
5d3111977623
Change a polygcd test to one that is more numerically accurate (bug #42742).
Daniel J Sebald <daniel.sebald@ieee.org>
parents:
17744
diff
changeset
|
97 %! for ii=1:100 |
5d3111977623
Change a polygcd test to one that is more numerically accurate (bug #42742).
Daniel J Sebald <daniel.sebald@ieee.org>
parents:
17744
diff
changeset
|
98 %! ## Exhibits numerical problems for multipliers of ~4 and greater. |
5d3111977623
Change a polygcd test to one that is more numerically accurate (bug #42742).
Daniel J Sebald <daniel.sebald@ieee.org>
parents:
17744
diff
changeset
|
99 %! p = (unique (randn (10, 1)) * 3).'; |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
100 %! p1 = p(3:end); |
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
101 %! p2 = p(1:end-2); |
19258
5d3111977623
Change a polygcd test to one that is more numerically accurate (bug #42742).
Daniel J Sebald <daniel.sebald@ieee.org>
parents:
17744
diff
changeset
|
102 %! assert (polygcd (poly (-p1), poly (-p2)), |
5d3111977623
Change a polygcd test to one that is more numerically accurate (bug #42742).
Daniel J Sebald <daniel.sebald@ieee.org>
parents:
17744
diff
changeset
|
103 %! poly (- intersect (p1, p2)), sqrt (eps)); |
13045
92cb87addf25
codesprint: Add tests to polygcd
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11523
diff
changeset
|
104 %! endfor |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
105 |