Mercurial > hg > octave-nkf

comparison scripts/specfun/isprime.m @ 19238:89e275a4f6f6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Allow isprime to handle complex numbers (bug #43041) There's no reason why we can't check for prime numbers in the Gaussian integers. Both the lcm and gcd functions already use the Gaussian integers, so primality testing can also happen in that domain. * isprime.m: Modify documentation and add new tests. (isprime): Call isgaussian prime for complex inputs. Relax the check on input arguments, and corresponding error message. (isgaussianprime): New function.
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sat, 23 Aug 2014 17:53:48 -0400
parents 920a400929ca
children 2f117c4b5cb0
comparison
equal deleted inserted replaced
19237:920a400929ca 19238:89e275a4f6f6
30 ## @group 30 ## @group
31 ## isprime (1:6) 31 ## isprime (1:6)
32 ## @result{} [0, 1, 1, 0, 1, 0] 32 ## @result{} [0, 1, 1, 0, 1, 0]
33 ## @end group 33 ## @end group
34 ## @end example 34 ## @end example
35 ##
36 ## If @code{class (@var{x})} is complex, then primality will be tested
37 ## in the domain of Gaussian integers. This means that some real prime
38 ## numbers will fail the primality test in this case. For example, 5 =
39 ## (1+2i)*(1-2i) and 2 = i*(1-i)^2.
40 ##
41 ## @example
42 ## @group
43 ## isprime ([i, 2, 3, 5])
44 ## @result{} [0, 0, 1, 0]
45 ## @end group
46 ## @end example
35 ## @seealso{primes, factor, gcd, lcm} 47 ## @seealso{primes, factor, gcd, lcm}
36 ## @end deftypefn 48 ## @end deftypefn
37 49
38 function t = isprime (x) 50 function t = isprime (x)
39 51
40 if (nargin != 1) 52 if (nargin != 1)
41 print_usage (); 53 print_usage ();
42 elseif (! isreal (x) || any ((x < 0 | x != fix (x))(:))) 54 elseif (any (fix (x) != x))
43 error ("isprime: X must be a non-negative integer"); 55 error ("isprime: X contains non-integer entries");
44 endif 56 endif
45 57
46 if (isempty (x)) 58 if (isempty (x))
47 t = x; 59 t = x;
48 return; 60 return;
49 endif 61 endif
62
63 if (iscomplex (x))
64 t = isgaussianprime (x);
65 return;
66 endif
67
68 ## Handle negative entries
69 x = abs (x);
50 70
51 maxn = max (x(:)); 71 maxn = max (x(:));
52 ## generate prime table of suitable length. 72 ## generate prime table of suitable length.
53 maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized. 73 maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized.
54 pr = primes (maxp); 74 pr = primes (maxp);
81 endif 101 endif
82 endif 102 endif
83 103
84 endfunction 104 endfunction
85 105
106 function t = isgaussianprime (z)
107 ## Assumed prime unless proven otherwise
108 t = true (size (z));
109
110 x = real (z);
111 y = imag (z);
112
113 ## If purely real or purely imaginary, ordinary prime test for
114 ## that complex part if that part is 3 mod 4.
115 xidx = y==0 & mod (x, 4) == 3;
116 yidx = x==0 & mod (y, 4) == 3;
117
118 t(xidx) &= isprime (x(xidx));
119 t(yidx) &= isprime (y(yidx));
120
121 ## Otherwise, prime if x^2 + y^2 is prime and NOT 3 mod 4.
122 zabs = x.^2 + y.^2;
123 zidx = mod (zabs, 4) != 3 & !xidx & !yidx;
124
125 t(zidx) &= isprime (zabs (zidx));
126 endfunction
127
128
86 129
87 %!assert (isprime (3), true) 130 %!assert (isprime (3), true)
88 %!assert (isprime (4), false) 131 %!assert (isprime (4), false)
132 %!assert (isprime (5i), false)
133 %!assert (isprime (7i), true)
134 %!assert (isprime ([1+2i, (2+3i)*(-1+2i)]), [true, false])
135 %!assert (isprime (-2), true)
136 %!assert (isprime (complex (-2)), false)
137 %!assert (isprime (2i), false)
138 %!assert (isprime ([i, 2, 3, 5]), [false, false, true, false])
139 %!assert (isprime (0), false)
89 %!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1])) 140 %!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1]))
90 141
91 %% Test input validation 142 %% Test input validation
92 %!error isprime () 143 %!error isprime ()
93 %!error isprime (1, 2) 144 %!error isprime (1, 2)
94 %!error <X must be a non-negative integer> isprime (i) 145 %!error <X contains non-integer entries> isprime (0.5i)
95 %!error <X must be a non-negative integer> isprime ([1 2; -3 4]) 146 %!error <X contains non-integer entries> isprime (0.5)
96 %!error <X must be a non-negative integer> isprime ([1 2; 3.1 4])
97