Problem 266: Pseudo Square Root

Project Euler Problem 266 Solution

EulerSolve provides an optimized solution for Project Euler Problem 266, Pseudo Square Root, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$P=\prod_{p<190} p.$$ We want the largest divisor $d\mid P$ such that $d\le \sqrt P$. This divisor is the pseudo square root of $P$. The implementation finally prints the value modulo $10^{16}$, but the search itself is done exactly. Mathematical Approach 1. Divisors as Subset Products Because $P$ is a product of distinct primes, every divisor corresponds to a subset of the prime list. If $S$ is the chosen subset, then $$d=\prod_{p\in S} p.$$ The condition $d\le \sqrt P$ is equivalent to $$d^2\le P.$$ So the problem is a subset-product maximization problem under an upper bound. 2. Logarithms Turn the Product Constraint into a Sum Constraint Multiplication is awkward to compare directly, but logarithms convert products into sums: $$\log d=\sum_{p\in S}\log p,\qquad \log d\le \tfrac12\log P.$$ That means we want a subset whose log-sum is as close as possible to $\tfrac12\log P$ from below. This is exactly the shape of a knapsack-style search. 3. Meet-in-the-Middle The prime list below 190 has 42 elements, so a brute-force scan over all $2^{42}$ subsets is too large. The code splits the primes into two halves, enumerates all subset sums in each half, sorts those lists, and then combines them efficiently....

Detailed mathematical approach

Problem Summary

The divisors of $12$ are $1,2,3,4,6$ and $12$, and the largest of them that does not exceed $\sqrt{12}\approx3.46$ is $3$. For a positive integer $n$, the largest divisor of $n$ that does not exceed $\sqrt n$ is called the pseudo square root $\operatorname{PSR}(n)$; the statement also gives $\operatorname{PSR}(3102)=47$. Let $p$ be the product of the $42$ primes below $190$, a number with $73$ digits. We must find $\operatorname{PSR}(p)\bmod10^{16}$. ^[1]

The whole solution in one sentence: split the $42$ primes into two halves, sort the products of each half by their logarithms, sweep through both lists to pair every product of the first half with its best partner below $\sqrt p$, and let exact integers decide among the closest pairs, because rounded logarithms cannot tell the closest divisors apart.

Mathematical Approach

1. Divisors are sets of primes

By the fundamental theorem of arithmetic, $3102=2\cdot3\cdot11\cdot47$ is a product of primes in only one way. ^[3] A divisor of $3102$ therefore contains each of these four primes at most once, and different sets of primes give different divisors. The $2^4=16$ sets give the divisors $1,2,3,6,11,22,33,47$ and $66,94,141,282,517,1034,1551,3102$. Since $\sqrt{3102}\approx55.70$, exactly the first eight do not exceed the square root, and the largest of them is $47$. In the same way, the divisors of $p$ correspond one to one to the subsets $S$ of the $42$ primes below $190$:

$$d=\prod_{q\in S}q,\qquad S\subseteq\{2,3,5,\dots,181\}.$$

So $p$ has $2^{42}=4\,398\,046\,511\,104$ divisors, and since all the numbers involved are integers, the condition $d\le\sqrt p$ can be tested exactly as $d^2\le p$.

2. Divisors come in pairs around the square root

If $d$ divides $p$, so does its partner $p/d$. Two numbers whose product is $p$ cannot both lie below $\sqrt p$ or both above it, and neither equals $\sqrt p$, since $p$ is not a perfect square: each prime occurs in it only once. The pairing therefore matches the divisors below $\sqrt p$ with those above it, so exactly $2^{41}$ lie below; and since $d\mapsto p/d$ reverses the order, the largest divisor below, $\operatorname{PSR}(p)$, is the partner of the smallest divisor above. For $3102$, the partner of $47$ is $66$. The OEIS lists the pseudo square roots of all $n$ as A033676 and these partners as A033677. ^[2] In logarithms, with $\tau=\tfrac12\ln p$, two partners lie equally far from $\tau$, one on each side:

$$\ln\frac{p}{d}-\tau=\tau-\ln d.$$

So no divisor of $p$ lies strictly between $\operatorname{PSR}(p)$ and its partner, a symmetry that Section 11 uses.

3. A balanced split of the primes

A divisor $d$ splits the primes into two groups, those in $d$ and those in its partner, and $\operatorname{PSR}(p)$ is the smaller product in the most balanced split. Since $\ln(ab)=\ln a+\ln b$, the task is to choose some of the numbers $\ln2,\ln3,\dots,\ln181$ with a sum as large as possible but at most

$$\tau=\tfrac12\ln p\approx83.736.$$

This is the optimisation form of the subset sum problem, a subset whose sum is at most a target $T$ and as close to $T$ as possible, ^[4] and of the partition problem; both are NP-hard in general. Simple rules fail: the greedy rule, which puts each number into the group with the smaller current sum, ^[5] splits the primes below $20$, taken from the largest down, into $2926\cdot3315$, worse than the split $3094\cdot3135$ of Section 6.

