Eratosthenes 篩選求質數
December 2, 2021質數（prime number）是大於 1 的自然數（目前數學上公認 1 不是質數，因為 1 只有一個因數，也就是 1 本身），除了自身之外，無法被其他整數整除的數，至於非質數的自然數稱為合（成）數（Composite number）。
質數沒有公式可求，如何快速地求出質數則一直是程式設計人員與數學家努力的課題，在這邊介紹著名的 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 的倍數篩去：
-
23456789101112131415161718192021 …….. N - 接著要篩去 3 的倍數，此時留下的數中，小於 3 的就是質數：
-
23456789101112131415161718192021…….. N
再來篩去 5 的倍數，此時留下的數中，小於 5 的就是質數，依序地將 7、11 … 的倍數篩去，最後留下的數就是質數，這就是 Eratosthenes 篩選（Sieve of Eratosthenes）。
檢查的次數還可以再減少，事實上，只要檢查 6n + 1 與 6n + 5 就可以了，也就是 2 與 3 不用檢查，後續數字檢查直接跳過 2 與 3 的倍數。
自古至今對於質數，數學家有著許多的研究，其中值得留意的是，歐幾里得提出有些質數可以寫成 2ᵖ - 1 的形式，馬蘭·梅森 針對 p 為質數且 p <= 257，2ᵖ - 1 是個質數列出了梅森質數（Mersenne prime），為了紀念他，2ᵖ - 1 形式的數稱為梅森數（Mersenne number）。
除了求質數之外，還有另一個面向是給定一個數，判定它是否為質數，如果該數是個很大的數，會是個艱難的任務，因為不存在快速、有效率的分解方法，退而求其次地，存在一些非確定性的質數測試方式。
例如費馬小定理，若 p 為質數，與 p 互質的數 a，則 aᵖ⁻¹ ≡ 1 (mod p)。
a ≡ b (mod N) 代表同餘式，也就是以 N 為模數，a mod N 與 b mod N 得到的餘數恒等；費馬小定理說的就是，若 p 為質數，與 p 互質的數 a，aᵖ⁻¹ mod p 會餘 1，相對地，若有個數 n，與 n 互質的數 a，aⁿ⁻¹ mod n 不餘 1，n 不是質數。
使用費馬小定理時，要注意的是，若 p 為質數，與 p 互質的數 a，aᵖ⁻¹ ≡ 1 (mod p) 成立，但不表示有個數 n，與 n 互質的數 a，滿足 aⁿ⁻¹ ≡ 1 (mod n)，n 就是質數，例如 3215031751 並不是質數，然而可以找到 a 為 2、3、5、7 滿足 aⁿ⁻¹ ≡ 1 (mod n)。
Miller–Rabin primality test 基於費馬小定理加以改良，是不少程式庫提供的測試函式的實現原理（例如 Java 的 BigInteger 提供的 isProbablePrime 等方法），通過 Miller–Rabin 測試的數是質數的機率，比費馬小定理高許多（isProbablePrime 可以指定 certainty，傳回數是質數的機率是 1 - 1/2ᶜᵉʳᵗᵃⁱⁿᵗʸ），雖然通過測試並不保證是質數，不過這類測試的另一個價值是，未通過測試的一定不是質數。
長久以來，質數僅僅只是數學家的研究對象，直到 Ron Rivest、Adi Shamir、Leonard Adleman 提出了 RSA 加密演算出現，人們才意識到質數的重要性（以及商機），RSA 的出發點簡單來說，找出更大的質數很艱難，判定很大的數是否為質數也很難，如果兩個很大的質數相乘得到一個數，若只知道這個數，想反過來求出這個數是哪兩個質數相乘，更是一大數學難題。
雖然理論上可以透過一些演算加上暴力法破解，不過需要的運算資源與時間成本過高而不可行，說是這麼說，演算法中兩個質數的選擇還必須滿足一些條件（例如，最基本的不能選 2，畢竟相乘就會是偶數，一看就知道其 中一個質因數是 2），也沒有理論能證明，破解 RSA 的難度等價於質因數分解，也就是說，目前 RSA 的安全性建立在實際經驗，也就是目前尚未找到有效的演算可以破解 RSA，至於運算力更為的強大量子電腦呢？目前量子電腦尚未成熟，未來尚不可知，也許會演變成超大質數與量子電腦速度間的較勁？
程式實作
#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)))
