Mercurial > hg > octave-nkf
annotate scripts/specfun/primes.m @ 9141:c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Split section into "Exponents and Logarithms" and "Utility Functions"
Use Tex in many more of the doc strings for pretty printing in pdf format.
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Mon, 20 Apr 2009 17:16:09 -0700 |
| parents | eb63fbe60fab |
| children | 8c71a86c4bf4 |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2000, 2006, 2007, 2009 Paul Kienzle |
| 5827 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 5827 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5827 | 18 |
| 19 ## -*- texinfo -*- | |
| 20 ## @deftypefn {Function File} {} primes (@var{n}) | |
| 21 ## | |
| 22 ## Return all primes up to @var{n}. | |
| 23 ## | |
|
9141
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## The algorithm used is the Sieve of Erastothenes. |
|
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
25 ## |
| 5827 | 26 ## Note that if you need a specific number of primes, you can use the |
| 27 ## fact the distance from one prime to the next is on average | |
| 28 ## proportional to the logarithm of the prime. Integrating, you find | |
| 8507 | 29 ## that there are about @math{k} primes less than @math{k \log (5 k)}. |
| 5827 | 30 ## @end deftypefn |
| 31 | |
| 32 ## Author: Paul Kienzle | |
| 33 ## Author: Francesco Potort� | |
| 34 ## Author: Dirk Laurie | |
| 35 | |
| 36 function x = primes (p) | |
| 37 | |
| 38 if (nargin != 1) | |
| 39 print_usage (); | |
| 40 endif | |
| 41 | |
| 42 if (! isscalar (p)) | |
| 43 error ("primes: n must be a scalar"); | |
| 44 endif | |
| 45 | |
| 46 if (p > 100000) | |
| 8506 | 47 ## Optimization: 1/6 less memory, and much faster (asymptotically) |
| 5827 | 48 ## 100000 happens to be the cross-over point for Paul's machine; |
| 49 ## below this the more direct code below is faster. At the limit | |
| 50 ## of memory in Paul's machine, this saves .7 seconds out of 7 for | |
| 8507 | 51 ## p = 3e6. Hardly worthwhile, but Dirk reports better numbers. |
| 5827 | 52 lenm = floor ((p+1)/6); # length of the 6n-1 sieve |
| 53 lenp = floor ((p-1)/6); # length of the 6n+1 sieve | |
| 54 sievem = ones (1, lenm); # assume every number of form 6n-1 is prime | |
| 55 sievep = ones (1, lenp); # assume every number of form 6n+1 is prime | |
| 56 | |
| 8506 | 57 for i = 1:(sqrt(p)+1)/6 # check up to sqrt(p) |
| 5827 | 58 if (sievem(i)) # if i is prime, eliminate multiples of i |
| 59 sievem(7*i-1:6*i-1:lenm) = 0; | |
| 60 sievep(5*i-1:6*i-1:lenp) = 0; | |
| 61 endif # if i is prime, eliminate multiples of i | |
| 62 if (sievep(i)) | |
| 63 sievep(7*i+1:6*i+1:lenp) = 0; | |
| 64 sievem(5*i+1:6*i+1:lenm) = 0; | |
| 65 endif | |
| 66 endfor | |
| 67 x = sort([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]); | |
| 8506 | 68 elseif (p > 352) # nothing magical about 352; must be >2 |
| 5827 | 69 len = floor ((p-1)/2); # length of the sieve |
| 70 sieve = ones (1, len); # assume every odd number is prime | |
| 71 for i = 1:(sqrt(p)-1)/2 # check up to sqrt(p) | |
| 72 if (sieve(i)) # if i is prime, eliminate multiples of i | |
| 73 sieve(3*i+1:2*i+1:len) = 0; # do it | |
| 74 endif | |
| 75 endfor | |
| 76 x = [2, 1+2*find(sieve)]; # primes remaining after sieve | |
| 77 else | |
| 78 a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ... | |
| 79 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ... | |
| 80 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ... | |
| 81 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ... | |
| 82 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ... | |
| 83 293, 307, 311, 313, 317, 331, 337, 347, 349]; | |
| 84 x = a(a <= p); | |
| 85 endif | |
| 86 | |
| 87 endfunction |