4. Meet in the middle

Trying all $2^{42}\approx4.4\cdot10^{12}$ subsets is far too slow. Instead we split the primes into a first half $2,3,\dots,73$ and a second half $79,83,\dots,181$ of $21$ primes each. Every divisor is $d=a\cdot b$ in exactly one way, where $a$ is a product of primes from the first half and $b$ a product of primes from the second half, and each half has $2^{21}=2\,097\,152$ such products. For a fixed $a$, the best partner is the largest $b$ with

$$b\le\frac{\sqrt p}{a}\iff a^2b^2\le p.$$

If the products $b$ are sorted, a binary search finds this partner in about $21$ steps, so $2^{21}$ searches replace $2^{42}$ trials. This is the meet-in-the-middle algorithm of Horowitz and Sahni for the subset sum problem. ^[4]

5. One pointer instead of many searches

With the $a$ sorted increasingly too, the searches merge into one sweep: as $a$ grows, the bound $\sqrt p/a$ shrinks, so a pointer into the sorted list of the $b$ only moves down. It starts at the largest $b$, and for each $a$ it moves down while $a\,b\gt\sqrt p$; if even $b=1$ is too large, this $a$ has no candidate. An entry that the pointer has passed failed for an earlier $a'\le a$, so $a\,b\ge a'b\gt\sqrt p$ and it fails for $a$ too; the pointer therefore stops at the best $b$. It moves at most $2^{21}$ times in total, so after sorting the sweep takes only about $2\cdot2^{21}$ steps.

For $3102$, the halves $\{2,3\}$ and $\{11,47\}$ give $a=1,2,3,6$ and $b=1,11,47,517$, and the bound is $55.70/a$. For $a=1$ the pointer moves from $517$ down to $47$, giving $47$. For $a=2$ the product $94$ is too large, so it moves to $11$, giving $22$. For $a=3$ it stays at $11$, giving $33$, and for $a=6$ it moves to $1$, giving $6$. The largest candidate is $47$.

6. A worked example: the primes below 20

For the eight primes below $20$ we have $p=9\,699\,690$ and $\sqrt p\approx3114.43$. The first half $2,3,5,7$ has the products

$1,2,3,5,6,7,10,14,15,21,$
$30,35,42,70,105,210,$

and the second half $11,13,17,19$ has the products

$1,11,13,17,19,143,187,209,$
$221,247,323,2431,$
$2717,3553,4199,46189.$

The sweep pairs $a=1$ with $b=2717$; $a=2,3,5,6,7$ with $323$; $a=10$ with $247$; $a=14$ with $221$; $a=15$ with $187$; $a=21$ with $143$; $a=30,35,42,70,105$ with $19$; and $a=210$ with $13$, while the pointer moves $13$ times. The three largest candidates are $3094=14\cdot221$, $3003=21\cdot143$ and $2805=15\cdot187$. So $\operatorname{PSR}(9\,699\,690)=3094$${}=2\cdot7\cdot13\cdot17$, whose partner is $3135=3\cdot5\cdot11\cdot19$, and a check of all $256$ divisors confirms it.

7. Logarithms as sort keys

For the real problem the products of the second half reach $79\cdot83\cdots181$${}\approx1.3\cdot10^{44}$, far beyond the $64$-bit integers. The program therefore stores, for every subset, its natural logarithm as a floating-point number of type double together with a bit mask of the chosen primes; the logarithm is cheap to build, sort and add, and the test $a\,b\le\sqrt p$ becomes $\ln a+\ln b\le\tau$. Sorting by logarithms is sorting by products, so the sweep of Section 5 works unchanged, as long as the logarithms are exact. They are not: a double has a significand of $53$ bits, so between $64$ and $128$ neighbouring doubles are $2^{-46}\approx1.4\cdot10^{-14}$ apart, and every computed value is rounded to one of them. ^[6]

8. How crowded the divisors are near the square root

Write $\ln d=\sum_i x_i\ln q_i$, where $q_1,\dots,q_{42}$ are the primes and $x_i$ is $1$ if $q_i$ divides $d$ and $0$ otherwise. If every $x_i$ were chosen by a fair coin toss, the term $x_i\ln q_i$ would have mean $\tfrac12\ln q_i$ and variance $\tfrac14(\ln q_i)^2$, so $\ln d$ would have mean $\tau$ and standard deviation

$$\sigma=\frac12\sqrt{\sum_{i=1}^{42}(\ln q_i)^2}\approx13.44.$$

The central limit theorem suggests approximating such a sum of many independent terms by a normal distribution, ^[7] whose density at its centre is $1/(\sigma\sqrt{2\pi})$. For intervals as narrow as ours this is only a heuristic, so we check it below. Each of the $2^{42}$ divisors is one outcome of the coin tosses, so near $\tau$ the heuristic predicts about

