#include <iostream>
#include <vector>
#include <chrono>
#include "mpirxx.h"       

using namespace std;

int primeLimit = 10000000;

mpz_class powm(mpz_class base, mpz_class exp, mpz_class mod)
{
	mpz_class result = 0;
	mpz_powm (result.get_mpz_t(), base.get_mpz_t(), exp.get_mpz_t(), mod.get_mpz_t());
	return result;
}

mpz_class gcd(mpz_class bu1, mpz_class bu2)
{
	mpz_class result = 0;
	mpz_gcd (result.get_mpz_t(), bu1.get_mpz_t(), bu2.get_mpz_t());
	return result;
}


int main()
{	
	cout << "Carmichael Numbers:" << endl << endl;

	////////////////////////// Generate primes lookup table ///////////////////////////////////////////////////////////

	vector<bool> primes(primeLimit + 1, true); // calculated primes
	primes.at(0) = primes.at(1) = false;
	int iMax = (int)sqrt((double)primeLimit) + 1; 
	for (int i = 2; i < iMax; ++i)
		if (primes.at(i))
			for (int j = 2 * i; j < primeLimit + 1; j += i)
				primes.at(j) = false; // Erastothenes sieve

	//////////////////////////  Fermat's little theorem ///////////////////////////////////////////////////////////////
	
	auto start_time = chrono::high_resolution_clock::now();

	for (int p = 2; p <= primeLimit; ++p)
	{
		uint64_t contra = 0;

		for (mpz_class a = 2; a <= 19; ++a)
		{
			mpz_class d = gcd(a, p);
			if (d != 1 || ((primes.at(a.get_ui())) == false))
				continue; 
			if (powm(a, p - 1, p) != 1)
				++contra;
		}//for a

		if (contra == 0)
		{
			if (primes.at(p) == false)
			{	
				auto end_time = chrono::high_resolution_clock::now();
				cout << p;								
				cout << "\ttime: " << chrono::duration_cast<chrono::milliseconds>(end_time - start_time).count() << " ms" << endl;
			}			
		}
	}//for p
}