Is probable prime

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

include/functions/math.hpp

@@ -7,9 +7,13 @@
 #ifndef BIG_INT_MATH_FUNCTIONS_HPP
 #define BIG_INT_MATH_FUNCTIONS_HPP
 
-#include <string>
+
+#include <tuple>
+#include <random>
+#include <functional>
 
 #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<int> 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)

include/functions/utility.hpp

@@ -8,7 +8,7 @@
 #define BIG_INT_UTILITY_FUNCTIONS_HPP
 
 #include <tuple>
-
+#include <math.h>
 
 /*
     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<int, int> 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

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);
+
+}