跳转到内容

算法实现/数学/素数生成

Add links
来自维基教科书，开放的书本，开放的世界

C

[编辑 | 编辑源代码]
#include <stdbool.h>
#include <stdio.h>

#define MAX 10000

/*
 * Prints all primes less than MAX using the Sieve of Eratosthenes.
 */
int
main(void)
{
	bool sieve[MAX];
	int i, j, primecount = 0, prime[MAX];

	for (i = 0; i < MAX; i++)
		sieve[i] = true;
	sieve[0] = sieve[1] = false;
	for (i = 2; i * i <= MAX; i++) {
		if (!sieve[i])
			continue;
		prime[primecount++] = i;
		for (j = i * i; j < MAX; j += i)
			sieve[j] = false;
	}
	for (i = 0; i < primecount; i++)
		printf("%d\n", prime[i]);
}

D

[编辑 | 编辑源代码]
int[] sieve(int n) {
    int[] sieve = new int[](n + 1); 
    int m = cast(int) sqrt(n);

    sieve[1..3] = 1;

    for (int i = 3; i <= m; i+=2) {
        for ( size_t j = 2; j < cast(int)sqrt(i); ++j) {
            if (sieve[i] == 1 || i % j != 0) {
                 continue;
            }
        }
        sieve[i] = 1;
    }
    
    return sieve;
}

Haskell

[编辑 | 编辑源代码]

对于每个数字，检查它是否可以被其平方根以下的任何素数整除。这是最优的试除法算法。

isPrime primes n = foldr (\p r -> p*p > n || (rem n p /= 0 && r))
                         True primes
primes = 2 : filter (isPrime primes) [3..]

注意，此代码利用 Haskell 的惰性求值，允许 primes 和 isPrime 相互定义。

埃拉托斯特尼筛法具有比试除法更好的理论复杂度。

primesTo m = 2 : sieve [3,5..m] where
  sieve s@(p:xs) | p*p > m = s
                 | True    = p : sieve (xs `minus` [p*p, p*p+2*p..m])

如果 (minus a b) 操作的复杂度为 ，正如它确实在命令式语言中的可变数组上一样，此代码将实现埃拉托斯特尼筛法的理论复杂度。不幸的是，它在链表中是 。

重新定义为无界定义，但由于树形折叠，速度更快，复杂度也更低。

primesU = 2 : _Y ((3 :) . minus [5,7..] . _U . map (\p-> [p*p, p*p+2*p..]))

_Y g = g (_Y g)                                      -- g . g . g . g . .... 
_U ((x:xs):t) = x : (union xs . _U . pairs) t        -- big_U, or unionAll 
   where  pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t

这用一个无界右深 union 树 代替了 primesTo 函数中 minus 操作的线性 左倾嵌套链。可以进一步通过轮分解实现相当大的常数倍速率提升。上面使用的工具函数（也在 data-ordlist 包的 Data.List.Ordered 模块中）是

minus (x:xs) (y:ys) = case compare x y of LT -> x : minus  xs (y:ys)
                                          EQ ->     minus  xs    ys
                                          GT ->     minus (x:xs) ys
minus  a      _     = a
union (x:xs) (y:ys) = case compare x y of LT -> x : union  xs (y:ys)
                                          EQ -> x : union  xs    ys
                                          GT -> y : union (x:xs) ys
union  a      []    = a
union  []     b     = b

使用数组很容易移除合数。处理自然地被分割成连续素数的平方之间的段。

import Data.List (inits, tails)
import Data.Array
 
primesSA = _Y (\ps -> 
                2 : [n | (r:q:_, px) <- (zip . tails . (2:) . map (^2)) ps (inits ps),
                         (n,True)    <- assocs (
                                          accumArray (\_ _ -> False) True (r+1,q-1)
                                    [(m,()) | p <- px, 
                                              let s=(r+p)`div`p*p, m <- [s,s+p..q-1]] )])

Java

[编辑 | 编辑源代码]

简单的朴素算法

