Mercurial > octave
view scripts/specfun/factor.m @ 21546:f7f97d7e9294
doc: Wrap m-file docstrings to 79 characters + newline (80 total).
* isrecording.m, soundsc.m, delaunay3.m, cell2mat.m, cumtrapz.m, del2.m,
inputParser.m, interp1.m, interp3.m, narginchk.m, profile.m,
validateattributes.m, delaunayn.m, tsearchn.m, voronoin.m, brighten.m,
cmunique.m, colorcube.m, imfinfo.m, imshow.m, edit.m, orderfields.m, run.m,
warning_ids.m, ode23.m, ode45.m, odeget.m, integrate_adaptive.m, kahan.m,
ode_struct_value_check.m, runge_kutta_23.m, fminunc.m, fsolve.m, fzero.m,
pkg.m, build.m, specular.m, view.m, bar.m, barh.m, contour3.m, isosurface.m,
line.m, pie.m, pie3.m, quiver3.m, scatter.m, scatter3.m, stem3.m, stemleaf.m,
surfl.m, tetramesh.m, isfigure.m, mkpp.m, pchip.m, residue.m, splinefit.m,
rmpref.m, unique.m, eigs.m, ilu.m, factor.m, factorial.m, gallery.m, hankel.m,
histc.m, ols.m, finv.m, fpdf.m, kruskal_wallis_test.m, weekday.m:
Wrap m-file docstrings to 79 characters + newline (80 total).
author | Rik <rik@octave.org> |
---|---|
date | Sun, 27 Mar 2016 15:50:01 -0700 |
parents | 516bb87ea72e |
children | b571fc85953f |
line wrap: on
line source
## Copyright (C) 2000-2015 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 {} {@var{pf} =} factor (@var{q}) ## @deftypefnx {} {[@var{pf}, @var{n}] =} factor (@var{q}) ## Return the prime factorization of @var{q}. ## ## The prime factorization is defined as @code{prod (@var{pf}) == @var{q}} ## where every element of @var{pf} is a prime number. If @code{@var{q} == 1}, ## return 1. ## ## With two output arguments, return the unique prime factors @var{pf} and ## their multiplicities. That is, ## @code{prod (@var{pf} .^ @var{n}) == @var{q}}. ## ## Implementation Note: The input @var{q} must be less than ## @code{flintmax} (9.0072e+15) in order to factor correctly. ## @seealso{gcd, lcm, isprime, primes} ## @end deftypefn ## Author: Paul Kienzle ## 2002-01-28 Paul Kienzle ## * remove recursion; only check existing primes for multiplicity > 1 ## * return multiplicity as suggested by Dirk Laurie ## * add error handling function [pf, n] = factor (q) if (nargin != 1) print_usage (); endif if (! isreal (q) || ! isscalar (q) || q != fix (q)) error ("factor: Q must be a real integer"); endif ## Special case of no primes less than sqrt(q). if (q < 4) pf = q; n = 1; return; endif q = double (q); # For the time being, calcs rely on double precision var. qorig = q; pf = []; ## There is at most one prime greater than sqrt(q), and if it exists, ## it has multiplicity 1, so no need to consider any factors greater ## than sqrt(q) directly. [If there were two factors p1, p2 > sqrt(q), ## then q >= p1*p2 > sqrt(q)*sqrt(q) == q. Contradiction.] p = primes (sqrt (q)); while (q > 1) ## Find prime factors in remaining q. p = p(rem (q, p) == 0); if (isempty (p)) ## Can't be reduced further, so q must itself be a prime. p = q; endif pf = [pf, p]; ## Reduce q. q /= prod (p); endwhile pf = sort (pf); ## Verify algorithm was succesful q = prod (pf); if (q != qorig) error ("factor: Q too large to factor"); elseif (q >= flintmax ()) warning ("factor: Q too large. Answer is unreliable"); endif ## Determine muliplicity. if (nargout > 1) idx = find ([0, pf] != [pf, 0]); pf = pf(idx(1:length (idx)-1)); n = diff (idx); endif endfunction %!assert (factor (1), 1) %!test %! for i = 2:20 %! pf = factor (i); %! assert (prod (pf), i); %! assert (all (isprime (pf))); %! [pf, n] = factor (i); %! assert (prod (pf.^n), i); %! assert (all ([0,pf] != [pf,0])); %! endfor ## Test input validation %!error factor () %!error factor (1,2) %!error <Q must be a real integer> factor (6i) %!error <Q must be a real integer> factor ([1,2]) %!error <Q must be a real integer> factor (1.5)