$$\frac{2^{42}}{\sigma\sqrt{2\pi}}\approx1.3\cdot10^{11}$$

divisors per unit of $\ln d$, one every $7.7\cdot10^{-12}$ on average. An exact count confirms this: $2534$ divisors of $p$ satisfy $|\ln d-\tau|\le10^{-8}$, where the estimate predicts about $2600$, and $1267$ of them lie on each side, as the pairing of Section 2 demands. The two closest of them differ by only $5.5\cdot10^{-15}$, less than the spacing of doubles there. Only one divisor on each side lies within $10^{-11}$ of $\tau$, namely $\operatorname{PSR}(p)$ and its partner, and the next divisor below $\operatorname{PSR}(p)$ is only about $2\cdot10^{-12}$ further down.

9. How large the rounding errors can be

The logarithm of each product is a sum of logarithms of primes, added one at a time, and the target is half the sum of all $42$ of them. A pair sum therefore uses at most $42$ logarithms and $41$ inexact additions: $20$ in each half and one to combine the halves, since adding to $0$ is exact. All these sums stay below $256$, where neighbouring doubles are at most $2^{-45}$ apart, so each addition, rounded to the nearest double, is off by at most $2^{-46}$. ^[6] The logarithm of each prime is below $8$, where doubles are at most $2^{-50}$ apart, and within one such unit of the exact value: Java guarantees this, ^[8] and on the test machine the library logarithm of the C++ and Python programs is even correctly rounded for all $42$ primes. Hence every computed pair sum and the computed target are within

$$\varepsilon=41\cdot2^{-46}+42\cdot2^{-50}\approx6.2\cdot10^{-13}$$

of the exact values. The measured errors are much smaller, about $10^{-14}$ for the target and $2\cdot10^{-15}$ for the sum of $\operatorname{PSR}(p)$, but the first already exceeds the smallest gap of Section 8. Logarithms can only find a short list of candidates.

10. Two passes: a filter and an exact decision

Write $\tilde s$ for the computed value of $\ln a+\ln b$ and $\tilde\tau$ for the computed target. The program sweeps twice.

Pass 1 runs the sweep of Section 5 with the test $\tilde s\le\tilde\tau+10^{-12}$ and records the largest accepted value $\beta$; it only tells pass 2 how high the best sums go.
Pass 2 sweeps again with the looser test $\tilde s\le\tilde\tau+10^{-10}$. For each $a$ it looks at the entry just above the pointer and at up to $13$ entries from the pointer down, moves on as soon as one of the latter has $\tilde s\lt\beta-10^{-9}$, and skips entries with $\tilde s\gt\tilde\tau+10^{-8}$. Every remaining pair is multiplied out exactly, and the largest $d=a\,b$ with $d^2\le p$ is kept.

Near the pointer, neighbouring entries of the second list are typically about $7\cdot10^{-5}$ apart, so for almost every $a$ the best sum falls short of $\beta$ by far more than $10^{-9}$, and pass 2 looks at two entries and moves on. It multiplies out only $1389$ candidates: $1249$ just above the pointer, all rejected because their squares exceed $p$, and $140$ from the pointer down.

11. Why the exact pass finds the pseudo square root

Let $s^{*}=\ln\operatorname{PSR}(p)$, write $\operatorname{PSR}(p)=a^{*}b^{*}$ with $a^{*}$ from the first half and $b^{*}$ from the second, and let $\tilde s^{*}$ be the computed sum of this pair. By Section 9, every computed sum and $\tilde\tau$ lie within $\varepsilon\lt10^{-12}$ of the exact values.

Step 1: pass 1 cannot overshoot by much. The value $\beta$ is the computed sum of some pair; let $s'$ be its exact sum. If $s'\le\tau$, the product of the pair is at most $\sqrt p$, hence at most $\operatorname{PSR}(p)$, and $\beta\le s'+\varepsilon\le s^{*}+\varepsilon$. If $s'\gt\tau$, the product exceeds $\sqrt p$, so by Section 2 it is at least the partner $p/\operatorname{PSR}(p)$, which means $s'-\tau\ge\tau-s^{*}$. Since $\beta\le\tilde\tau+10^{-12}$, this gives

$$\begin{aligned}\tau-s^{*}&\le s'-\tau\le\beta+\varepsilon-\tau\\&\le\tilde\tau+10^{-12}+\varepsilon-\tau\le10^{-12}+2\varepsilon,\end{aligned}$$

so the pseudo square root is then very close to the target, and

$$\begin{aligned}\beta&\le\tilde\tau+10^{-12}\le\tau+10^{-12}+\varepsilon\\&\le s^{*}+2\cdot10^{-12}+3\varepsilon.\end{aligned}$$

Since $s^{*}\le\tilde s^{*}+\varepsilon$, both cases together give

