素數計算
素數相關的計算，主要有這幾個方面：
- 列出某個范圍內的所有素數；
- 判斷某個數是否為素數；
- 其實是2)的擴展，快速獲取一個大素數
列出某個范圍的所有素數
這個可以分成兩種類型，一種是列出從1至N的所有素數，另一個是在一個較大數值的區間，列出所有素數。
列出1至N的所有素數
1) 普通計算方式, 校驗每個數字
優化的幾處:
- 判斷是否整除時, 除數使用小於自身的平方根的素數
- 大於3的素數, 都在6的整數倍兩側, 即 6m - 1 和 6m + 1
public class DemoPrime {
private int[] primes;
private int max;
private int pos;
private int total;
public DemoPrime(int max) {
this.max = max;
int length = max / 3;
primes = new int[length];
pos = 0;
total = 0;
}
private void put(int prime) {
primes[pos] = prime;
if (pos < primes.length - 1) {
pos++;
} else {
throw new RuntimeException("Length exceed");
}
}
private boolean isPrime(int num) {
int limit = (int)Math.sqrt(num);
for (int i = 0; i < pos; i++) {
if (primes[i] > limit) {
break;
}
total++;
if (num % primes[i] == 0) return false;
}
return true;
}
public void calculate() {
put(2);
put(3);
int val = 1;
for (int i = 0; val <= max - 6;) {
val += 4;
if (isPrime(val)) {
put(val);
}
val += 2;
if (isPrime(val)) {
put(val);
}
}
System.out.println("Tried: " + total);
}
public void print() {
System.out.println("Total: " + pos);
/*for (int i = 0; i < pos; i++) {
System.out.println(primes[i]);
}*/
}
public static void main(String[] args) {
DemoPrime dp = new DemoPrime(10000000);
dp.calculate();
dp.print();
}
}
2) 使用數組填充的方式，即Sieve of Eratosthenes 埃拉托色尼篩法
一次性創建大小為N的int數組的方式, 在每得到一個素數時, 將其整數倍的下標(除自身以外)的元素都置位, 並且只需要遍歷到N的平方根處. 最后未置位的元素即為素數, 在過程中可以統計素數個數. 這種方法比前一種效率高一個數量級.
public class DemoPrime2 {
private int[] cells;
private int max;
private int total;
public DemoPrime2(int max) {
this.max = max;
cells = new int[max];
total = max;
}
private void put(int prime) {
int i = prime + prime;
while (i < max) {
if (cells[i] == 0) {
cells[i] = 1;
total--;
}
i += prime;
}
}
public void calculate() {
total -= 2; // Exclude 0 and 1
put(2);
put(3);
int limit = (int)Math.sqrt(max);
for (int i = 4; i <= limit; i++) {
if (cells[i] == 0) {
put(i);
}
}
}
public void print() {
System.out.println("Total: " + total);
/*for (int i = 2; i < max; i++) {
if (cells[i] == 0) {
System.out.println(i);
}
}*/
}
public static void main(String[] args) throws InterruptedException {
DemoPrime2 dp = new DemoPrime2(10000000);
Thread.sleep(1000L);
long ts = System.currentTimeMillis();
dp.calculate();
dp.print();
long elapse = System.currentTimeMillis() - ts;
System.out.println("Time: " + elapse);
}
}
在一個較大數值的區間，列出所有素數
這個問題等價於，在這個區間里對每一個數判斷是否為素數
判斷大數是否為素數
對大數進行素數判斷，常用的是Miller Rabin算法
https://en.wikibooks.org/wiki/Algorithm_Implementation/Mathematics/Primality_Testing
在JDK中，BigInteger有一個isProbablePrime(int certainty)方法用於判斷大數是否為素數，里面聯合使用了Miller-Rabin和Lucas-Lehmer算法。后者盧卡斯萊默算法僅用於檢測值為2p - 1的數的素性。
Miller-Rabin算法
對於大數的素性判斷，目前Miller-Rabin算法應用最廣泛。Miller Rabin算法基於費馬小定理和二次探測定理，其中
費馬小定理：若P為素數，且有0<a<p，則a^(p-1) % p = 1
二次探測定理：x*x % p == 1, 若P為素數, 則x的解只能是x = 1或者x = p - 1
一般底數為隨機選取，但當待測數不太大時，選擇測試底數就有一些技巧了。比如，如果被測數小於4 759 123 141，那么只需要測試三個底數2, 7和61就足夠了。當然，你測試的越多，正確的范圍肯定也越大。如果你每次都用前7個素數(2, 3, 5, 7, 11, 13和17)進行測試，所有不超過341 550 071 728 320的數都是正確的。如果選用2, 3, 7, 61和24251作為底數，那么10^16內唯一的強偽素數為46 856 248 255 981。這樣的一些結論使得Miller-Rabin算法在OI中非常實用。通常認為，Miller-Rabin素性測試的正確率可以令人接受，隨機選取k個底數進行測試算法的失誤率大概為4^(-k)。
/**
* * Java Program to Implement Miller Rabin Primality Test Algorithm
**/
import java.util.Scanner;
import java.util.Random;
import java.math.BigInteger;
/** Class MillerRabin **/
public class MillerRabin {
/** Function to check if prime or not **/
public boolean isPrime(long n, int iteration) {
/** base case **/
if (n == 0 || n == 1)
return false;
/** base case - 2 is prime **/
if (n == 2)
return true;
/** an even number other than 2 is composite **/
if (n % 2 == 0)
return false;
long s = n - 1;
while (s % 2 == 0)
s /= 2;
Random rand = new Random();
for (int i = 0; i < iteration; i++) {
long r = Math.abs(rand.nextLong());
long a = r % (n - 1) + 1, temp = s;
long mod = modPow(a, temp, n);
while (temp != n - 1 && mod != 1 && mod != n - 1) {
mod = mulMod(mod, mod, n);
temp *= 2;
}
if (mod != n - 1 && temp % 2 == 0)
return false;
}
return true;
}
/** Function to calculate (a ^ b) % c **/
public long modPow(long a, long b, long c) {
long res = 1;
for (int i = 0; i < b; i++) {
res *= a;
res %= c;
}
return res % c;
}
/** Function to calculate (a * b) % c **/
public long mulMod(long a, long b, long mod) {
return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
}
/** Main function **/
public static void main(String[] args) {
Scanner scan = new Scanner(System.in);
System.out.println("Miller Rabin Primality Algorithm Test\n");
/** Make an object of MillerRabin class **/
MillerRabin mr = new MillerRabin();
/** Accept number **/
System.out.println("Enter number\n");
long num = scan.nextLong();
/** Accept number of iterations **/
System.out.println("\nEnter number of iterations");
int k = scan.nextInt();
/** check if prime **/
boolean prime = mr.isPrime(num, k);
if (prime)
System.out.println("\n" + num + " is prime");
else
System.out.println("\n" + num + " is composite");
}
}
Miller-Rabin算法是一個RP算法。RP是時間復雜度的一種，主要針對判定性問題。一個算法是RP算法表明它可以在多項式的時間里完成，對於答案為否定的情形能夠准確做出判斷，但同時它也有可能把對的判成錯的（錯誤概率不能超過1/2）。RP算法是基於隨機化的，因此多次運行該算法可以降低錯誤率。還有其它的素性測試算法也是概率型的，比如Solovay-Strassen算法。
/**
* Class SolovayStrassen
**/
public class SolovayStrassen {
/**
* Function to calculate jacobi (a/b)
**/
public long Jacobi(long a, long b) {
if (b <= 0 || b % 2 == 0)
return 0;
long j = 1L;
if (a < 0) {
a = -a;
if (b % 4 == 3)
j = -j;
}
while (a != 0) {
while (a % 2 == 0) {
a /= 2;
if (b % 8 == 3 || b % 8 == 5)
j = -j;
}
long temp = a;
a = b;
b = temp;
if (a % 4 == 3 && b % 4 == 3)
j = -j;
a %= b;
}
if (b == 1)
return j;
return 0;
}
/**
* Function to check if prime or not
**/
public boolean isPrime(long n, int iteration) {
/** base case **/
if (n == 0 || n == 1)
return false;
/** base case - 2 is prime **/
if (n == 2)
return true;
/** an even number other than 2 is composite **/
if (n % 2 == 0)
return false;
Random rand = new Random();
for (int i = 0; i < iteration; i++) {
long r = Math.abs(rand.nextLong());
long a = r % (n - 1) + 1;
long jacobian = (n + Jacobi(a, n)) % n;
long mod = modPow(a, (n - 1) / 2, n);
if (jacobian == 0 || mod != jacobian)
return false;
}
return true;
}
/**
* Function to calculate (a ^ b) % c
**/
public long modPow(long a, long b, long c) {
long res = 1;
for (int i = 0; i < b; i++) {
res *= a;
res %= c;
}
return res % c;
}
/**
* Main function
**/
public static void main(String[] args) {
Scanner scan = new Scanner(System.in);
System.out.println("SolovayStrassen Primality Algorithm Test\n");
/** Make an object of SolovayStrassen class **/
SolovayStrassen ss = new SolovayStrassen();
/** Accept number **/
System.out.println("Enter number\n");
long num = scan.nextLong();
/** Accept number of iterations **/
System.out.println("\nEnter number of iterations");
int k = scan.nextInt();
/** check if prime **/
boolean prime = ss.isPrime(num, k);
if (prime)
System.out.println("\n" + num + " is prime");
else
System.out.println("\n" + num + " is composite");
}
}
AKS算法
AKS最關鍵的重要性在於它是第一個被發表的一般的、多項式的、確定性的和無仰賴的素數判定算法。先前的算法至多達到了其中三點，但從未達到全部四個。
- AKS算法可以被用於檢測任何一般的給定數字是否為素數。很多已知的高速判定算法只適用於滿足特定條件的素數。例如，盧卡斯-萊默檢驗法僅對梅森素數適用，而Pépin測試僅對費馬數適用。
- 算法的最長運行時間可以被表為一個目標數字長度的多項式。ECPP和APR能夠判斷一個給定數字是否為素數，但無法對所有輸入給出多項式時間范圍。
- 算法可以確定性地判斷一個給定數字是否為素數。隨機測試算法，例如米勒-拉賓檢驗和Baillie–PSW，可以在多項式時間內對給定數字進行校驗，但只能給出概率性的結果。
- AKS算法並未“仰賴”任何未證明猜想。一個反例是確定性米勒檢驗：該算法可以在多項式時間內對所有輸入給出確定性結果，但其正確性卻基於尚未被證明的廣義黎曼猜想。
AKS算法的時間復雜度是 O(log(n)), 比Miller-Rabin要慢
/***************************************************************************
* Team
**************
* Arijit Banerjee
* Suchit Maindola
* Srikanth Manikarnike
*
**************
* This is am implementation of Agrawal–Kayal–Saxena primality test in java.
*
**************
* The algorithm is -
* 1. l <- log n
* 2. for i<-2 to l
* a. if an is a power fo l
* return COMPOSITE
* 3. r <- 2
* 4. while r < n
* a. if gcd( r, n) != 1
* return COMPSITE
* b. if sieve marked n as PRIME
* q <- largest factor (r-1)
* o < - r-1 / q
* k <- 4*sqrt(r) * l
* if q > k and n <= r
* return PRIME
* c. x = 2
* d. for a <- 1 to k
* if (x + a) ^n != x^n + mod (x^r - 1, n)
* return COMPOSITE
* e. return PRIME
*/
public class DemoAKS {
private int log;
private boolean sieveArray[];
private int SIEVE_ERATOS_SIZE = 100000000;
/* aks constructor */
public DemoAKS(BigInteger input) {
sieveEratos();
boolean result = checkIsPrime(input);
if (result) {
System.out.println("1");
} else {
System.out.println("0");
}
}
/* function to check if a given number is prime or not */
public boolean checkIsPrime(BigInteger n) {
BigInteger lowR, powOf, x, leftH, rightH, fm, aBigNum;
int totR, quot, tm, aCounter, aLimit, divisor;
log = (int) logBigNum(n);
if (findPower(n, log)) {
return false;
}
lowR = new BigInteger("2");
x = lowR;
totR = lowR.intValue();
for (lowR = new BigInteger("2");
lowR.compareTo(n) < 0;
lowR = lowR.add(BigInteger.ONE)) {
if ((lowR.gcd(n)).compareTo(BigInteger.ONE) != 0) {
return false;
}
totR = lowR.intValue();
if (checkIsSievePrime(totR)) {
quot = largestFactor(totR - 1);
divisor = (int) (totR - 1) / quot;
tm = (int) (4 * (Math.sqrt(totR)) * log);
powOf = mPower(n, new BigInteger("" + divisor), lowR);
if (quot >= tm && (powOf.compareTo(BigInteger.ONE)) != 0) {
break;
}
}
}
fm = (mPower(x, lowR, n)).subtract(BigInteger.ONE);
aLimit = (int) (2 * Math.sqrt(totR) * log);
for (aCounter = 1; aCounter < aLimit; aCounter++) {
aBigNum = new BigInteger("" + aCounter);
leftH = (mPower(x.subtract(aBigNum), n, n)).mod(n);
rightH = (mPower(x, n, n).subtract(aBigNum)).mod(n);
if (leftH.compareTo(rightH) != 0) return false;
}
return true;
}
/* function that computes the log of a big number*/
public double logBigNum(BigInteger bNum) {
String str;
int len;
double num1, num2;
str = "." + bNum.toString();
len = str.length() - 1;
num1 = Double.parseDouble(str);
num2 = Math.log10(num1) + len;
return num2;
}
/*function that computes the log of a big number input in string format*/
public double logBigNum(String str) {
String s;
int len;
double num1, num2;
len = str.length();
s = "." + str;
num1 = Double.parseDouble(s);
num2 = Math.log10(num1) + len;
return num2;
}
/* function to compute the largest factor of a number */
public int largestFactor(int num) {
int i;
i = num;
if (i == 1) return i;
while (i > 1) {
while (sieveArray[i] == true) {
i--;
}
if (num % i == 0) {
return i;
}
i--;
}
return num;
}
/*function given a and b, computes if a is power of b */
public boolean findPowerOf(BigInteger bNum, int val) {
int l;
double len;
BigInteger low, high, mid, res;
low = new BigInteger("10");
high = new BigInteger("10");
len = (bNum.toString().length()) / val;
l = (int) Math.ceil(len);
low = low.pow(l - 1);
high = high.pow(l).subtract(BigInteger.ONE);
while (low.compareTo(high) <= 0) {
mid = low.add(high);
mid = mid.divide(new BigInteger("2"));
res = mid.pow(val);
if (res.compareTo(bNum) < 0) {
low = mid.add(BigInteger.ONE);
} else if (res.compareTo(bNum) > 0) {
high = mid.subtract(BigInteger.ONE);
} else if (res.compareTo(bNum) == 0) {
return true;
}
}
return false;
}
/* creates a sieve array that maintains a table for COMPOSITE-ness
* or possibly PRIME state for all values less than SIEVE_ERATOS_SIZE
*/
public boolean checkIsSievePrime(int val) {
if (sieveArray[val] == false) {
return true;
} else {
return false;
}
}
long mPower(long x, long y, long n) {
long m, p, z;
m = y;
p = 1;
z = x;
while (m > 0) {
while (m % 2 == 0) {
m = (long) m / 2;
z = (z * z) % n;
}
m = m - 1;
p = (p * z) % n;
}
return p;
}
/* function, given a and b computes if a is a power of b */
boolean findPower(BigInteger n, int l) {
int i;
for (i = 2; i < l; i++) {
if (findPowerOf(n, i)) {
return true;
}
}
return false;
}
BigInteger mPower(BigInteger x, BigInteger y, BigInteger n) {
BigInteger m, p, z, two;
m = y;
p = BigInteger.ONE;
z = x;
two = new BigInteger("2");
while (m.compareTo(BigInteger.ZERO) > 0) {
while (((m.mod(two)).compareTo(BigInteger.ZERO)) == 0) {
m = m.divide(two);
z = (z.multiply(z)).mod(n);
}
m = m.subtract(BigInteger.ONE);
p = (p.multiply(z)).mod(n);
}
return p;
}
/* array to populate sieve array
* the sieve array looks like this
*
* y index -> 0 1 2 3 4 5 6 ... n
* x index 1
* | 2 T - T - T ...
* \/ 3 T - - T ...
* 4 T - - ...
* . T - ...
* . T ...
* n
*
*
*
*
*/
public void sieveEratos() {
int i, j;
sieveArray = new boolean[SIEVE_ERATOS_SIZE + 1];
sieveArray[1] = true;
for (i = 2; i * i <= SIEVE_ERATOS_SIZE; i++) {
if (!sieveArray[i]) {
for (j = i * i; j <= SIEVE_ERATOS_SIZE; j += i) {
sieveArray[j] = true;
}
}
}
}
public static void main(String[] args) {
new DemoAKS(new BigInteger("100000217"));
}
}