素数相关的算法
素数计算
素数相关的计算，主要有这几个方面：
- 列出某个范围内的所有素数；
- 判断某个数是否为素数；
- 其实是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"));
}
}
