Mercurial > octave-nkf
changeset 19131:e7bffcea619f
factor.m: Overhaul function.
* factor.m: Rewrite docstring. Use same variable names between documentation
and code. Add input validation that Q is a real number. Shorten error
messages. Add input validation tests for Q being real.
author | Rik <rik@octave.org> |
---|---|
date | Fri, 19 Sep 2014 14:21:46 -0700 |
parents | 90dd85369c2e |
children | 4591a1238ee0 |
files | scripts/specfun/factor.m |
diffstat | 1 files changed, 35 insertions(+), 33 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/factor.m Fri Sep 19 12:51:12 2014 -0700 +++ b/scripts/specfun/factor.m Fri Sep 19 14:21:46 2014 -0700 @@ -17,20 +17,20 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{p} =} factor (@var{q}) -## @deftypefnx {Function File} {[@var{p}, @var{n}] =} factor (@var{q}) -## -## Return the prime factorization of @var{q}. That is, -## @code{prod (@var{p}) == @var{q}} and every element of @var{p} is a prime -## number. If @code{@var{q} == 1}, return 1. +## @deftypefn {Function File} {@var{pf} =} factor (@var{q}) +## @deftypefnx {Function File} {[@var{pf}, @var{n}] =} factor (@var{q}) +## Return the prime factorization of @var{q}. ## -## With two output arguments, return the unique primes @var{p} and -## their multiplicities. That is, @code{prod (@var{p} .^ @var{n}) == -## @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 not be greater than +## Implementation Note: The input @var{q} must be less than than ## @code{bitmax} (9.0072e+15) in order to factor correctly. -## @seealso{gcd, lcm, isprime} +## @seealso{gcd, lcm, isprime, primes} ## @end deftypefn ## Author: Paul Kienzle @@ -40,29 +40,29 @@ ## * return multiplicity as suggested by Dirk Laurie ## * add error handling -function [x, n] = factor (q) +function [pf, n] = factor (q) - if (nargin < 1) + if (nargin != 1) print_usage (); endif - if (! isscalar (q) || q != fix (q)) - error ("factor: Q must be a scalar integer"); + 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) - x = q; + pf = q; n = 1; return; endif q = double (q); # For the time being, calcs rely on double precision var. qorig = q; - x = []; + 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), + ## 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) @@ -72,24 +72,24 @@ ## Can't be reduced further, so q must itself be a prime. p = q; endif - x = [x, p]; + pf = [pf, p]; ## Reduce q. q /= prod (p); endwhile - x = sort (x); + pf = sort (pf); ## Verify algorithm was succesful - q = prod (x); + q = prod (pf); if (q != qorig) - error ("factor: Input Q too large to factor"); + error ("factor: Q too large to factor"); elseif (q > bitmax) - warning ("factor: Input Q too large. Answer is unreliable"); + warning ("factor: Q too large. Answer is unreliable"); endif ## Determine muliplicity. if (nargout > 1) - idx = find ([0, x] != [x, 0]); - x = x(idx(1:length (idx)-1)); + idx = find ([0, pf] != [pf, 0]); + pf = pf(idx(1:length (idx)-1)); n = diff (idx); endif @@ -99,16 +99,18 @@ %!assert (factor (1), 1) %!test %! for i = 2:20 -%! p = factor (i); -%! assert (prod (p), i); -%! assert (all (isprime (p))); -%! [p,n] = factor (i); -%! assert (prod (p.^n), i); -%! assert (all ([0,p] != [p,0])); +%! 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 <Q must be a scalar integer> factor ([1,2]) -%!error <Q must be a scalar integer> factor (1.5) +%!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)