octave_packages/m/miscellaneous/list_primes.m

   1 ## Copyright (C) 1993-2012 John W. Eaton
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} list_primes ()
  21 ## @deftypefnx {Function File} {} list_primes (@var{n})
  22 ## List the first @var{n} primes.  If @var{n} is unspecified, the first
  23 ## 25 primes are listed.
  24 ##
  25 ## The algorithm used is from page 218 of the @TeX{}book.
  26 ## @seealso{primes, isprime}
  27 ## @end deftypefn
  28
  29 ## Author: jwe
  30
  31 function retval = list_primes (n)
  32
  33   if (nargin > 0)
  34     if (! isscalar (n))
  35       error ("list_primes: argument must be a scalar");
  36     endif
  37   endif
  38
  39   if (nargin == 0)
  40     n = 25;
  41   endif
  42
  43   if (n == 1)
  44     retval = 2;
  45     return;
  46   endif
  47
  48   if (n == 2)
  49     retval = [2; 3];
  50     return;
  51   endif
  52
  53   retval = zeros (1, n);
  54   retval (1) = 2;
  55   retval (2) = 3;
  56
  57   n = n - 2;
  58   i = 3;
  59   p = 5;
  60   while (n > 0)
  61
  62     is_prime = 1;
  63     is_unknown = 1;
  64     d = 3;
  65     while (is_unknown)
  66       a = fix (p / d);
  67       if (a <= d)
  68         is_unknown = 0;
  69       endif
  70       if (a * d == p)
  71         is_prime = 0;
  72         is_unknown = 0;
  73       endif
  74       d = d + 2;
  75     endwhile
  76
  77     if (is_prime)
  78       retval (i++) = p;
  79       n--;
  80     endif
  81     p = p + 2;
  82
  83   endwhile
  84
  85 endfunction
  86
  87 %!test
  88 %! assert (list_primes(), [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,\
  89 %!                        43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]);
  90 %! assert (list_primes(5), [2, 3, 5, 7, 11]);
  91