Eratosthenes 篩選求質數

December 2, 2021

質數是大於 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 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 …….. N

接著要篩去 3 的倍數,此時留下的數中,小於 3 的就是質數:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 …….. N

再來篩去 5 的倍數,此時留下的數中,小於 5 的就是質數,依序地將 7、11 … 的倍數篩去,最後留下的數就是質數,這就是 Eratosthenes 篩選(Sieve of Eratosthenes)。

檢查的次數還可以再減少,事實上,只要檢查 6n + 1 與 6n + 5 就可以了,也就是 2 與 3 不用檢查,後續數字檢查直接跳過 2 與 3 的倍數。

程式實作

#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);
        }
    }
}
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="")
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
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
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)))

分享到 LinkedIn 分享到 Facebook 分享到 Twitter