Mercurial > hg > octave-nkf

comparison scripts/specfun/factor.m @ 19289:e7bffcea619f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 d63878346099
children 0f9c5a15c8fa
comparison
equal deleted inserted replaced
19288:90dd85369c2e 19289:e7bffcea619f
15 ## You should have received a copy of the GNU General Public License 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {@var{p} =} factor (@var{q}) 20 ## @deftypefn {Function File} {@var{pf} =} factor (@var{q})
21 ## @deftypefnx {Function File} {[@var{p}, @var{n}] =} factor (@var{q}) 21 ## @deftypefnx {Function File} {[@var{pf}, @var{n}] =} factor (@var{q})
22 ## Return the prime factorization of @var{q}.
22 ## 23 ##
23 ## Return the prime factorization of @var{q}. That is, 24 ## The prime factorization is defined as @code{prod (@var{pf}) == @var{q}}
24 ## @code{prod (@var{p}) == @var{q}} and every element of @var{p} is a prime 25 ## where every element of @var{pf} is a prime number. If @code{@var{q} == 1},
25 ## number. If @code{@var{q} == 1}, return 1. 26 ## return 1.
26 ## 27 ##
27 ## With two output arguments, return the unique primes @var{p} and 28 ## With two output arguments, return the unique prime factors @var{pf} and
28 ## their multiplicities. That is, @code{prod (@var{p} .^ @var{n}) == 29 ## their multiplicities. That is, @code{prod (@var{pf} .^ @var{n}) == @var{q}}.
29 ## @var{q}}.
30 ## 30 ##
31 ## Implementation Note: The input @var{q} must not be greater than 31 ## Implementation Note: The input @var{q} must be less than than
32 ## @code{bitmax} (9.0072e+15) in order to factor correctly. 32 ## @code{bitmax} (9.0072e+15) in order to factor correctly.
33 ## @seealso{gcd, lcm, isprime} 33 ## @seealso{gcd, lcm, isprime, primes}
34 ## @end deftypefn 34 ## @end deftypefn
35 35
36 ## Author: Paul Kienzle 36 ## Author: Paul Kienzle
37 37
38 ## 2002-01-28 Paul Kienzle 38 ## 2002-01-28 Paul Kienzle
39 ## * remove recursion; only check existing primes for multiplicity > 1 39 ## * remove recursion; only check existing primes for multiplicity > 1
40 ## * return multiplicity as suggested by Dirk Laurie 40 ## * return multiplicity as suggested by Dirk Laurie
41 ## * add error handling 41 ## * add error handling
42 42
43 function [x, n] = factor (q) 43 function [pf, n] = factor (q)
44 44
45 if (nargin < 1) 45 if (nargin != 1)
46 print_usage (); 46 print_usage ();
47 endif 47 endif
48 48
49 if (! isscalar (q) || q != fix (q)) 49 if (! isreal (q) || ! isscalar (q) || q != fix (q))
50 error ("factor: Q must be a scalar integer"); 50 error ("factor: Q must be a real integer");
51 endif 51 endif
52 52
53 ## Special case of no primes less than sqrt(q). 53 ## Special case of no primes less than sqrt(q).
54 if (q < 4) 54 if (q < 4)
55 x = q; 55 pf = q;
56 n = 1; 56 n = 1;
57 return; 57 return;
58 endif 58 endif
59 59
60 q = double (q); # For the time being, calcs rely on double precision var. 60 q = double (q); # For the time being, calcs rely on double precision var.
61 qorig = q; 61 qorig = q;
62 x = []; 62 pf = [];
63 ## There is at most one prime greater than sqrt(q), and if it exists, 63 ## There is at most one prime greater than sqrt(q), and if it exists,
64 ## it has multiplicity 1, so no need to consider any factors greater 64 ## it has multiplicity 1, so no need to consider any factors greater
65 ## than sqrt(q) directly. [If there were two factors p1, p2 > sqrt(q), 65 ## than sqrt(q) directly. [If there were two factors p1, p2 > sqrt(q),
66 ## then q >= p1*p2 > sqrt(q)*sqrt(q) == q. Contradiction.] 66 ## then q >= p1*p2 > sqrt(q)*sqrt(q) == q. Contradiction.]
67 p = primes (sqrt (q)); 67 p = primes (sqrt (q));
68 while (q > 1) 68 while (q > 1)
69 ## Find prime factors in remaining q. 69 ## Find prime factors in remaining q.
70 p = p(rem (q, p) == 0); 70 p = p(rem (q, p) == 0);
71 if (isempty (p)) 71 if (isempty (p))
72 ## Can't be reduced further, so q must itself be a prime. 72 ## Can't be reduced further, so q must itself be a prime.
73 p = q; 73 p = q;
74 endif 74 endif
75 x = [x, p]; 75 pf = [pf, p];
76 ## Reduce q. 76 ## Reduce q.
77 q /= prod (p); 77 q /= prod (p);
78 endwhile 78 endwhile
79 x = sort (x); 79 pf = sort (pf);
80 80
81 ## Verify algorithm was succesful 81 ## Verify algorithm was succesful
82 q = prod (x); 82 q = prod (pf);
83 if (q != qorig) 83 if (q != qorig)
84 error ("factor: Input Q too large to factor"); 84 error ("factor: Q too large to factor");
85 elseif (q > bitmax) 85 elseif (q > bitmax)
86 warning ("factor: Input Q too large. Answer is unreliable"); 86 warning ("factor: Q too large. Answer is unreliable");
87 endif 87 endif
88 88
89 ## Determine muliplicity. 89 ## Determine muliplicity.
90 if (nargout > 1) 90 if (nargout > 1)
91 idx = find ([0, x] != [x, 0]); 91 idx = find ([0, pf] != [pf, 0]);
92 x = x(idx(1:length (idx)-1)); 92 pf = pf(idx(1:length (idx)-1));
93 n = diff (idx); 93 n = diff (idx);
94 endif 94 endif
95 95
96 endfunction 96 endfunction
97 97
98 98
99 %!assert (factor (1), 1) 99 %!assert (factor (1), 1)
100 %!test 100 %!test
101 %! for i = 2:20 101 %! for i = 2:20
102 %! p = factor (i); 102 %! pf = factor (i);
103 %! assert (prod (p), i); 103 %! assert (prod (pf), i);
104 %! assert (all (isprime (p))); 104 %! assert (all (isprime (pf)));
105 %! [p,n] = factor (i); 105 %! [pf, n] = factor (i);
106 %! assert (prod (p.^n), i); 106 %! assert (prod (pf.^n), i);
107 %! assert (all ([0,p] != [p,0])); 107 %! assert (all ([0,pf] != [pf,0]));
108 %! endfor 108 %! endfor
109 109
110 %% Test input validation 110 %% Test input validation
111 %!error factor () 111 %!error factor ()
112 %!error <Q must be a scalar integer> factor ([1,2]) 112 %!error factor (1,2)
113 %!error <Q must be a scalar integer> factor (1.5) 113 %!error <Q must be a real integer> factor (6i)
114 %!error <Q must be a real integer> factor ([1,2])
115 %!error <Q must be a real integer> factor (1.5)
114 116