Primality test and helper functions

include/BigInt.hpp

@@ -102,6 +102,10 @@ class BigInt {
 
         // Random number generating functions:
         friend BigInt big_random(size_t);
+        friend BigInt n_random(std::string);
+
+        //Primality test
+        bool is_probable_prime(size_t);
 };
 
 #endif  // BIG_INT_HPP

include/functions/math.hpp

@@ -10,7 +10,7 @@
 #include <string>
 
 #include "functions/conversion.hpp"
-
+#include "functions/random.hpp"
 
 /*
     abs
@@ -64,6 +64,34 @@ BigInt pow(const BigInt& base, int exp) {
     return result * result_odd;
 }
 
+/* pow (BigInt, BigInt)
+   ----------------
+   Returns a BigInt equal to base^exp where exp is a BigInt
+*/
+
+BigInt pow(const BigInt& base, BigInt exp) {
+    if (exp < 0) {
+        if (base == 0)
+            throw std::logic_error("Cannot divide by zero");
+        return abs(base) == 1 ? base : 0;
+    }
+    if (exp == 0) {
+        if (base == 0)
+            throw std::logic_error("Zero cannot be raised to zero");
+        return 1;
+    }
+
+    BigInt result = base, result_odd = 1;
+    while (exp > 1){
+        if (exp % 2 == 1)
+            result_odd *= result;
+        result *= result;
+        exp /= 2;
+    }
+
+    return result * result_odd;
+}
+
 
 /*
     pow (Integer)
@@ -247,5 +275,66 @@ BigInt lcm(const std::string& num1, const BigInt& num2){
     return lcm(BigInt(num1), num2);
 }
 
+/*
+    is_probable_prime(size_t)
+    ------------------------
+    Uses the Miller-Rabin primality test to return if the BigInt
+    is prime with probablity ( 1 - (4^-certainty) ) * 100%
+*/
+
+bool BigInt::is_probable_prime(size_t certainty){
+    //treats 1, 2, and 3 as prime numbers
+    if (*this == BigInt(1) || *this == BigInt(2) || *this == BigInt(3)){
+        return true;
+    }
+
+    //even numbers cannot be prime
+    if (*this % BigInt(2) == 0){
+        return false;
+    }
+
+    const BigInt maxRand = *this - 2; //we later choose random value between 0 to n-1
+    const BigInt one = 1;
+    const BigInt two = 2;
+    BigInt randNum;
+
+    //need to compute d and r such that d*2^r = n - 1.  where n = this
+    BigInt d;
+    BigInt x;
+    int r;
+    d = maxRand;
+    ++d;
+    r = 0;
+    int continueWhile = 1;
+    while( d % two == 0){
+        ++r;
+        d /= two;
+    }
+    while ( certainty-- > 0 ){
+        //pick a random number 
+        randNum = n_random(maxRand.value);
+
+        x = pow(randNum, d);
+        x = x % *this;
+
+        if (x == one || x == *this - one){
+            continue;
+        }    
+        continueWhile = 0;
+        for( int i=0; i < r-1; i++){
+            x = pow(x, 2) % *this;
+            if (x == *this-one){
+                continueWhile = 1;
+                break;
+            }
+        }
+        if(continueWhile) continue;
+        return false;
+
+    }
+    return true;
+}
+
+
 
 #endif  // BIG_INT_MATH_FUNCTIONS_HPP

include/functions/random.hpp

@@ -9,6 +9,8 @@
 
 #include <random>
 #include <climits>
+#include <stdlib.h>
+#include <string>
 
 #include "BigInt.hpp"
 
@@ -44,5 +46,54 @@ BigInt big_random(size_t num_digits = 0) {
     return big_rand;
 }
 
+/*
+    n_random (n)
+    --------------
+    Returns a random BigInt from the range of 2 to n inclusive
+*/
+
+BigInt n_random(std::string val){
+    std::string randVal = "";
+    int pushVal;
+    int freeRandom = 0; //to check if current digit is constrained
+    for(size_t i=0; i<val.size(); i++){
+        //i.e if all the previous randomly generated most significant digit
+        //was equal to the previous most significant digit of n
+        if(freeRandom == 0){
+            pushVal = std::rand() % (val[i] - '0' + 1);
+            randVal += std::to_string(pushVal);
+            if(pushVal < (val[i] - '0')){
+                freeRandom = 1;
+            }
+        }
+        else{
+            pushVal = std::rand() % 10;
+            randVal += std::to_string(pushVal);
+        }
+    }
+    //delete 0's if it is in the most significant digit
+    while(*(randVal.begin()) == '0' && randVal.length() > 1){
+        randVal.erase(randVal.begin());
+    }
+
+    BigInt randomNum = randVal; 
+
+    //if a 0 or 1 has been generated
+    if(randomNum == 0 || randomNum == 1){
+        if(val.length() == 1){
+            pushVal = rand() % (val[0] - 1) + 2;
+        }
+        else{
+            pushVal = rand() % 7 + 2;
+        }
+        randVal += std::to_string(pushVal);
+        randVal.erase(randVal.begin());
+        randomNum = randVal;
+    }
+
+
+
+    return randomNum;
+}
 
 #endif  // BIG_INT_RANDOM_FUNCTIONS_HPP

