Mercurial > hg > octave-nkf

comparison scripts/specfun/primes.m @ 8507:cadc73247d65

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
style fixes
author John W. Eaton <jwe@octave.org>
date Tue, 13 Jan 2009 14:08:36 -0500
parents bc982528de11
children eb63fbe60fab
comparison
equal deleted inserted replaced
8506:bc982528de11 8507:cadc73247d65
22 ## Return all primes up to @var{n}. 22 ## Return all primes up to @var{n}.
23 ## 23 ##
24 ## Note that if you need a specific number of primes, you can use the 24 ## Note that if you need a specific number of primes, you can use the
25 ## fact the distance from one prime to the next is on average 25 ## fact the distance from one prime to the next is on average
26 ## proportional to the logarithm of the prime. Integrating, you find 26 ## proportional to the logarithm of the prime. Integrating, you find
27 ## that there are about @math{k} primes less than @math{k \log ( 5 k )}. 27 ## that there are about @math{k} primes less than @math{k \log (5 k)}.
28 ## 28 ##
29 ## The algorithm used is called the Sieve of Erastothenes. 29 ## The algorithm used is called the Sieve of Erastothenes.
30 ## @end deftypefn 30 ## @end deftypefn
31 31
32 ## Author: Paul Kienzle 32 ## Author: Paul Kienzle
46 if (p > 100000) 46 if (p > 100000)
47 ## Optimization: 1/6 less memory, and much faster (asymptotically) 47 ## Optimization: 1/6 less memory, and much faster (asymptotically)
48 ## 100000 happens to be the cross-over point for Paul's machine; 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 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 50 ## of memory in Paul's machine, this saves .7 seconds out of 7 for
51 ## p=3e6. Hardly worthwhile, but Dirk reports better numbers. 51 ## p = 3e6. Hardly worthwhile, but Dirk reports better numbers.
52 lenm = floor ((p+1)/6); # length of the 6n-1 sieve 52 lenm = floor ((p+1)/6); # length of the 6n-1 sieve
53 lenp = 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 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 55 sievep = ones (1, lenp); # assume every number of form 6n+1 is prime
56 56