Mercurial > octave-nkf

comparison scripts/specfun/primes.m @ 5827:1fe78adb91bc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-22 06:25:14 by jwe]
author jwe
date Mon, 22 May 2006 06:25:14 +0000
parents
children 93c65f2a5668
comparison
equal deleted inserted replaced
5826:6c6ff9b82577 5827:1fe78adb91bc
1 ## Copyright (C) 2000 Paul Kienzle
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 2, or (at your option)
8 ## 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, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} primes (@var{n})
22 ##
23 ## Return all primes up to @var{n}.
24 ##
25 ## Note that if you need a specific number of primes, you can use the
26 ## fact the distance from one prime to the next is on average
27 ## proportional to the logarithm of the prime. Integrating, you find
28 ## that there are about @math{k} primes less than @math{k \log ( 5 k )}.
29 ##
30 ## The algorithm used is called the Sieve of Erastothenes.
31 ## @end deftypefn
32
33 ## Author: Paul Kienzle
34 ## Author: Francesco Potort́
35 ## Author: Dirk Laurie
36
37 function x = primes (p)
38
39 if (nargin != 1)
40 print_usage ();
41 endif
42
43 if (! isscalar (p))
44 error ("primes: n must be a scalar");
45 endif
46
47 if (p > 100000)
48 ## optimization: 1/6 less memory, and much faster (asymptotically)
49 ## 100000 happens to be the cross-over point for Paul's machine;
50 ## below this the more direct code below is faster. At the limit
51 ## of memory in Paul's machine, this saves .7 seconds out of 7 for
52 ## p=3e6. Hardly worthwhile, but Dirk reports better numbers.
53 lenm = floor ((p+1)/6); # length of the 6n-1 sieve
54 lenp = floor ((p-1)/6); # length of the 6n+1 sieve
55 sievem = ones (1, lenm); # assume every number of form 6n-1 is prime
56 sievep = ones (1, lenp); # assume every number of form 6n+1 is prime
57
58 for i = 1:(sqrt(p)+1)/6 # check up to sqrt(p)
59 if (sievem(i)) # if i is prime, eliminate multiples of i
60 sievem(7*i-1:6*i-1:lenm) = 0;
61 sievep(5*i-1:6*i-1:lenp) = 0;
62 endif # if i is prime, eliminate multiples of i
63 if (sievep(i))
64 sievep(7*i+1:6*i+1:lenp) = 0;
65 sievem(5*i+1:6*i+1:lenm) = 0;
66 endif
67 endfor
68 x = sort([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]);
69 elseif (p > 352) # nothing magical about 352; just has to be greater than 2
70 len = floor ((p-1)/2); # length of the sieve
71 sieve = ones (1, len); # assume every odd number is prime
72 for i = 1:(sqrt(p)-1)/2 # check up to sqrt(p)
73 if (sieve(i)) # if i is prime, eliminate multiples of i
74 sieve(3*i+1:2*i+1:len) = 0; # do it
75 endif
76 endfor
77 x = [2, 1+2*find(sieve)]; # primes remaining after sieve
78 else
79 a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ...
80 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ...
81 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ...
82 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ...
83 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ...
84 293, 307, 311, 313, 317, 331, 337, 347, 349];
85 x = a(a <= p);
86 endif
87
88 endfunction