octave_packages/m/specfun/isprime.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ## Copyright (C) 2010 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn {Function File} {} isprime (@var{x})
  22 ## Return a logical array which is true where the elements of @var{x} are
  23 ## prime numbers and false where they are not.
  24 ##
  25 ## If the maximum value in @var{x} is very large, then you should be using
  26 ## special purpose factorization code.
  27 ##
  28 ## @example
  29 ## @group
  30 ## isprime (1:6)
  31 ##     @result{} [0, 1, 1, 0, 1, 0]
  32 ## @end group
  33 ## @end example
  34 ## @seealso{primes, factor, gcd, lcm}
  35 ## @end deftypefn
  36
  37 function t = isprime (x)
  38
  39   if (nargin == 1)
  40     if (any ((x != floor (x) | x < 0)(:)))
  41       error ("isprime: needs positive integers");
  42     endif
  43     maxn = max (x(:));
  44     ## generate prime table of suitable length.
  45     maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized.
  46     pr = primes (maxp);
  47     ## quick search for table matches.
  48     t = lookup (pr, x, "b");
  49     ## take the rest.
  50     m = x(x > maxp);
  51     if (! isempty (m))
  52       ## there are still possible primes. filter them out by division.
  53       if (maxn <= intmax ("uint32"))
  54         m = uint32 (m);
  55       elseif (maxn <= intmax ("uint64"))
  56         m = uint64 (m);
  57       else
  58         warning ("isprime: too large integers being tested");
  59       endif
  60       pr = cast (pr(pr <= sqrt (maxn)), class (m));
  61       for p = pr
  62         m = m(rem (m, p) != 0);
  63         if (length (m) < length (pr) / 10)
  64           break;
  65         endif
  66       endfor
  67       pr = pr(pr > p);
  68       mm = arrayfun (@(x) all (rem (x, pr)), m);
  69       m = m(mm);
  70       if (! isempty (m))
  71         m = cast (sort (m), class (x));
  72         t |= lookup (m, x, "b");
  73       endif
  74     endif
  75
  76   else
  77     print_usage ();
  78   endif
  79
  80 endfunction
  81
  82
  83 %!assert (isprime (4), logical (0));
  84 %!assert (isprime (3), logical (1));
  85 %!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1]));
  86 %!error isprime ()
  87 %!error isprime (1, 2)