Miller–Rabin primality test: Difference between revisions
New post.
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
(New post.) |
||
Line 903:
}
}</syntaxhighlight>
=={{header|C++}}==
This test is deterministic for n < 3'474'749'660'383
<syntaxhighlight lang="c++">
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>
std::vector<uint32_t> small_primes{2, 3};
uint64_t add_modulus(const uint64_t& a, const uint64_t& b, const uint64_t& modulus) {
uint64_t am = ( a < modulus ) ? a : a % modulus;
if ( b == 0 ) {
return am;
}
uint64_t bm = ( b < modulus ) ? b : b % modulus;
uint64_t b_from_m = modulus - bm;
if ( am >= b_from_m ) {
return am - b_from_m;
}
return am + bm;
}
uint64_t multiply_modulus(const uint64_t& a, const uint64_t& b, const uint64_t& modulus) {
uint64_t am = ( a < modulus ) ? a : a % modulus;
uint64_t bm = ( b < modulus ) ? b : b % modulus;
if ( bm > am ) {
std::swap(am, bm);
}
uint64_t result;
while ( bm > 0 ) {
if ( ( bm & 1) == 1 ) {
result = add_modulus(result, am, modulus);
}
am = ( am << 1 ) - ( am >= ( modulus - am ) ? modulus : 0 );
bm >>= 1;
}
return result;
}
uint64_t exponentiation_modulus(const uint64_t& base, const uint64_t& exponent, const uint64_t& modulus) {
uint64_t b = base;
uint64_t e = exponent;
uint64_t result = 1;
while ( e > 0 ) {
if ( ( e & 1 ) == 1 ) {
result = multiply_modulus(result, b, modulus);
}
e >>= 1;
b = multiply_modulus(b, b, modulus);
}
return result;
}
bool is_composite(const uint32_t& a, const uint64_t& d, const uint64_t& n, const uint32_t& s) {
if ( exponentiation_modulus(a, d, n) == 1 ) {
return false;
}
for ( uint64_t i = 0; i < s; ++i ) {
if ( exponentiation_modulus(a, pow(2, i) * d, n) == n - 1 ) {
return false;
}
}
return true;
}
bool composite_test(const std::vector<uint32_t>& primes, const uint64_t& d, const uint64_t& n, const uint32_t& s) {
for ( const uint32_t& prime : primes ) {
if ( is_composite(prime, d, n, s) ) {
return true;
}
}
return false;
}
bool is_prime(const uint64_t& n) {
if ( n == 0 || n == 1 ) {
return false;
}
if ( std::find(small_primes.begin(), small_primes.end(), n) != small_primes.end() ) {
return true;
}
if ( std::any_of(small_primes.begin(), small_primes.end(), [n](uint32_t p) { return n % p == 0; }) ) {
return false;
}
uint64_t d = n - 1;
uint32_t s = 0;
while ( ! d % 2 ) {
d >>= 1;
s++;
}
if ( n < 1'373'653 ) {
return composite_test({ 2, 3 }, d, n, s);
}
if ( n < 25'326'001 ) {
return composite_test({ 2, 3, 5 }, d, n, s);
}
if ( n < 118'670'087'467 ) {
if ( n == 3'215'031'751 ) {
return false;
}
return composite_test({ 2, 3, 5, 7 }, d, n, s);
}
if ( n < 2'152'302'898'747 ) {
return composite_test({ 2, 3, 5, 7, 11 }, d, n, s);
}
if ( n < 3'474'749'660'383 ) {
return composite_test({ 2, 3, 5, 7, 11, 13 }, d, n, s);
}
if ( n < 341'550'071'728'321 ) {
return composite_test({ 2, 3, 5, 7, 11, 13, 17 }, d, n, s);
}
const std::vector<uint32_t> test_primes(small_primes.begin(), small_primes.begin() + 16);
return composite_test(test_primes, d, n, s);
}
void create_small_primes() {
for ( uint32_t i = 5; i < 1'000; i += 2 ) {
if ( is_prime(i) ) {
small_primes.emplace_back(i);
}
}
}
int main() {
create_small_primes();
for ( const uint64_t number : { 1'234'567'890'123'456'733, 1'234'567'890'123'456'737 } ) {
std::cout << "is_prime(" << number << ") = " << std::boolalpha << is_prime(number) << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
is_prime(1234567890123456733) = false
is_prime(1234567890123456737) = true
</pre>
=={{header|Clojure}}==
| |||