test/functions/math.cpp

@@ -416,3 +416,83 @@ 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);
+
+}
+
+
+TEST_CASE("Even number cases for is_probable_prime()", "[functions][math][is_probable_prime]"){
+    BigInt num = 22;
+    REQUIRE(num.is_probable_prime(25) == 0);
+    num = "2342349048751934651982342934622123757987072";
+    REQUIRE(num.is_probable_prime(25) == 0);
+    num = "489275109347659813465918246912837641923746189234619268541924";
+    REQUIRE(num.is_probable_prime(25) == 0);
+
+}
+
+TEST_CASE("Definitely composite numbers for is_probable_prime()", "[functions][math][is_probable_prime]"){
+    BigInt num = 576308207413;
+
+//Commented out due to very significant runtime (hours)
+/*    REQUIRE(num.is_probable_prime(2) == 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);
+*/
+    num = "33";
+    REQUIRE(num.is_probable_prime(25) == 0);
+
+    num = 500067;
+//    REQUIRE(num.is_probable_prime(5) == 0);
+
+    num = 20097;
+//    REQUIRE(num.is_probable_prime(5) == 0);
+}
+
+
+TEST_CASE("Prime numbers for is_probable_prime()", "[functions][math][is_probable_prime]"){
+    BigInt num = 4361161843811;
+
+//Commented out due to very significant runtime (hours)
+/*    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);
+*/    
+    num = 31;
+    REQUIRE(num.is_probable_prime(25) == 1);
+
+    num = 24862048;
+//    REQUIRE(num.is_probable_prime(5) == 1);
+
+    num = 500057;
+//    REQUIRE(num.is_probable_prime(5) == 1);
+
+    num = 19069;
+    REQUIRE(num.is_probable_prime(10) == 1);
+
+}

test/functions/random.cpp

@@ -4,9 +4,9 @@
 #include "functions/conversion.hpp"
 #include "functions/random.hpp"
 #include "operators/io_stream.hpp"
-
+#include "operators/relational.hpp"
 #include "third_party/catch.hpp"
-
+#include <string>
 
 TEST_CASE("Generate random BigInts", "[functions][random]") {
     for (int i = 0; i < 50; i++)
@@ -22,3 +22,56 @@ TEST_CASE("Generate random BigInts having a specified number of digits",
         REQUIRE(big_random(num_digits).to_string().size() == num_digits);
     }
 }
+
+TEST_CASE("Test n_random", "[functions][random]"){
+    std::string maxRand = "10";
+    BigInt randVal;
+    int zeroA = 0;
+    int oneA = 0;
+    int twoA = 0;
+    int threeA = 0;
+    int fourA = 0;
+    int fiveA = 0;
+    int sixA = 0;
+    int sevenA = 0;
+    int eightA = 0;
+    int nineA = 0;
+    int tenA = 0;
+    for(int i=0; i<300; i++){
+        randVal = n_random(maxRand);
+        REQUIRE(randVal <= maxRand);
+        REQUIRE(randVal >= 2);
+        if(randVal == 0) zeroA = 1;
+        if(randVal == 1) oneA = 1;
+        if(randVal == 2) twoA = 1;
+        if(randVal == 3) threeA = 1;
+        if(randVal == 4) fourA = 1;
+        if(randVal == 5) fiveA = 1;
+        if(randVal == 6) sixA = 1;
+        if(randVal == 7) sevenA = 1;
+        if(randVal == 8) eightA = 1;
+        if(randVal == 9) nineA = 1;
+        if(randVal == 10) tenA = 1;
+    }
+    REQUIRE(zeroA == 0);
+    REQUIRE(oneA == 0);
+    REQUIRE(twoA == 1);
+    REQUIRE(threeA == 1);
+    REQUIRE(fourA == 1);
+    REQUIRE(fiveA == 1);
+    REQUIRE(sixA == 1);
+    REQUIRE(sevenA == 1);
+    REQUIRE(eightA == 1);
+    REQUIRE(nineA == 1);
+    REQUIRE(tenA == 1);
+    maxRand = "239";
+    for(int i=0; i<300; i++){
+        randVal = n_random(maxRand);
+        REQUIRE(randVal <= maxRand );
+        REQUIRE(randVal >= 2);
+    }
+
+
+
+
+}