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)