Mercurial > octave
view scripts/polynomial/polygcd.m @ 14327:4d917a6a858b stable
doc: Use Octave coding conventions in @example blocks of docstrings.
* accumarray.m, accumdim.m, bar.m, base2dec.m, bincoeff.m, bitcmp.m, bitset.m,
celldisp.m, chop.m, clabel.m, cloglog.m, colon.m, compass.m, computer.m,
contour3.m, contourc.m, corr.m, cstrcat.m, ctime.m, cylinder.m, date.m,
dec2base.m, demo.m, dir.m, dlmwrite.m, expm.m, ezcontourf.m, ezcontour.m,
ezmeshc.m, ezmesh.m, ezplot.m, ezsurfc.m, ezsurf.m, feather.m, findobj.m,
flipdim.m, fplot.m, genvarname.m, getfield.m, hankel.m, hilb.m, hist.m,
idivide.m, index.m, int2str.m, interp1.m, is_leap_year.m, ismember.m,
isocolors.m, isonormals.m, isosurface.m, kurtosis.m, legendre.m, linkprop.m,
logit.m, logm.m, __makeinfo__.m, __marching_cube__.m, median.m, mkoctfile.m,
moment.m, mpoles.m, orderfields.m, pcg.m, pcr.m, plot3.m, plotmatrix.m,
polyaffine.m, polygcd.m, poly.m, polyout.m, print.m, qp.m, quadgk.m, qzhess.m,
randi.m, rat.m, refreshdata.m, residue.m, rose.m, rot90.m, saveas.m, saveobj.m,
shiftdim.m, skewness.m, spaugment.m, spdiags.m, sqp.m, stem.m, str2num.m,
strcat.m, strjust.m, strread.m, strsplit.m, structfun.m, subplot.m,
subsindex.m, substruct.m, surfl.m, surfnorm.m, svds.m, uimenu.m, union.m,
voronoi.m, warning_ids.m, wblpdf.m: Use Octave coding conventions in
@example blocks of docstrings.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sat, 04 Feb 2012 22:12:50 -0800 |
parents | 72c96de7a403 |
children | f3d52523cde1 |
line wrap: on
line source
## Copyright (C) 2000-2012 Paul Kienzle ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{q} =} polygcd (@var{b}, @var{a}) ## @deftypefnx {Function File} {@var{q} =} polygcd (@var{b}, @var{a}, @var{tol}) ## ## Find the greatest common divisor of two polynomials. This is equivalent ## to the polynomial found by multiplying together all the common roots. ## Together with deconv, you can reduce a ratio of two polynomials. ## The tolerance @var{tol} defaults to @code{sqrt(eps)}. ## ## @strong{Caution:} This is a numerically unstable algorithm and should not ## be used on large polynomials. ## ## Example code: ## ## @example ## @group ## polygcd (poly (1:8), poly (3:12)) - poly (3:8) ## @result{} [ 0, 0, 0, 0, 0, 0, 0 ] ## deconv (poly (1:8), polygcd (poly (1:8), poly (3:12))) - poly(1:2) ## @result{} [ 0, 0, 0 ] ## @end group ## @end example ## @seealso{poly, roots, conv, deconv, residue} ## @end deftypefn function x = polygcd (b, a, tol) if (nargin == 2 || nargin == 3) if (nargin == 2) if (isa (a, "single") || isa (b, "single")) tol = sqrt (eps ("single")); else tol = sqrt (eps); endif endif if (length (a) == 1 || length (b) == 1) if (a == 0) x = b; elseif (b == 0) x = a; else x = 1; endif else a /= a(1); while (1) [d, r] = deconv (b, a); nz = find (abs (r) > tol); if (isempty (nz)) x = a; break; else r = r(nz(1):length(r)); endif b = a; a = r / r(1); endwhile endif else print_usage (); endif endfunction %!test %! poly1 = [1 6 11 6]; % (x+1)(x+2)(x+3) %! poly2 = [1 3 2]; % (x+1)(x+2) %! poly3 = polygcd (poly1, poly2); %! assert (poly3, poly2, sqrt (eps)) %!test %! assert (polygcd (poly(1:8), poly(3:12)), poly(3:8), sqrt (eps)) %!test %! assert (deconv (poly(1:8), polygcd (poly(1:8), poly(3:12))), poly(1:2), sqrt (eps)) %!test %! for ii=1:10 %! p = (unique (randn (10, 1)) * 10).'; %! p1 = p(3:end); %! p2 = p(1:end-2); %! assert (polygcd (poly (-p1), poly (-p2)), poly (- intersect (p1, p2)), sqrt (eps)) %! endfor