import java.util.*;

// finds all prime numbers up to max
public static List<Integer> generatePrimes(int max)
{
    List<Integer> primes = new ArrayList<Integer>();

    // start from 2
    OUTERLOOP:
    for (int i = 2; i <= max; i++) {
        // try to divide it by all known primes
        for (Integer p : primes)
            if (i % p == 0)
                continue OUTERLOOP; // i is not prime

        // i is prime
        primes.add(i);
    }
    return primes;
}

埃拉托斯特尼筛法，来自这里

import java.util.*;

public class Sieve
{
    private BitSet sieve;

    private Sieve() {}

    private Sieve(int size) {
        sieve = new BitSet((size+1)/2);
    }

    private boolean is_composite(int k)
    {
        assert k >= 3 && (k % 2) == 1;
        return sieve.get((k-3)/2);
    }

    private void set_composite(int k)
    {
        assert k >= 3 && (k % 2) == 1;
        sieve.set((k-3)/2);
    }

    public static List<Integer> sieve_of_eratosthenes(int max)
    {
        Sieve sieve = new Sieve(max + 1); // +1 to include max itself
        for (int i = 3; i*i <= max; i += 2) {
            if (sieve.is_composite(i))
                continue;

            // We increment by 2*i to skip even multiples of i
            for (int multiple_i = i*i; multiple_i <= max; multiple_i += 2*i)
                sieve.set_composite(multiple_i);

        }
        List<Integer> primes = new ArrayList<Integer>();
        primes.add(2);
        for (int i = 3; i <= max; i += 2)
            if (!sieve.is_composite(i))
                primes.add(i);

        return primes;
    }
}

Javascript

[编辑 | 编辑源代码]
<!DOCTYPE html>
<html>
<body>

<button onclick="myFunction()">Sieve</button>

<script>
var sieve = [false, false];

function myFunction(){
  var n = sieve.length, nn = n * n;
  for(var i = n; i < nn; i++){ sieve.push(true); }

  for(var p = 0; p < n; p++){
    if(sieve[p]){
      for(var m = p * p; m < nn; m += p){
        sieve[m] = false;
      }
    }
  }

  var txt = document.body
   .appendChild(document.createElement('p'))
   .appendChild(document.createTextNode(''));
  for(var i = n; i < nn; i++){
    if(sieve[i]){ txt.nodeValue += ' ' + i; }
  }
}
</script>

</body>
</html>

Perl

[编辑 | 编辑源代码]
$MAX = 10000;

for($i=2; $i<$MAX; $i++){
	$sieve[$i]=1;
}
$sieve[0]=0;
$sieve[1]=0;
$primecount=0;
for($i=2; $i<$MAX; $i++){
	while($sieve[$i]==0 && $i<$MAX){$i++;}
	$prime[$primecount]=$i;
	for($j=2*$i; $j<$MAX; $j+=$i){$sieve[$j]=0;}
	$primecount++;
}

Phix

[编辑 | 编辑源代码]
constant limit = 1000
sequence primes = {}
sequence flags = repeat(1, limit)
for i=2 to floor(sqrt(limit)) do
    if flags[i] then
        for k=i*i to limit by i do
            flags[k] = 0
        end for
    end if
end for
for i=2 to limit do
    if flags[i] then
        primes &= i
    end if
end for
? primes

你可以从以下代码获得完全相同的输出（前 168 个素数）

?get_primes(-168)

Python

[编辑 | 编辑源代码]
def prime_sieve(n):
    """Generate the primes less than ``n`` using the Sieve of Eratosthenes."""
    a = [True] * n
    a[0] = a[1] = False
    for i, isprime in enumerate(a):
        if isprime:
            yield i
            for j in range(i * i, n, i):
                a[j] = False

def primes_wilson(n):
    """Generate the primes less than ``n`` using Wilson's theorem."""
    fac = 1
    for i in range(2, n):
        fac *= i - 1
        if (fac+1) % i == 0:
            yield i

Ruby

