Eratosthenes篩選求質數


說明

大於1的數,除了自身之外,無法被其它整數整除的數稱之為質數,要求質數很簡單,但如何快速的求出質數則一直是程式設計人員與數學家努力的課題,在這邊介紹一個著名的 Eratosthenes求質數方法。

解法

首先知道這個問題可以使用迴圈來求解,將一個指定的數除以所有小於它的數,若可以整除就不是質數,然而如何減少迴圈的檢查次數?如何求出小於N的所有質數?

首先假設要檢查的數是N好了,則事實上只要檢查至N的開根號就可以了,道理很簡單,假設A*B = N,如果A大於N的開根號,則事實上在小於A之前的檢查就可以先檢查到B這個數可以整除N。不過在程式中使用開根號會精確度的問題,所以可以使用 i * i <= N進行檢查,且執行更快。

再來假設有一個篩子存放1∼N,例如:
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ........ N

先將2的倍數篩去:
2 3 5 7 9 11 13 15 17 19 21 ........ N

再將3的倍數篩去:
2 3 5 7 11 13 17 19 ........ N

再來將5的倍數篩去,再來將7的質數篩去,再來將11的倍數篩去........,如此進行到最後留下的數就都是質數,這就是Eratosthenes篩選方法(Eratosthenes Sieve Method)。

檢查的次數還可以再減少,事實上,只要檢查6n + 1與6n + 5就可以了,也就是直接跳過2與3的倍數,使得程式中的if的檢查動作可以減少。

實作:Toy    C    Java    Python    Scala    Ruby    JavaScript    Haskell

from '/lib/math' import pow
N = 1000

def sieve(primes, i) {
    if primes.get(i) {
        iterate(2, j -> j * i <= N).forEach(j -> primes.set(j * i, 0)) 
    }
}

def sixPlus1Or5(primes) {
    r = range(1, i -> 6 * i + 5 <= N)
    r.forEach(i -> primes.set(6 * i + 1, 1))
    r.forEach(i -> primes.set(6 * i + 5, 1))
    if 6 * r.length() + 1 <= N { 
        primes.set(6 * r.length() + 1, 1)
    }
}

def sieveSixPlus1Or5(primes) {
    r = range(0, i -> pow(6 * i + 5, 2) <= N)
    r.forEach(i -> sieve(primes, 6 * i + 1))
    r.forEach(i -> sieve(primes, 6 * i + 5))
    n = r.length()
    if pow(6 * n + 1, 2) <= N { 
        sieve(primes, 6 * n + 1)
    }  
}

def create() {
    primes = List.create(N + 1, 0)
    primes.set(2, 1)
    primes.set(3, 1) 
    primes.set(5, 1)

    sixPlus1Or5(primes)
    sieveSixPlus1Or5(primes)

    return primes
}

primes = create()
println(iterate(2, N).select(i -> primes.get(i)))

#include <stdio.h> 
#include <stdlib.h>

#define N 1000

void create(int*);
void filter(int*, int);

int main(void) {
int primes[N + 1] = {0};
create(primes);

int i;
for(i = 2; i < N; i++) if(primes[i]) {
printf("%4d", i);
}
 
return 0;
}

void create(int* primes) {
primes[2] = primes[3] = primes[5] = 1;

int i;
for(i = 1;6 * i + 5 <= N; i++) {
primes[6 * i + 1] = primes[6 * i + 5] = 1;
}
if(6 * i + 1 <= N) { primes[6 * i + 1] = 1; }

int n;
for(n = 0;(6 * n + 5) * (6 * n + 5) <= N; n++) {
filter(primes, 6 * n + 1);
filter(primes, 6 * n + 5);
}
if((6 * n + 1) * (6 * n + 1) <= N) { filter(primes, 6 * n + 1); }
}

void filter(int* primes, int i) {
if(primes[i]) {
int j;
for(j = 2; j * i <= N; j++) {
primes[j * i] = 0;
}
}
}

import java.util.*;

