#include <iostream>
#include <algorithm>
template <typename R, typename S, typename T>
T POW(R base, S exponent, const T mod){
T result = 1;
while (exponent){
if (exponent & 1)
result = (result * base) % mod;
exponent >>= 1;
base = (base * base) % mod;
}
return result;
}
// used uint64_t to prevent overflow, but only testing with small numbers for now
bool MillerRabin_FIPS186(uint64_t w, unsigned int iterations = 50){
srand(time(0));
unsigned int a = 0;
uint64_t W = w - 1; // dont want to keep calculating w - 1
uint64_t m = W;
while (!(m & 1)){
m >>= 1;
a++;
}
// skipped getting wlen
// when i had this function using my custom arbitrary precision integer class,
// and could get len(w), getting it and using it in an actual RBG
// made no difference
for(unsigned int i = 0; i < iterations; i++){
uint64_t b = (rand() % (W - 3)) + 2; // 2 <= b <= w - 2
uint64_t z = POW(b, m, w);
if ((z == 1) || (z == W))
continue;
else
for(unsigned int j = 1; j < a; j++){
z = POW(z, 2, w);
if (z == W)
break;
if (z == 1)
return 0;// Composite
}
}
return 1;// Probably Prime
}
int main() {
std::cout << MillerRabin_FIPS186(33) << std::endl;
std::cout << MillerRabin_FIPS186(35) << std::endl;
std::cout << MillerRabin_FIPS186(37) << std::endl;
std::cout << MillerRabin_FIPS186(39) << std::endl;
std::cout << MillerRabin_FIPS186(45) << std::endl;
std::cout << MillerRabin_FIPS186(49) << std::endl;
return 0;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8YWxnb3JpdGhtPgoKdGVtcGxhdGUgPHR5cGVuYW1lIFIsIHR5cGVuYW1lIFMsIHR5cGVuYW1lIFQ+ClQgUE9XKFIgYmFzZSwgUyBleHBvbmVudCwgY29uc3QgVCBtb2QpewogICAgVCByZXN1bHQgPSAxOwogICAgd2hpbGUgKGV4cG9uZW50KXsKICAgICAgICBpZiAoZXhwb25lbnQgJiAxKQogICAgICAgICAgICByZXN1bHQgPSAocmVzdWx0ICogYmFzZSkgJSBtb2Q7CiAgICAgICAgZXhwb25lbnQgPj49IDE7CiAgICAgICAgYmFzZSA9IChiYXNlICogYmFzZSkgJSBtb2Q7CiAgICB9CiAgICByZXR1cm4gcmVzdWx0Owp9CgovLyB1c2VkIHVpbnQ2NF90IHRvIHByZXZlbnQgb3ZlcmZsb3csIGJ1dCBvbmx5IHRlc3Rpbmcgd2l0aCBzbWFsbCBudW1iZXJzIGZvciBub3cKYm9vbCBNaWxsZXJSYWJpbl9GSVBTMTg2KHVpbnQ2NF90IHcsIHVuc2lnbmVkIGludCBpdGVyYXRpb25zID0gNTApewogICAgc3JhbmQodGltZSgwKSk7CiAgICB1bnNpZ25lZCBpbnQgYSA9IDA7CiAgICB1aW50NjRfdCBXID0gdyAtIDE7IC8vIGRvbnQgd2FudCB0byBrZWVwIGNhbGN1bGF0aW5nIHcgLSAxCiAgICB1aW50NjRfdCBtID0gVzsKICAgIHdoaWxlICghKG0gJiAxKSl7CiAgICAgICAgbSA+Pj0gMTsKICAgICAgICBhKys7CiAgICB9CgogICAgLy8gc2tpcHBlZCBnZXR0aW5nIHdsZW4KICAgIC8vIHdoZW4gaSBoYWQgdGhpcyBmdW5jdGlvbiB1c2luZyBteSBjdXN0b20gYXJiaXRyYXJ5IHByZWNpc2lvbiBpbnRlZ2VyIGNsYXNzLAogICAgLy8gYW5kIGNvdWxkIGdldCBsZW4odyksIGdldHRpbmcgaXQgYW5kIHVzaW5nIGl0IGluIGFuIGFjdHVhbCBSQkcgCiAgICAvLyBtYWRlIG5vIGRpZmZlcmVuY2UgCgogICAgZm9yKHVuc2lnbmVkIGludCBpID0gMDsgaSA8IGl0ZXJhdGlvbnM7IGkrKyl7CiAgICAgICAgdWludDY0X3QgYiA9IChyYW5kKCkgJSAoVyAtIDMpKSArIDI7IC8vIDIgPD0gYiA8PSB3IC0gMgogICAgICAgIHVpbnQ2NF90IHogPSBQT1coYiwgbSwgdyk7CiAgICAgICAgaWYgKCh6ID09IDEpIHx8ICh6ID09IFcpKQogICAgICAgICAgICBjb250aW51ZTsKICAgICAgICBlbHNlCiAgICAgICAgICAgIGZvcih1bnNpZ25lZCBpbnQgaiA9IDE7IGogPCBhOyBqKyspewogICAgICAgICAgICAgICAgeiA9IFBPVyh6LCAyLCB3KTsKICAgICAgICAgICAgICAgIGlmICh6ID09IFcpCiAgICAgICAgICAgICAgICAgICAgYnJlYWs7CiAgICAgICAgICAgICAgICBpZiAoeiA9PSAxKQogICAgICAgICAgICAgICAgICAgIHJldHVybiAwOy8vIENvbXBvc2l0ZQogICAgICAgICAgICB9CiAgICB9CiAgICByZXR1cm4gMTsvLyBQcm9iYWJseSBQcmltZQp9CgppbnQgbWFpbigpIHsKCXN0ZDo6Y291dCA8PCBNaWxsZXJSYWJpbl9GSVBTMTg2KDMzKSA8PCBzdGQ6OmVuZGw7CglzdGQ6OmNvdXQgPDwgTWlsbGVyUmFiaW5fRklQUzE4NigzNSkgPDwgc3RkOjplbmRsOwoJc3RkOjpjb3V0IDw8IE1pbGxlclJhYmluX0ZJUFMxODYoMzcpIDw8IHN0ZDo6ZW5kbDsKCXN0ZDo6Y291dCA8PCBNaWxsZXJSYWJpbl9GSVBTMTg2KDM5KSA8PCBzdGQ6OmVuZGw7CglzdGQ6OmNvdXQgPDwgTWlsbGVyUmFiaW5fRklQUzE4Nig0NSkgPDwgc3RkOjplbmRsOwoJc3RkOjpjb3V0IDw8IE1pbGxlclJhYmluX0ZJUFMxODYoNDkpIDw8IHN0ZDo6ZW5kbDsKCXJldHVybiAwOwp9