[编辑 | 编辑源代码]
def sieve_of_eratosthenes(max)
  arr=(2..max).to_a
  (2..Math::sqrt(max)).each do |i|
     arr.delete_if {|a|a % i == 0 && a!=i}
   end
   arr
end

# Example Call
puts sieve_of_eratosthenes(20)

注意，Ruby 有一个内置的埃拉托斯特尼筛法素数生成器。

require 'mathn'
for p in Prime.new
  puts p
end

F# / OCaml

[编辑 | 编辑源代码]
let prime_series_sieve (limit:bigint) = 
    let series = List.to_array [0I..limit]
    series.SetValue(0I, 1)

    let rec eliminate_multiples (n:bigint) (i:bigint) = 
        let index = (i * n)
        if index < bigint.Parse(series.Length.ToString()) then 
            series.SetValue(0I, (bigint.ToInt64(index)))
            eliminate_multiples n (i + 1I)
    
    for n in [2I..(limit/2I)] do
        eliminate_multiples n 2I
    
    series;;

C (位运算)

[编辑 | 编辑源代码]
#include <math.h>
#define MAX 1000000
const int S=(int)sqrt((double)MAX);

#define rP(n) (sieve[n>>6]|=(1<<((n>>1)&31)))
#define gP(n) sieve[n>>6]&(1<<((n>>1)&31))

unsigned sieve[(MAX>>6)+1]={0};
int i, j,k,l=0 , prime[MAX];
prime[0]=2;
for(i=3;i<=S;i+=2) if(!(gP(i))) {
	k=(i<<1);
	prime[l++]=i;
	for(j=i*i;j<=MAX;j+=k) rP(j);
}

for(;i<=MAX;i++) if(!(gP(i))) prime[l++]=i;

C (快速位运算)

[编辑 | 编辑源代码]

埃拉托斯特尼筛法的快速位运算实现。有点作弊，因为它只计算奇素数（即：你需要手动添加第一个素数 - "2"）。

在代码中还有 OpenMP #pragma，因此你可以尝试使用 -fopenmp 编译它。

#include <stdio.h>
#include <stdlib.h>
#include <string.h> 
 
#define MAX_N 1000000000
#define MAX MAX_N/2

/*
Good way to compile is:
gcc --std=c99 -O4 -o sieve_fast sieve_fast.c  
*/

int main(int argc, char * args[])
{
  size_t mem_sz = MAX / 8 + 1;
  size_t superpage_size = 2 * 1024 * 1024; // 2MB
  unsigned long cnt = 0;
  unsigned char * arr;

  if (posix_memalign((void **)&arr, superpage_size, mem_sz)) {
    perror("memalign");
    exit(1);
  }

  bzero(arr, mem_sz);

  #pragma omp parallel 
  for(register unsigned long n=1; ((n*(n+1))<<1)<MAX; n++)
  {
    register unsigned long t = n >> 3;
    if((arr[t] & (1 << (n & 7)))){
      //printf("Skipped: %lu\n", n*2+1);
      continue;
    }
    //printf("prime: %lu; n: %lu\n", n<<1+1, n);
    #pragma omp for schedule(static, 65536)
    for(register unsigned long c=(n*(n+1))<<1; c<MAX; c+=(n<<1)+1){
      t = c >> 3;
      arr[t] |= (1 << (c & 7));
      //printf("Sieving %lu with step %lu\n", c*2+1, n*2+1);
    }
  }
 
  #pragma omp parallel for reduction(+:cnt)
  for(register unsigned long n=1; n<MAX; n++)
  {
    register unsigned long t = n >> 3;
    if(!(arr[t] & (1 << (n & 7))))
    {
      //printf("%lu, ", n*2+1);
      cnt++;
    }  
  }
 
  printf("Found total: %lu\n", cnt);
}


取自 "https://wikibooks.cn/w/index.php?title=Algorithm_Implementation/Mathematics/Prime_number_generation&oldid=3796374"
华夏公益教科书