add Miller–Rabin_primality_tes

Author: xtaci

include/imath.h

@@ -14,6 +14,7 @@
 
 #include <stdint.h>
 #include <string.h>
+#include <limits.h>
 
 namespace alg
 {
@@ -67,6 +68,27 @@ namespace alg
 }
 return result;
 }
+
+/**
+ * Count the consecutive zero bits (trailing) on the right linearly
+ */
+static int ZerosR(unsigned int v) {
+int c; // output: c will count v's trailing zero bits,
+// so if v is 1101000 (base 2), then c will be 3
+if (v)
+{
+ v = (v ^ (v - 1)) >> 1; // Set v's trailing 0s to 1s and zero rest
+ for (c = 0; v; c++) {
+v >>= 1;
+ }
+}
+else
+{
+ c = CHAR_BIT * sizeof(v);
+}
+
+return c;
+}
 }
 
 #endif //

include/prime.h

@@ -8,25 +8,81 @@
  * PRIME TEST FUNCTION
  *
  * http://en.wikipedia.org/wiki/Primality_test
+ * http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
  *
  ******************************************************************************/
 #ifndef __PRIME_H__
 #define __PRIME_H__
 
+#include <stdlib.h>
+#include <math.h>
+#include "imath.h"
+
 namespace alg
 {
 /**
  * check whether a given number is a prime number.
+ * using naive method.
  */
-static inline int is_prime(unsigned int n)
+static bool test_prime(unsigned int n)
 {
-unsigned int p;
-if (!(n & 1) || n < 2 ) return n == 2;
- 
-/* comparing p*p <= n can overflow */
-for (p = 3; p <= n/p; p += 2)
-if (!(n % p)) return 0;
-return 1;
+if (n <=1) {
+return false;
+}
+
+if (n==2) {
+return true;
+}
+
+unsigned sqrtn = sqrt(n);
+ for (unsigned int i = 2; i <= sqrtn; ++i) {
+ if (n % i == 0) {
+ return false;
+}
+}
+return true;
+}
+
+/**
+ * http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test 
+ */
+static bool miller_rabin_test(unsigned int n) {
+switch (n) {
+case 0:
+case 1:
+return false;
+}
+
+unsigned s = ZerosR(n-1);
+unsigned d = (n-1) >> s;
+
+// test 3-times
+for (int k=0;k<3;k++){
+unsigned a = rand()%(n-1) + 1;
+
+unsigned x = Exp(a, d, n);
+if (x == 1 || x == n - 1) {
+continue;
+}
+
+for (unsigned i=0;i<s;i++) {
+x = Exp(x, 2, n);
+if (x == 1) return false;
+if (x == n-1) break;
+}
+}
+
+return true;
+}
+
+/**
+ * mixed implementation
+ */
+static bool is_prime(unsigned int n) {
+if (miller_rabin_test(n)) {
+return test_prime(n);
+}
+return false;
 }
 }
 

src/imath_demo.cpp

@@ -8,4 +8,7 @@ int main() {
 printf("%d^%d mod %d = %d\n", i,j,997, alg::Exp(i, j, 997));
 }
 }
+
+unsigned int n = 104;
+printf("%x has %d zeros on the right side\n", n, alg::ZerosR(n));
 }

src/prime_test.cpp

@@ -1,12 +1,18 @@
 #include <stdio.h>
 #include <prime.h>
+#include <time.h>
 
 using namespace alg;
 int main()
 {
-int i;
-for (i = 0; i<100; i++) {
-printf("number: %d, prime %d\n", i, is_prime(i));
+unsigned i;
+int count = 0;
+time_t t1 = time(NULL);
+for (i = 1000000000; i<1000100000; i++) {
+if (is_prime(i)) count++;
 }
+
+time_t t2 = time(NULL);
+printf("found %d primes, cost %ld secs\n", count, t2-t1);
 return 0;
 }