素数相关的算法

素数计算

素数相关的计算，主要有这几个方面：

1. 列出某个范围内的所有素数；
2. 判断某个数是否为素数；
3. 其实是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 1put(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算法。后者卢卡斯莱默算法仅用于检测值为2^p - 1的数的素性。

Miller-Rabin算法

对于大数的素性判断，目前Miller-Rabin算法应用最广泛。Miller Rabin算法基于费马小定理和二次探测定理，其中

费马小定理：若P为素数，且有0

二次探测定理：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");elseSystem.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");elseSystem.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"));}
}

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！