C++ and way-too-big numbers (Miller-Rabin primality test; calculating a power with pow() and then applying the modulo operator)

Question

I had to write a code, that determines whether a number is prime or not via Miller-Rabin primality test.

The code here:

#include <iostream>
#include <time.h>
#include <cmath>
#include <stdio.h>  

using namespace std;

int main()
{
    srand(time(0));
    int m, r, st, a, x, t,st2;
    cout << "Введите число m: ";
    cin >> m;
    cout << "Введите число r: ";
    cin >> r;
    t = m - 1;
    st = 0;

    while (t % 2 == 0)
    {
        st = st + 1;
        t = t / 2;
    }

    for (int i = 0; i < r; i++)
    {
        a = rand() % (m - 1) + 1;
        int y = int(pow(a, t));
        if ((int(pow(a, t))) == (1 % m))
        {
            cout << "Number " << m << " is probably prime" << endl;
            exit(0);
        }
        for (int j = 0; j < st - 1; j++)
        {
            st2 = int(pow(2, i));
            x = int(pow(a, st2 * t));
            if (x == -1 % m)
            {
                continue;
            }
            else
            {
                cout << "Number " << m << " is composite" << endl;
                exit(0);
            }
        }
    }
    cout << "Number " << m << " probably prime" << endl;
}

It did pretty well with small numbers like 221, 12 etc. But i had a bunch of numbers and on one of them, number 2633, something went wrong: this number is prime, but program says it's not. I did a line by line run and noticed it is stucks on line if ((double(pow(a, t))) == (1 % m)) The problem here is that it needs to evaluate power 329 for some random number between 1 and 2632. For example we'll take 1595 (one of the numbers taken randomly during debug). Now we evaluate int y = int(pow(a, t)) (took it for debug to see what it evaluates) and working with Windows Calculator we're getting this number:

5,1081795295996211569794748579036e+1053

and the program gets -2147483648... Every time...

I don't know exactly where is the problem - in pow function or in number types. Also tried double instead of int but it didn't work (i guess because 5 * 10^1053 is way bigger than double's 1.5 * 10^308)

How do i solve that problem and get it to work?

Answer 1

This

if ((int(pow(a, t))) == (1 % m))

is very wrong. First of all 1 % m is just 1 for any m > 1. And indeed, Miller–Rabin needs

if ( int(pow(a, t)) % m == 1 )

kind of tests. This is how you perform modular arithmetic in programming languages and you should modify your algorithm to use this approach everywhere, e.g. (x == -1 % m) should be replaced with (x % m == m-1) (because % always returns a nonnegative value for nonnegative input and -1 corresponds to m-1 in modular arithmetic).

Secondly, yes, this will go out of range very quickly. Calculating a power and then applying modulo operator is not a valid approach for big numbers. What you should do is to implement one of modular exponentiation algorithms instead. The simplest solution is

int modular_pow(int base, int exponent, int modulus) {
    if (modulus == 1) return 0;
    
    int c = 1;
    for (int i = 0; i < exponent; i++)
        c = (1LL * c * base) % modulus;  // 1LL to avoid overflow

    return c;
}

Not very efficient but at least will do the job. And then you check

if (modular_pow(a,t,m) == 1)

Thirdy, you use int(pow(...)) conversion, which can lead to numerical errors, since pow returns a floating point value. But you need exact values to perform exact modular arithmetic. Do not use pow in such cases. Read more about integer pow here: The most efficient way to implement an integer based power function pow(int, int) Although let me remind you that modular exponentiation is what you really need here.

It is possible that there are other issues with your Miller-Rabin's implementation, but these three are the biggest I see.