$$\begin{aligned}\beta&\le\tilde s^{*}+2\cdot10^{-12}+4\varepsilon\\&\lt\tilde s^{*}+10^{-11}.\end{aligned}$$

Step 2: pass 2 examines the right pair. When pass 2 handles $a^{*}$, the pointer stops at the last entry whose computed sum with $a^{*}$ is at most $\tilde\tau+10^{-10}$. Because $\tilde s^{*}\le s^{*}+\varepsilon$ and $s^{*}\lt\tau\le\tilde\tau+\varepsilon$, we have $\tilde s^{*}\lt\tilde\tau+2\varepsilon$, so $b^{*}$ passes this test. No entry above $b^{*}$ passes it. The exact logarithms of the second-half products are at least $1.77\cdot10^{-9}$ apart, far more than $2\varepsilon$, so the computed order is the exact order, and an entry $b\gt b^{*}$ gives the exact sum $u\ge s^{*}+1.77\cdot10^{-9}$. Its product exceeds $\operatorname{PSR}(p)$, hence $\sqrt p$, so Section 2 gives $u-\tau\ge\tau-s^{*}$; passing the test would mean $u\le\tau+1.02\cdot10^{-10}$, allowing for the rounding of the sum, the target and the tolerance, and then $u-s^{*}\le2.04\cdot10^{-10}$, a contradiction. So the pointer stops exactly on $b^{*}$, which the window examines at $t=0$. It is not skipped, since $\tilde s^{*}\lt\tilde\tau+10^{-8}$, and the stopping test does not fire, since $\tilde s^{*}\gt\beta-10^{-11}$; the additions of the tolerances are rounded too, by at most $2^{-46}$, which these margins easily absorb.

Step 3: the exact test decides. A candidate is kept only if its exact square is at most $p$, so the result is a divisor that does not exceed $\sqrt p$, hence at most $\operatorname{PSR}(p)$. The pseudo square root itself is examined by Step 2 and passes, so the result is exactly $\operatorname{PSR}(p)$.

The spacing $1.77\cdot10^{-9}$ is a property of these $21$ primes, checked with exact integer ratios, and an instrumented run confirms that the pointer stops on $b^{*}$. For other inputs the window of $13$ entries below the pointer is only a safety margin, not a proven bound.

The second case of Step 1 does occur. For the $35$ primes below $150$, the pseudo square root lies only about $6\cdot10^{-13}$ below the target, so $\beta$ belongs to its partner, whose square exceeds the product; pass 2 rejects the partner exactly and finds the pseudo square root about $1.2\cdot10^{-12}$ lower. The same happens for the $40$ primes below $175$.

How the Code Works

The C++ program ^[9] starts in main. parse_arguments accepts only --skip-checkpoints, without which run_checkpoints runs first. primes_below(190) lists the $42$ primes with the sieve of Eratosthenes, ^[10] solve_for_primes returns $\operatorname{PSR}(p)\lt\sqrt p\lt10^{37}$ as an exact integer, and mod_10_16 gives the printed remainder. The modulus is applied only at the end, because the search compares sizes and remainders do not preserve order.

solve_for_primes builds product_all, total_log and its half target_log, splits the primes at split $=21$ and calls enumerate_subsets on each half. A bit mask encodes a subset, bit $i$ standing for the $i$-th prime of the half, counted from $0$: for the primes $2,3,5,7$, the mask $10=1010_2$ selects $3$ and $7$. Each mask gets a Subset ($16$ bytes on the test machine) holding the logarithm as a double and the mask, at one addition per nonzero mask: mask & (~mask + 1U) isolates the lowest set bit, __builtin_ctz gives its position, and the entry of the smaller mask without that bit, already filled, plus that prime's logarithm gives the new entry. In the example, the mask $8$ selects only $7$, and adding $\ln3$ to its entry gives $\ln21$. std::sort orders the entries by logarithm. The sweeps of Section 10 use the pointer j, the window offset t from $-1$ to $12$, best_log for $\beta$ and the tolerances 1e-12, 1e-10, 1e-9 and 1e-8; build_from_masks multiplies a candidate out, square_leq compares its exact square with product_all, and best_value keeps the largest valid one.

BigInt has six $64$-bit words, ample for $p$ with its $242$ bits and for candidates whose squares exceed $p$ by a factor of at most about $1+2\cdot10^{-8}$; its schoolbook multiplication uses $128$-bit intermediate products, and div_mod divides by a machine word from the top word down. run_checkpoints compares the results for the $10$, $12$ and $15$ primes below $30$, $40$ and $50$, whose pseudo square roots are $79534$, $2714690$ and $783152070$, with brute_small, which tries every subset in unsigned $128$-bit arithmetic; on a mismatch it prints a message and main returns the exit code $2$.

