Mercurial > octave-nkf

changeset 17716:e48f5a52e838

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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	Mon Oct 21 22:57:34 2013 +0100
+++ b/scripts/specfun/factor.m	Fri Dec 07 10:03:33 2012 -0500
@@ -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)
+