Mercurial > hg > octave-nkf
changeset 17716:e48f5a52e838
factor.m: Warn when the input is too large to calculate reliable answer.
* scripts/specfun/factor.m: After factorization, test that product of
calculated factors equals original input. Add note to docstring about
limitation to inputs less than bitmax(). Add %!error input validation.
author | Doug Stewart <doug.dastew@gmail.com> |
---|---|
date | Fri, 07 Dec 2012 10:03:33 -0500 |
parents | 8dd280b64de1 |
children | 76f448d8089d |
files | scripts/specfun/factor.m |
diffstat | 1 files changed, 22 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/factor.m +++ b/scripts/specfun/factor.m @@ -20,14 +20,17 @@ ## @deftypefn {Function File} {@var{p} =} factor (@var{q}) ## @deftypefnx {Function File} {[@var{p}, @var{n}] =} factor (@var{q}) ## -## Return prime factorization of @var{q}. That is, +## 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. ## ## With two output arguments, return the unique primes @var{p} and ## their multiplicities. That is, @code{prod (@var{p} .^ @var{n}) == ## @var{q}}. -## @seealso{gcd, lcm} +## +## Implementation Note: The input @var{q} must not be greater than +## @code{bitmax} (9.0072e+15) in order to factor correctly. +## @seealso{gcd, lcm, isprime} ## @end deftypefn ## Author: Paul Kienzle @@ -54,6 +57,8 @@ return; endif + q = double (q); # For the time being, calcs rely on double precision var. + qorig = q; x = []; ## 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 @@ -62,17 +67,25 @@ p = primes (sqrt (q)); while (q > 1) ## Find prime factors in remaining q. - p = p (rem (q, p) == 0); + p = p(rem (q, p) == 0); if (isempty (p)) ## Can't be reduced further, so q must itself be a prime. p = q; endif x = [x, p]; ## Reduce q. - q = q / prod (p); + q /= prod (p); endwhile x = sort (x); + ## Verify algorithm was succesful + q = prod (x); + if (q != qorig) + error ("factor: Input Q too large to factor"); + elseif (q > bitmax) + warning ("factor: Input Q too large. Answer is unreliable"); + endif + ## Determine muliplicity. if (nargout > 1) idx = find ([0, x] != [x, 0]); @@ -94,3 +107,8 @@ %! assert (all ([0,p] != [p,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) +