Python and Java use the same sweeps, tolerances and window, without checkpoints or options. Python's enumerate_subsets tests all $21$ bits of every mask and adds the logarithm of each chosen prime, in a different order than C++, so some sums may differ from the C++ values in the last bits, within the bound of Section 9; its exact products are unbounded Python integers, and the residue is best_value % MOD. Java uses the same lowest-bit recurrence, with mask & (-mask) and Integer.numberOfTrailingZeros, sorts with Arrays.sort and uses BigInteger for the exact products and the final remainder.

Complexity and Verification

With $n=42$ primes and $H=2^{21}$ products per half, C++ and Java build each list with $H-1$ additions, while Python makes $21H$ bit tests and $21H/2$ additions per half, since each bit is set in half of the masks. Sorting costs about $H\log_2H=21H$ comparisons per half, and each sweep makes at most $H$ pointer moves; pass 2 looks at no more than $14$ window entries per product of the first half, and in fact at $4\,194\,444$ entries in all, about two per product. The time is therefore $O(2^{n/2}\,n)$, dominated by sorting in C++ and Java and by sorting and enumeration alike in Python, and the memory is $O(2^{n/2})$, compared with $2^{n}\approx4.4\cdot10^{12}$ subsets for brute force. On the test machine (an Apple M4 Max; g++ -O2, OpenJDK 17, Python 3.11) the C++ program runs in about $0.3$ seconds with a peak memory of about $70$ MB, mostly its two arrays of $2^{21}$ entries of $16$ bytes, $64$ MiB in all; Java needs about $0.7$ seconds, and Python about $6$ to $8$ seconds and $600$ MB, because each entry is a tuple of Python objects.

The following checks were made for this article and are not part of the published solutions. Trial division gives $\operatorname{PSR}(12)=3$ and $\operatorname{PSR}(3102)=47$. A Python copy of solve_for_primes agrees with a brute force over all subsets of the first $k$ primes for every $k\le24$, and for every $k\le41$ with an exact meet in the middle without logarithms, which finds for each $a$ by binary search the largest product $b$ with

$$b\le\Bigl\lfloor\sqrt{\lfloor p/a^2\rfloor}\Bigr\rfloor,$$

which for integers is equivalent to $a^2b^2\le p$. The C++ program run with the prime limits $100$, $150$, $160$, $175$ and $180$ agrees with the exact method, and for the primes below $190$ all three programs print the same residue as the exact method. The counts of Sections 8 to 11 come from the same exact lists and from an instrumented copy of the C++ program.

Footnotes and References

The references below supply the background facts used above; the pairing argument of Section 11, the error bound of Section 9 and all counts are derived or computed in this article.

Project Euler 266 — Pseudo Square Root. The official statement defines PSR(n) as the largest divisor of n that does not exceed √n, gives PSR(12) = 3 and PSR(3102) = 47, and asks for PSR(p) mod 10¹⁶ for the product p of the primes below 190.
OEIS A033676 — Largest divisor of n <= sqrt(n). The entry lists these values up to n = 10000 and gives a(n) = n/A033677(n), where A033677(n) is the smallest divisor of n that is at least √n.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 can be written as a product of primes in exactly one way, up to the order of the factors.
Wikipedia — Subset sum problem. The optimisation version asks for a subset whose sum is at most T and as close to T as possible; the problem is NP-hard. Horowitz and Sahni list the 2^(n/2) subset sums of each half in sorted order and pass through one list downwards and the other upwards, in time and memory O(2^(n/2)) when the lists are sorted by merging.
Wikipedia — Partition problem. Splitting numbers into two parts whose sums differ as little as possible is NP-hard, and the greedy algorithm, which puts each number into the part with the smaller current sum, is only an approximation.
Wikipedia — Double-precision floating-point format. A double has a 53-bit significand, and consecutive doubles between 2^k and 2^(k+1) are 2^(k−52) apart.
Wikipedia — Central limit theorem. Under Lyapunov's moment condition, standardised sums of independent, not necessarily identically distributed random variables converge in distribution to the standard normal distribution; counts in very narrow intervals need stronger local limit theorems.
Java SE 17 — Math.log. The computed natural logarithm must be within 1 ulp of the exact result and semi-monotonic.
C++ — Euler266.cpp. The linked immutable C++ revision contains enumerate_subsets, solve_for_primes, BigInt and run_checkpoints; it is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. The multiples of each prime may be crossed out from its square on, and the sieve may stop once the square of the current prime exceeds the limit.

Problem 266 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u128 = unsigned __int128;

struct Options {
  bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
  for (int i = 1; i < argc; ++i) {
    std::string arg(argv[i]);
    if (arg == "--skip-checkpoints") {
      options.run_checkpoints = false;
      continue;
    }
    std::cerr << "Unknown argument: " << arg << '\n';
    return false;
  }
  return true;
}

struct BigInt {
  static constexpr int WORDS = 6;
  u64 data[WORDS];

  BigInt() { std::memset(data, 0, sizeof(data)); }

  explicit BigInt(u64 val) {
    std::memset(data, 0, sizeof(data));
    data[0] = val;
  }

