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?

If you want to do it yourself, the simple and naive solution is to use arrays of digits, and do the math as taught in basis school digit by digit. Otherwise there are nice libraries for arbitrary precision numbers and arithmetic. — Some programmer dude, Apr 26 '22 at 09:07
If you use integer, then there isn't a chance you can represent `10^1053` in any of the built in types. (As a very naive guess I would say you'd need somewhere like 4kb of precision for an int to suffice). `double` isn't great either since at that point I'd expect that's precision to degrade as well. — Lala5th, Apr 26 '22 at 09:10
The `1%m` looks VERY suspect. That's just `1` for every value of m. Did you mean "a equal to 1, modulo m"? Because that's `(a-1) %m == 0`. — MSalters, Apr 26 '22 at 09:12
First of all `(1 % m)`? I'm sure you meant `double(pow(a, t)) % m == 1` right? Probably cast to `int`, not `double` btw. If that's the case then you don't need to calculate `pow` and do modulo then. This will exceed range. You should apply one of the [modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) algorithms instead. — freakish, Apr 26 '22 at 09:12
@MSalters well, it did work for python. And for small numbers as i said in the beginning of the post. I wrote that code using wiki's algorithm — Max Buldakov, Apr 26 '22 at 09:16
@freakish yeah, made typo for `int` and `double` :D Tried changing `1 % m` to `pow % m == 1` and it didn't make sense. I guess i understood Miller-Rabin algorithm wrong and did mistake in a code then. — Max Buldakov, Apr 26 '22 at 09:21
@MaxBuldakov The Wiki page for Miller–Rabin primality test does say what @MSalters assumes that you need `a^t = 1 (mod n)`, which in code translates to `a^t mod n = 1`. The `(mod n)` in the definition only signals that you are working on a specific modular arithmetic, not that you need the modulo there (the original expression reads as `a^t` congruent to `1` under modulo `n`). Also since you are woking under some modular arithmetic it might be better to implement a modular arithmetic pow function somehow — Lala5th, Apr 26 '22 at 09:21
You're using `pow()` which always returns a floating point value. Floating point variables, by their nature, do not represent all integral values exactly, so using floating point functions/operations to produce integral results will - for some set of arguments - introduce numerical error. Your claim that it worked in Python is irrelevant - C++ is not Python, so some things work differently in C++ than Python - and this is one of them. Work out a way (e.g. algorithm) to raise an integral value to an integral power and produce an integral result without intermediate use of floating point — Peter, Apr 26 '22 at 09:24
The line `int y = int(pow(a, t));` is useless, you call again `pow` on the following line. — Fabio says Reinstate Monica, Apr 26 '22 at 09:31
@FabiosaysReinstateMonica did it just to see what value it gets while did debug — Max Buldakov, Apr 26 '22 at 09:56

freakish · Answer 1 · 2022-04-27T07:24:13.830

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.

Thanks for explanation. I got the main idea of modular exponentiation. Will try your fixes and think about other parts of code. Thank you a lot! — Max Buldakov, Apr 26 '22 at 09:38
@MaxBuldakov: when you are ready to improve the efficiency of your modular_pow method consider [these algorithms](https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method) which are only slightly more complicated but potentially much faster. — President James K. Polk, Apr 26 '22 at 19:11

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

1 Answers1