public class Prime {
public static List<Integer> create(int max) {
LinkedList<Integer> primes = new LinkedList<>();
primes.add(2); primes.add(3); primes.add(5);
for(int i = 1; i <= (max - 5) / 6; i++) {
primes.add(6 * i + 1); primes.add(6 * i + 5);
}
if(primes.getLast() + 2 <= max) {
primes.add(primes.getLast() + 2);
}

for(int i = 2; i < primes.size(); i++) {
Integer p = primes.get(i);
for(int j = 2; j * p <= primes.getLast(); j++) {
primes.remove(Integer.valueOf(j * p));
}
}
return primes;
}

public static void main(String[] args) {
for(Integer p : create(1000)) {
System.out.printf("%4d", p);
}
}
}

  • Python
def flatten(list):
return [item for subList in list for item in subList]

def iterate(counter, primes):
if counter < len(primes):
newps = primes[0:counter + 1] + list(
filter(lambda p: p % primes[counter] != 0, primes[counter + 1:]))
return iterate(counter + 1, newps)
else:
return primes

def create(max):
comp = flatten(
[[6 * i + 1, 6 * i + 5] for i in range(1, (max - 5) // 6 + 1)])
primes = [2, 3, 5] + comp + \
([comp[-1] + 2] if comp[-1] + 2 <= max else [])
return iterate(2, primes)

for p in create(1000):
print("%4d" % p, end="")

  • Scala
def create(max: Int) = {
val comp = (for(i <- 1 to (max - 5) / 6)
yield List(6 * i + 1, 6 * i + 5)).toList.flatten
val primes = List(2, 3, 5) ++ comp ++
(if(comp.last + 2 <= max) List(comp.last + 2) else Nil)
iterate(2, primes)
}

def iterate(counter: Int, primes: List[Int]): List[Int] = {
if(counter < primes.size) {
val newps = primes.take(counter + 1) ++
primes.drop(counter + 1).filter(p => p % primes(counter) != 0)
iterate(counter + 1, newps)
} else primes
}

create(1000).foreach(printf("%4d", _))

def iterate(counter, primes)
if counter < primes.size
newps = primes.take(counter + 1) + primes.drop(counter + 1)
.select { |p| p % primes[counter] != 0}
iterate(counter + 1, newps)
else
primes
end
end

def create(max)
comp = (1..(max - 5) / 6).map {|i| [6 * i + 1, 6 * i + 5] }.flatten
primes = [2, 3, 5] + comp + (comp[-1] + 2 <= max ? [comp[-1] + 2] : [])
iterate(2, primes)
end

create(1000).each do |p|
printf("%4d", p)
end

  • JavaScript
Array.prototype.reduce = function(init, f) {
var value = init;
for(var i = 0; i < this.length; i++) {
value = f(value, this[i]);
}
return value;
};

function range(from, to) {
var r = [];
for(var c = 0, i = from; i < to; c++, i++) { r[i] = i; }
return r;
}

function iterate(counter, primes) {
if(counter < primes.length) {
var newps = primes.slice(0, counter + 1).concat(
primes.slice(counter + 1, primes.length)
.filter(function(p) {
return p % primes[counter] != 0;
}));
return iterate(counter + 1, newps);
} else {
return primes;
}
}

function create(max) {
var comp = range(1, parseInt((max - 5) / 6) + 1).map(function(i) {
return [6 * i + 1, 6 * i + 5];
}).reduce([], function(ac, list) {
return ac.concat(list);
});
var last = comp.pop();
var primes = [2, 3, 5].concat(
comp, last + 2 <= max ? [last, last + 2] : [last]);
return iterate(2, primes);
}

print(create(1000).join(' '));

flatten list = [item | subList <- list, item <- subList]

iter counter primes =
if counter < length primes then
let newps = take (counter + 1) primes ++
(filter (\p -> p `mod` (primes !! counter) /= 0) \$
drop (counter + 1) primes)
in iter (counter + 1) newps
else primes

create max = iter 2 primes
where comp = flatten [[6 * i + 1, 6 * i + 5] | i <- [1..(max - 5) `div` 6]]
primes = [2, 3, 5] ++ comp ++
if last comp + 2 <= max then [last comp + 2] else []

main = print \$ create 1000