ProjectEuler48.cpp

#include "stdafx.h"
#include "boost/multiprecision/cpp_int.hpp"
#include <vector>
#include <ppl.h>
#include <numeric>
#include <chrono>
#include <type_traits>

template <typename Integral>
Integral modularProduct(Integral a, Integral b, Integral c) {
	return ((a%c) * (b%c)) % c;
}

template <typename Integral>
Integral modularSum(Integral a, Integral b, Integral c) {
	return ((a%c) + (b%c)) % c;
}

template <typename Integral>
Integral modularPow(Integral a, Integral b, Integral c) {
	if (b == 0) {
		return Integral(1);
	}
	else if (b == 1) {
		return a % c;
	}
	Integral const half = b / Integral(2);
	const auto t = modularPow(a, half, c);
	const auto tSquared = modularProduct(t, t, c);
	if (b % 2 == 0) {
		return tSquared;
	}
	else {
		return (tSquared * modularPow(a, Integral(1), c)) % c;
	}
}

int main()
{
	using namespace boost::multiprecision;
	auto const start = std::chrono::system_clock::now();

	cpp_int modulo = 10'000'000'000;

	//************************* Parallel Version ************************************
	std::vector<cpp_int> range = [] {
		std::vector<cpp_int> v(1'000'000);
		std::iota(std::begin(v), std::end(v), 1);
		return v;
	}();

	concurrency::parallel_transform(
		std::begin(range), std::end(range), std::begin(range),
		[&](auto const& v) { return modularPow(v, v, modulo); });

	auto const product = concurrency::parallel_reduce(
		std::begin(range), std::end(range), cpp_int(0), 
		[&](auto const& x, auto const& y) { return modularSum(x, y, modulo); });

	//************************ Sequential Version ************************************
	//cpp_int product = 0;
	//for (cpp_int i = 1; i <= 1'000'000; ++i) {
	//	product = modularSum(modularPow(i, i, modulo), product, modulo);
	//}

	auto const end = std::chrono::system_clock::now();

	std::cout
		<< product
		<< " in "
		<< std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
		<< " ms"
		<< std::endl;

    return 0;
}