#include <iostream>
#include <vector>
#include <bitset>

using namespace std;

int repeated_squaring_modular_exponentiation(int input, int power, int modulo)
{
    std::bitset<32> bits(power);
    std::vector<int> powers{ input };
    int pw = 2;
    int temp = input;
    while (pw <= power)
    {
        temp = (temp * temp) % modulo;
        powers.push_back(temp);
        pw *= 2;
    }

    int result = 1;
    for (unsigned int i = 0; i < powers.size(); ++i)
    {
        if (!bits.test(i))
            continue;

        result = (result * powers[i]) % modulo;
    }

    return result;
}

bool miller_rabin(int n, int s, int t, int b)
{
    // Two possible cases:
    // b^t = 1 mod n
    // or
    // there is a j: 0 <= j < s such that b^((2^j)*t) = -1 mod n

    if (repeated_squaring_modular_exponentiation(b, t, n) == 1)
        return true;

    for (int j = 0; j < s; j++)
        if (repeated_squaring_modular_exponentiation(b, repeated_squaring_modular_exponentiation(2, j, n) * t, n) == n - 1)
            return true;

    return false;
}

int gcd(int a, int b) {
    return b == 0 ? a : gcd(b, a % b);
}

int main() {
    int p, q, n;
    cout << "Input p: ";
    cin >> p;
    cout << "Input q: ";
    cin >> q;
    n = p * q;
    int max_bases = (p - 1) * (q - 1);

    int temp = n - 1;
    int s = 0;
    while (!(temp % 2))
    {
        s++;
        temp /= 2;
    }

    cout << "n = " << n << '\n';
    cout << "The equation is:\n";
    cout << n << " - 1 = 2^" << s << " * " << temp << '\n';
    cout << "Number of bases to check: " << max_bases;

    cout << "\nn is pseudoprime with respect to the following bases: \n";
    int bases_checked = 0;
    int b = 1;
    string bases = "";
    int bases_found = 0;
    while (bases_checked < max_bases)
    {
        if (gcd(n, b) == 1)
        {
            if (miller_rabin(n, s, temp, b))
            {
                bases = bases + to_string(b) + ", ";
                bases_found++;
            }
            bases_checked++;
        }
        b++;
    }

    bases.resize(bases.size() - 2);
    cout << bases << "\n";
    cout << "Statistics: " << n << " is pseudoprime to " << bases_found * 100.0 / max_bases << "% of the bases.\n";
}