Mercurial > octave-nkf
changeset 19080:89e275a4f6f6
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 | 1288a2f27769 |
| files | scripts/specfun/isprime.m |
| diffstat | 1 files changed, 55 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/isprime.m Sat Aug 23 17:38:56 2014 -0400 +++ b/scripts/specfun/isprime.m Sat Aug 23 17:53:48 2014 -0400 @@ -32,6 +32,18 @@ ## @result{} [0, 1, 1, 0, 1, 0] ## @end group ## @end example +## +## If @code{class (@var{x})} is complex, then primality will be tested +## in the domain of Gaussian integers. This means that some real prime +## numbers will fail the primality test in this case. For example, 5 = +## (1+2i)*(1-2i) and 2 = i*(1-i)^2. +## +## @example +## @group +## isprime ([i, 2, 3, 5]) +## @result{} [0, 0, 1, 0] +## @end group +## @end example ## @seealso{primes, factor, gcd, lcm} ## @end deftypefn @@ -39,8 +51,8 @@ if (nargin != 1) print_usage (); - elseif (! isreal (x) || any ((x < 0 | x != fix (x))(:))) - error ("isprime: X must be a non-negative integer"); + elseif (any (fix (x) != x)) + error ("isprime: X contains non-integer entries"); endif if (isempty (x)) @@ -48,6 +60,14 @@ return; endif + if (iscomplex (x)) + t = isgaussianprime (x); + return; + endif + + ## Handle negative entries + x = abs (x); + maxn = max (x(:)); ## generate prime table of suitable length. maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized. @@ -83,15 +103,44 @@ endfunction +function t = isgaussianprime (z) + ## Assumed prime unless proven otherwise + t = true (size (z)); + + x = real (z); + y = imag (z); + + ## If purely real or purely imaginary, ordinary prime test for + ## that complex part if that part is 3 mod 4. + xidx = y==0 & mod (x, 4) == 3; + yidx = x==0 & mod (y, 4) == 3; + + t(xidx) &= isprime (x(xidx)); + t(yidx) &= isprime (y(yidx)); + + ## Otherwise, prime if x^2 + y^2 is prime and NOT 3 mod 4. + zabs = x.^2 + y.^2; + zidx = mod (zabs, 4) != 3 & !xidx & !yidx; + + t(zidx) &= isprime (zabs (zidx)); +endfunction + + %!assert (isprime (3), true) %!assert (isprime (4), false) +%!assert (isprime (5i), false) +%!assert (isprime (7i), true) +%!assert (isprime ([1+2i, (2+3i)*(-1+2i)]), [true, false]) +%!assert (isprime (-2), true) +%!assert (isprime (complex (-2)), false) +%!assert (isprime (2i), false) +%!assert (isprime ([i, 2, 3, 5]), [false, false, true, false]) +%!assert (isprime (0), false) %!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1])) %% Test input validation %!error isprime () %!error isprime (1, 2) -%!error <X must be a non-negative integer> isprime (i) -%!error <X must be a non-negative integer> isprime ([1 2; -3 4]) -%!error <X must be a non-negative integer> isprime ([1 2; 3.1 4]) - +%!error <X contains non-integer entries> isprime (0.5i) +%!error <X contains non-integer entries> isprime (0.5)