From f6b6c4c8694cd398b67ac0c2b4ad4fdf0b782c58 Mon Sep 17 00:00:00 2001 From: JBatzill Date: Tue, 1 Dec 2015 13:50:20 +0100 Subject: Added miller rabin implementation tested with big primes <= 10^18 from wikipedia and solved problem: https://open.kattis.com/problems/primes2 --- math/miller_rabin.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 math/miller_rabin.cpp (limited to 'math') diff --git a/math/miller_rabin.cpp b/math/miller_rabin.cpp new file mode 100644 index 0000000..ad8c163 --- /dev/null +++ b/math/miller_rabin.cpp @@ -0,0 +1,20 @@ +//theoretical: n < 318,665,857,834,031,151,167,461 ( > 10^23) +//but: n ~<= 10^18 (because of MAX(ll)) +//O(logn) +bool isPrime(ll n) { + if(n == 2) return true; + if(n < 2 || n % 2 == 0) return false; + ll d=n-1,j=0; + while(d % 2 == 0) d >>= 1, j++; + for(int a = 2; a <= min((ll)37, n-1); a++) { + ll v = pow_mod(a, d, n); + if(v == 1 || v == n-1) continue; + for(int i = 1; i <= j; i++) { + v = mult_mod(v, v, n); + if(v == n-1 || v <= 1) break; + } + + if(v != n-1) return false; + } + return true; +} -- cgit v1.2.3 From 30ea30f2a7029c143aeb823a8f1a77c9f20d9e48 Mon Sep 17 00:00:00 2001 From: JBatzill Date: Tue, 1 Dec 2015 14:06:38 +0100 Subject: added new pow_mod method and mult_pow! Improves the old version since multiplication was able to overflow more easily! --- math/modExp.cpp | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) (limited to 'math') diff --git a/math/modExp.cpp b/math/modExp.cpp index f738268..863ff4e 100644 --- a/math/modExp.cpp +++ b/math/modExp.cpp @@ -1,7 +1,17 @@ -ll modPow(ll b, ll e, ll p) { - if (e == 0) return 1; - if (e == 1) return b; - ll half = modPow(b, e / 2, p), res = (half * half) % p; - if (e & 1) res *= b; res %= p; - return res; -} \ No newline at end of file +//0<=a,b <=n and n <= MAX(ll)/2 +ll mult_mod(ll a, ll b, ll n) { + if(a == 0 || b == 0) return 0; + if(b == 1) return a % n; + + if(b % 2 == 1) return (a + mult_mod(a, b-1, n)) % n; + else return mult_mod((a + a) % n, b / 2, n); +} + +//0<=a,b<=n and n <= MAX(ll)/2 +ll pow_mod(ll a, ll b, ll n) { + if(b == 0) return 1; + if(b == 1) return a % n; + + if(b % 2 == 1) return mult_mod(pow_mod(a, b-1, n), a, n); + else return pow_mod(mult_mod(a, a, n), b / 2, n); +} -- cgit v1.2.3