  bool operator<(const BigInt& other) const {
    for (int i = WORDS - 1; i >= 0; --i) {
      if (data[i] != other.data[i]) return data[i] < other.data[i];
    }
    return false;
  }

  bool operator<=(const BigInt& other) const { return !(other < *this); }

  bool operator>(const BigInt& other) const { return other < *this; }

  BigInt operator*(const BigInt& other) const {
    BigInt res;
    for (int i = 0; i < WORDS; ++i) {
      u128 carry = 0;
      for (int j = 0; j < WORDS - i; ++j) {
        const u128 prod = static_cast<u128>(data[i]) * other.data[j] + res.data[i + j] + carry;
        res.data[i + j] = static_cast<u64>(prod);
        carry = prod >> 64;
      }
    }
    return res;
  }

  void multiply_by(u64 p) {
    u128 carry = 0;
    for (int i = 0; i < WORDS; ++i) {
      const u128 prod = static_cast<u128>(data[i]) * p + carry;
      data[i] = static_cast<u64>(prod);
      carry = prod >> 64;
    }
  }

  u64 div_mod(u64 divisor) {
    u128 rem = 0;
    for (int i = WORDS - 1; i >= 0; --i) {
      const u128 cur = data[i] + (rem << 64);
      data[i] = static_cast<u64>(cur / divisor);
      rem = cur % divisor;
    }
    return static_cast<u64>(rem);
  }

  u64 mod_10_16() const {
    BigInt tmp = *this;
    return tmp.div_mod(10000000000000000ULL);
  }

  u128 low_u128() const {
    return static_cast<u128>(data[1]) << 64 | static_cast<u128>(data[0]);
  }
};

bool square_leq(const BigInt& x, const BigInt& bound) {
  const BigInt sq = x * x;
  return sq <= bound;
}

std::vector<u32> primes_below(int limit) {
  std::vector<unsigned char> is_prime(static_cast<std::size_t>(limit), 1);
  is_prime[0] = 0;
  is_prime[1] = 0;
  for (int p = 2; p * p < limit; ++p) {
    if (!is_prime[static_cast<std::size_t>(p)]) continue;
    for (int q = p * p; q < limit; q += p) is_prime[static_cast<std::size_t>(q)] = 0;
  }
  std::vector<u32> primes;
  for (int p = 2; p < limit; ++p) {
    if (is_prime[static_cast<std::size_t>(p)]) primes.push_back(static_cast<u32>(p));
  }
  return primes;
}

struct Subset {
  double log_value;
  u32 mask;

  bool operator<(const Subset& other) const { return log_value < other.log_value; }
};

std::vector<Subset> enumerate_subsets(const std::vector<u32>& primes) {
  const std::size_t n = primes.size();
  const std::size_t m = std::size_t{1} << n;
  std::vector<double> logs(n);
  for (std::size_t i = 0; i < n; ++i) logs[i] = std::log(static_cast<double>(primes[i]));

  std::vector<Subset> out(m);
  out[0] = {0.0, 0U};
  for (u32 mask = 1; mask < m; ++mask) {
    const u32 lsb = mask & (~mask + 1U);
    const int bit = __builtin_ctz(lsb);
    const u32 prev = mask ^ lsb;
    out[mask].log_value = out[prev].log_value + logs[static_cast<std::size_t>(bit)];
    out[mask].mask = mask;
  }

  std::sort(out.begin(), out.end());
  return out;
}

BigInt build_from_masks(const std::vector<u32>& left_primes, const std::vector<u32>& right_primes,
                        u32 left_mask, u32 right_mask) {
  BigInt v(1);
  u32 lm = left_mask;
  while (lm) {
    const int b = __builtin_ctz(lm);
    v.multiply_by(static_cast<u64>(left_primes[static_cast<std::size_t>(b)]));
    lm &= (lm - 1U);
  }
  u32 rm = right_mask;
  while (rm) {
    const int b = __builtin_ctz(rm);
    v.multiply_by(static_cast<u64>(right_primes[static_cast<std::size_t>(b)]));
    rm &= (rm - 1U);
  }
  return v;
}

