diff --git a/include/BigInt.hpp b/include/BigInt.hpp index e3be86c..1cb3774 100644 --- a/include/BigInt.hpp +++ b/include/BigInt.hpp @@ -102,6 +102,9 @@ class BigInt { // Random number generating functions: friend BigInt big_random(size_t); + + //Math + bool is_probable_prime(size_t); }; #endif // BIG_INT_HPP diff --git a/include/functions/math.hpp b/include/functions/math.hpp index 9dee74e..ceb6bc4 100644 --- a/include/functions/math.hpp +++ b/include/functions/math.hpp @@ -7,9 +7,13 @@ #ifndef BIG_INT_MATH_FUNCTIONS_HPP #define BIG_INT_MATH_FUNCTIONS_HPP -#include + +#include +#include +#include #include "functions/conversion.hpp" +#include "functions/utility.hpp" /* @@ -153,6 +157,61 @@ BigInt gcd(const BigInt &num1, const BigInt &num2){ return abs_num1; } +/* + BigInt::is_probable_prime(size_t) + --------------------------------- + Returns boolean representing the primality of the BigInt object that this method refers + to using Rabin-Miller Primality testing. + NOTE: This test is probablisitic so it is not failproof. A composite number passes the test + for at most 1/4 of the possible bases. This is true for each iteration of the test, so the + Probability a composite number passes is 1/4^N or less, where N is equal to certainty + */ + +bool BigInt::is_probable_prime(size_t certainty) { + // The BigInt object that this method refers to will be referred to as n in comments + + if (*this <= BigInt(3) || *this == BigInt(5)) { + return true; + } + + std::default_random_engine generator; + std::uniform_int_distribution distribution(2, this->to_long_long()-2); + auto get_random_number = std::bind(distribution, generator); + + const BigInt ONE = 1; + const BigInt TWO = 2; + BigInt rand_num; + while (certainty-- > 0) { + // 1 <= rand_num < n + rand_num = get_random_number(); + // If there exists an x > 1 such that rand_num % x == 0 && n % x == 0 + // then n is composite + if (gcd(rand_num, *this) != ONE) { + return false; + } + int s, m; + // Calculates needed variables that fit the equation + // n - 1 = 2^s*m such that s >= 1 and m is odd + std::tie(s, m) = calculate_vars(this->to_long_long()); + // x = rand_num^m%n + BigInt x = pow(rand_num, m)%*this; + if (x == ONE || x == *this-ONE) { + continue; + } + for (int i = 0; i < s-1; i++) { + // x = x^2%n + x = pow(x, 2)%*this; + if (x == 1) { + return false; + } else if (x == *this-ONE) { + continue; + } + return false; + } + } + return true; +} + /* gcd(BigInt, Integer) diff --git a/include/functions/utility.hpp b/include/functions/utility.hpp index a26ed45..33d7d57 100644 --- a/include/functions/utility.hpp +++ b/include/functions/utility.hpp @@ -8,7 +8,7 @@ #define BIG_INT_UTILITY_FUNCTIONS_HPP #include - +#include /* is_valid_number @@ -110,4 +110,36 @@ bool is_power_of_10(const std::string& num){ return true; // first digit is 1 and the following digits are all 0 } + +/* + calculate_vars + -------------- + Calculates the s and m variables needed to complete Rabin Miller Primality + test. + + s and m come from the formula n-1 = 2^s*m, where m is odd and s >= 1 +*/ + +std::tuple calculate_vars(long long n) { + int s = 1; + int m = 1; + const long long sentinel_value = n-1; + long long calculated_value = 2; + + while (sentinel_value >= calculated_value) { + if (sentinel_value > calculated_value) { + s++; + long long spower = pow(2, s); + if (spower > sentinel_value) { + return std::make_tuple(0,0); + } + m = sentinel_value/spower; + calculated_value = spower*m; + } else if (sentinel_value == calculated_value) { + return std::make_tuple(s, m); + } + } + return std::make_tuple(s, m); +} + #endif // BIG_INT_UTILITY_FUNCTIONS_HPP diff --git a/test/functions/math.cpp b/test/functions/math.cpp index c769652..3a40dfb 100644 --- a/test/functions/math.cpp +++ b/test/functions/math.cpp @@ -91,7 +91,7 @@ TEST_CASE("Base cases for pow()", "[functions][math][pow]") { } TEST_CASE("pow() with BigInt base", "[functions][math][pow]") { - BigInt num = 11; + BigInt num = 11; REQUIRE(pow(num, 9) == 2357947691); num = -27; REQUIRE(pow(num, 16) == "79766443076872509863361"); @@ -382,3 +382,47 @@ TEST_CASE("lcm()of big integers", "[functions][math][lcm][big]") { "533220692127153044311875258011747917053108027629278373174251200266431" "428784066739966"); } + +TEST_CASE("Base cases for is_probable_prime()", "[functions][math][is_probable_prime]") { + BigInt num = 1; + REQUIRE(num.is_probable_prime(25) == 1); + num = 2; + REQUIRE(num.is_probable_prime(25) == 1); + num = 3; + REQUIRE(num.is_probable_prime(25) == 1); + num = 5; + REQUIRE(num.is_probable_prime(25) == 1); +} + +TEST_CASE("is_probable_prime() for big integers true", "[functions][math][is_probable_prime]") { + BigInt num = 4361161843811; + REQUIRE(num.is_probable_prime(25) == 1); + num = 91584398720031; + REQUIRE(num.is_probable_prime(25) == 1); + num = "54362229927468991056799869539182953179604007"; + REQUIRE(num.is_probable_prime(25) == 1); + num = "1141606828476848812192797260322842016771684147"; + REQUIRE(num.is_probable_prime(25) == 1); + num = "237082482904158189833801188468727382999221896206963750677"; + REQUIRE(num.is_probable_prime(25) == 1); + num = "4978732140987321986509824957843275042983659820346238764217"; + REQUIRE(num.is_probable_prime(25) == 1); +} + +TEST_CASE("is_probable_prime() for big integers false", "[functions][math][is_probable_prime]") { + BigInt num = 576308207413; + REQUIRE(num.is_probable_prime(25) == 0); + num = 648273642634986; + REQUIRE(num.is_probable_prime(25) == 0); + num = "328964398726983264982"; + REQUIRE(num.is_probable_prime(25) == 0); + num = "4079327665427094820942557"; + REQUIRE(num.is_probable_prime(25) == 0); + num = "879654387682647825646328764"; + REQUIRE(num.is_probable_prime(25) == 0); + num = "98732243986019286982046325298743"; + REQUIRE(num.is_probable_prime(25) == 0); + num = "589658224987973265974369876397863298796328749"; + REQUIRE(num.is_probable_prime(25) == 0); + +}