// ////////////////////////////////////////////////////////
// # Title
// Pseudo Square Root
//
// # URL
// https://projecteuler.net/problem=266
// http://euler.stephan-brumme.com/266/
//
// # Problem
// The divisors of 12 are: 1,2,3,4,6 and 12.
// The largest divisor of 12 that does not exceed the square root of 12 is 3.
// We shall call the largest divisor of an integer n that does not exceed the square root of n the pseudo square root (PSR) of n.
// It can be seen that PSR(3102)=47.
//
// Let p be the product of the primes below 190.
// Find PSR(p) mod `10^16`.
//
// # Solved by
// Stephan Brumme
// July 2017
//
// # Algorithm
// The product of all primes below 190 is `2 * 3 * 5 * 7 * ... * 181 = 5397346292805549782720214077673687806275517530364350655459511599582614290`.
// That number is obviously too large for my poor C++ compiler ... but its logarithm isn't:
// `log(2) + log(3) + log(5) + ... + log(181) = log(53973462...2614290) approx 167.472`
//
// My program looks for an optimal subsets of those primes where
// `log(p_1) + log(p_2) + ... + log(p_n) approx 167.472 / 2 approx 83.736 = ` (because `log(sqrt{x}) = log(x) / 2`)
//
// There are `2^42 approx 4.4 * 10^12` subsets. And again, I have a number that is too large for my computer ...
// But there is hope - I split those 42 primes into two equally-sized sets:
// the container ''left'' will be "responsible" for the first 21 primes, and ''right'' will control the next 21 primes.
//
// Now there are `2^21 approx 2.1 million` subsets of each ''left'' and ''right''.
// For each of those subsets I compute the logarithm of the product of its primes and a bitmask that indicates which primes were used.
// If the n-th bit is set the n-th prime was used (that's the case for ''left'') or the (n+21)-th prime was used (''right'').
//
// For each value of ''left'' I need to find the largest value of ''right'' such that their sum doesn't exceed 83.736 (''logRoot'').
// ''std::upper_bound'' can perform this task very fast if the input data is sorted (but it returns an iterator that is "one step too far").
//
// Until now, all computations were based on logarithms. I could try `e^{best sum}` but rounding issues produce only garbage.
// That's why I have to look at the bitmasks of the best combination of values from ''left'' and ''right''.
//
// I multiply all used primes according to ''leftMask'' and ''rightMask'' and repeatedly apply modulo `10^16`.
//
// # Note
// Programming languages with proper support for large number probably don't need to deal with logarithms.
#include
#include
#include
#include
// all primes below 190
std::vector primes =
{
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181
};
// store a subsets logarithm and a bitmask of used primes
struct LogAndBitmask
{
// logarithm of product of used primes
double logarithm;
// bitmask of used primes
unsigned int bitmask;
// for std::sort and std::upper_bound
bool operator<(const LogAndBitmask& other) const
{
return logarithm < other.logarithm;
}
};
int main()
{
unsigned int maxPrime = 190;
std::cin >> maxPrime;
// keep only relevant primes
while (primes.back() > maxPrime)
primes.pop_back();
// compute logarithm of all primes
double logProduct = 0;
std::vector logPrimes;
for (auto p : primes)
{
auto current = log(p);
logPrimes.push_back(current);
logProduct += current;
}
// log(sqrt(x)) = log(x) / 2
double logRoot = logProduct / 2;
// subdivides primes into two equally-sized sets (21 elements each)
const auto half = primes.size() / 2;
// generate all subsets of the lower half of the primes (left)
// and all subsets of the upper half of the primes (right)
// key is the logarithm of product of the used primes, value is the bitmask
// (if n-th bit is set, then the n-th prime is used)
std::vector right;
for (unsigned int bitmask = 0; bitmask < (1ULL << half); bitmask++)
{
// find logarithm for the upper half of primes
double logRight = 0;
for (unsigned int pos = 0; pos < half; pos++)
if ((bitmask & (1 << pos)) != 0)
logRight += logPrimes[pos + half];
// store logarithms and their generating bitmasks
if (logRight <= logRoot) // skip a few of the very large products
right.push_back({ logRight, bitmask});
}
// sort container / required for fast binary search (the algorithm behind std::upper_bound)
std::sort(right.begin(), right.end());
// find "best" combination of left and right:
// for each value of "left" look for best match in "right"
// keep the highest sum which doesn't exceed logRoot
double best = 0;
unsigned int leftMask = 0; // optimal bitmask for lower primes
unsigned int rightMask = 0; // optimal bitmask for larger primes
for (unsigned int bitmask = 0; bitmask < (1ULL << half); bitmask++)
{
// find logarithm for the same bitmask of the lower half of primes
double logLeft = 0;
for (unsigned int pos = 0; pos < half; pos++)
if ((bitmask & (1 << pos)) != 0)
logLeft += logPrimes[pos];
// how much should come from the bigger primes ?
LogAndBitmask missing = { logRoot - logLeft, 0 }; // bitmask doesn't matter, set to zero
// find best match
auto iteLogRight = std::upper_bound(right.begin(), right.end(), missing);
// ... which will slightly overshoot, therefore go back one step
iteLogRight--;
auto logRight = iteLogRight->logarithm;
// higher than before ?
if (best < logLeft + logRight)
{
best = logLeft + logRight;
// store bitmasks
leftMask = bitmask;
rightMask = iteLogRight->bitmask;
}
}
// last 16 digits of the product
unsigned long long result = 1;
// multiply all relevant primes (according to leftMask and rightMask)
const auto Modulo = 10000000000000000ULL;
// index of current prime
unsigned int currentPrime = 0;
// the first 21 primes ...
for (; currentPrime < half; currentPrime++)
{
// use that prime ?
if ((leftMask & 1) != 0)
{
result *= primes[currentPrime];
result %= Modulo;
}
// next bit
leftMask >>= 1;
}
// ... and the same for the next 21 primes
for (; currentPrime < primes.size(); currentPrime++)
{
// use that prime ?
if ((rightMask & 1) != 0)
{
result *= primes[currentPrime];
result %= Modulo;
}
// next bit
rightMask >>= 1;
}
// print result
std::cout << result << std::endl;
return 0;
}