My code

… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <vector>
       
#include <algorithm>
       
#include <cmath>
       
 
       
// all primes below 190

       
std::vector<unsigned int> 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<double> 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<LogAndBitmask> 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;
       
}

This solution contains 16 empty lines, 35 comments and 4 preprocessor commands.