BigInt solve_for_primes(const std::vector<u32>& primes) {
  BigInt product_all(1);
  double total_log = 0.0;
  for (u32 p : primes) {
    product_all.multiply_by(static_cast<u64>(p));
    total_log += std::log(static_cast<double>(p));
  }
  const double target_log = total_log * 0.5;

  const std::size_t split = primes.size() / 2;
  std::vector<u32> left_primes(primes.begin(), primes.begin() + static_cast<std::ptrdiff_t>(split));
  std::vector<u32> right_primes(primes.begin() + static_cast<std::ptrdiff_t>(split), primes.end());

  const std::vector<Subset> left = enumerate_subsets(left_primes);
  const std::vector<Subset> right = enumerate_subsets(right_primes);

  double best_log = -1.0;
  std::int64_t j = static_cast<std::int64_t>(right.size()) - 1;
  for (const Subset& l : left) {
    while (j > 0 && l.log_value + right[static_cast<std::size_t>(j)].log_value > target_log + 1e-12) --j;
    const double cand = l.log_value + right[static_cast<std::size_t>(j)].log_value;
    if (cand <= target_log + 1e-12 && cand > best_log) best_log = cand;
  }

  BigInt best_value(0);
  j = static_cast<std::int64_t>(right.size()) - 1;
  for (const Subset& l : left) {
    while (j > 0 && l.log_value + right[static_cast<std::size_t>(j)].log_value > target_log + 1e-10) --j;

    for (int t = -1; t <= 12; ++t) {
      const std::int64_t k = j - t;
      if (k < 0 || k >= static_cast<std::int64_t>(right.size())) continue;

      const Subset& r = right[static_cast<std::size_t>(k)];
      const double cand_log = l.log_value + r.log_value;
      if (t >= 0 && cand_log + 1e-9 < best_log) break;
      if (cand_log > target_log + 1e-8) continue;

      const BigInt cand = build_from_masks(left_primes, right_primes, l.mask, r.mask);
      if (square_leq(cand, product_all) && cand > best_value) best_value = cand;
    }
  }

  return best_value;
}

u128 brute_small(const std::vector<u32>& primes) {
  const int n = static_cast<int>(primes.size());
  const u32 all = (n == 32 ? 0xFFFFFFFFU : ((1U << n) - 1U));
  (void)all;

  u128 P = 1;
  for (u32 p : primes) P *= static_cast<u128>(p);

  const u64 lim = (n >= 31) ? 0 : (1ULL << n);
  u128 best = 0;
  for (u64 mask = 0; mask < lim; ++mask) {
    u128 v = 1;
    for (int i = 0; i < n; ++i) {
      if (mask & (1ULL << i)) v *= static_cast<u128>(primes[static_cast<std::size_t>(i)]);
    }
    if (v * v <= P && v > best) best = v;
  }
  return best;
}

bool run_checkpoints() {
  {
    std::vector<u32> p = primes_below(30);
    const BigInt got = solve_for_primes(p);
    if (got.low_u128() != brute_small(p)) {
      std::cerr << "Checkpoint failed for primes<30" << '\n';
      return false;
    }
  }
  {
    std::vector<u32> p = primes_below(40);
    const BigInt got = solve_for_primes(p);
    if (got.low_u128() != brute_small(p)) {
      std::cerr << "Checkpoint failed for primes<40" << '\n';
      return false;
    }
  }
  {
    std::vector<u32> p = primes_below(50);
    const BigInt got = solve_for_primes(p);
    if (got.low_u128() != brute_small(p)) {
      std::cerr << "Checkpoint failed for primes<50" << '\n';
      return false;
    }
  }
  return true;
}

}  // namespace

int main(int argc, char** argv) {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);

  Options options;
  if (!parse_arguments(argc, argv, options)) return 1;
  if (options.run_checkpoints && !run_checkpoints()) return 2;

  const std::vector<u32> primes = primes_below(190);
  const BigInt ans = solve_for_primes(primes);
  std::cout << ans.mod_10_16() << '\n';
  return 0;
}

Python

import math

def solve():
    def primes_below(limit):
        is_p = bytearray(b'\x01' * limit)
        is_p[0] = 0
        is_p[1] = 0
        for p in range(2, limit):
            if p * p >= limit:
                break
            if is_p[p]:
                for q in range(p*p, limit, p):
                    is_p[q] = 0
        return [p for p in range(2, limit) if is_p[p]]

    primes = primes_below(190)
    n = len(primes)
    product_all = 1
    for p in primes:
        product_all *= p
    total_log = sum(math.log(p) for p in primes)
    target_log = total_log * 0.5

    split = n // 2
    left_primes = primes[:split]
    right_primes = primes[split:]

    def enumerate_subsets(ps):
        m = len(ps)
        logs = [math.log(p) for p in ps]
        subsets = []
        for mask in range(1 << m):
            log_val = 0.0
            for i in range(m):
                if mask & (1 << i):
                    log_val += logs[i]
            subsets.append((log_val, mask))
        subsets.sort()
        return subsets

    left = enumerate_subsets(left_primes)
    right = enumerate_subsets(right_primes)

    def build_value(lmask, rmask):
        v = 1
        for i in range(len(left_primes)):
            if lmask & (1 << i):
                v *= left_primes[i]
        for i in range(len(right_primes)):
            if rmask & (1 << i):
                v *= right_primes[i]
        return v

    # Two-pointer pass to find best log
    best_log = -1.0
    j = len(right) - 1
    for lv, lm in left:
        while j > 0 and lv + right[j][0] > target_log + 1e-12:
            j -= 1
        cand = lv + right[j][0]
        if cand <= target_log + 1e-12 and cand > best_log:
            best_log = cand

    # Second pass to find exact best value
    best_value = 0
    j = len(right) - 1
    for lv, lm in left:
        while j > 0 and lv + right[j][0] > target_log + 1e-10:
            j -= 1
        for t in range(-1, 13):
            k = j - t
            if k < 0 or k >= len(right):
                continue
            rv, rm = right[k]
            cand_log = lv + rv
            if t >= 0 and cand_log + 1e-9 < best_log:
                break
            if cand_log > target_log + 1e-8:
                continue
            cand = build_value(lm, rm)
            if cand * cand <= product_all and cand > best_value:
                best_value = cand

    MOD = 10**16
    return str(best_value % MOD)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;

public class Euler266 {

    static List<Integer> primesBelow(int limit) {
        boolean[] isPrime = new boolean[limit];
        Arrays.fill(isPrime, true);
        isPrime[0] = false;
        isPrime[1] = false;
        for (int p = 2; p * p < limit; p++) {
            if (isPrime[p]) {
                for (int q = p * p; q < limit; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int p = 2; p < limit; p++) {
            if (isPrime[p])
                primes.add(p);
        }
        return primes;
    }

    static class Subset implements Comparable<Subset> {
        double logValue;
        int mask;

        Subset(double logValue, int mask) {
            this.logValue = logValue;
            this.mask = mask;
        }

        @Override
        public int compareTo(Subset other) {
            return Double.compare(this.logValue, other.logValue);
        }
    }

    static Subset[] enumerateSubsets(List<Integer> primes) {
        int n = primes.size();
        int m = 1 << n;
        double[] logs = new double[n];
        for (int i = 0; i < n; i++) {
            logs[i] = Math.log(primes.get(i));
        }

        Subset[] out = new Subset[m];
        out[0] = new Subset(0.0, 0);
        for (int mask = 1; mask < m; mask++) {
            int lsb = mask & (-mask);
            int bit = Integer.numberOfTrailingZeros(lsb);
            int prev = mask ^ lsb;
            out[mask] = new Subset(out[prev].logValue + logs[bit], mask);
        }
        Arrays.sort(out);
        return out;
    }

    static BigInteger buildFromMasks(List<Integer> leftPrimes, List<Integer> rightPrimes, int leftMask, int rightMask) {
        BigInteger v = BigInteger.ONE;
        int lm = leftMask;
        while (lm != 0) {
            int b = Integer.numberOfTrailingZeros(lm);
            v = v.multiply(BigInteger.valueOf(leftPrimes.get(b)));
            lm &= (lm - 1);
        }
        int rm = rightMask;
        while (rm != 0) {
            int b = Integer.numberOfTrailingZeros(rm);
            v = v.multiply(BigInteger.valueOf(rightPrimes.get(b)));
            rm &= (rm - 1);
        }
        return v;
    }

    static BigInteger solveForPrimes(List<Integer> primes) {
        BigInteger productAll = BigInteger.ONE;
        double totalLog = 0.0;
        for (int p : primes) {
            productAll = productAll.multiply(BigInteger.valueOf(p));
            totalLog += Math.log(p);
        }
        double targetLog = totalLog * 0.5;

        int split = primes.size() / 2;
        List<Integer> leftPrimes = new ArrayList<>(primes.subList(0, split));
        List<Integer> rightPrimes = new ArrayList<>(primes.subList(split, primes.size()));

        Subset[] left = enumerateSubsets(leftPrimes);
        Subset[] right = enumerateSubsets(rightPrimes);

        double bestLog = -1.0;
        int j = right.length - 1;
        for (Subset l : left) {
            while (j > 0 && l.logValue + right[j].logValue > targetLog + 1e-12)
                j--;
            double cand = l.logValue + right[j].logValue;
            if (cand <= targetLog + 1e-12 && cand > bestLog) {
                bestLog = cand;
            }
        }

        BigInteger bestValue = BigInteger.ZERO;
        j = right.length - 1;
        for (Subset l : left) {
            while (j > 0 && l.logValue + right[j].logValue > targetLog + 1e-10)
                j--;

            for (int t = -1; t <= 12; t++) {
                int k = j - t;
                if (k < 0 || k >= right.length)
                    continue;

                Subset r = right[k];
                double candLog = l.logValue + r.logValue;
                if (t >= 0 && candLog + 1e-9 < bestLog)
                    break;
                if (candLog > targetLog + 1e-8)
                    continue;

                BigInteger cand = buildFromMasks(leftPrimes, rightPrimes, l.mask, r.mask);
                if (cand.multiply(cand).compareTo(productAll) <= 0 && cand.compareTo(bestValue) > 0) {
                    bestValue = cand;
                }
            }
        }
        return bestValue;
    }

    public static void main(String[] args) {
        List<Integer> primes = primesBelow(190);
        BigInteger ans = solveForPrimes(primes);
        System.out.println(ans.remainder(new BigInteger("10000000000000000")).toString());
    }
}