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

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.

For a left subset $a$ and a right subset $b$, the candidate is valid when

$$\log a+\log b\le \tfrac12\log P.$$

Among all valid pairs, we need the one with the greatest total logarithm.

4. Exact Validation After Log Filtering

Floating-point logs are used only to narrow the search. They do not decide the final answer. Once the algorithm finds a near-optimal pair of masks, it reconstructs the exact product with big integers and checks

$$d^2\le P$$

using exact arithmetic. This avoids rounding errors and guarantees the final divisor is truly maximal.

5. Why the Search Is Efficient

Each half has roughly $2^{21}$ subsets. Enumerating both halves is still feasible, especially because subset log sums can be built incrementally from bit masks:

if a mask removes its lowest set bit, the new subset sum is the previous sum plus the log of the corresponding prime.

After sorting, the code uses a monotone scan over the two lists, then inspects a small neighborhood of right-hand candidates around the current boundary to catch the best exact solution.

6. Small Toy Example

Suppose the primes were just $\{2,3,5,7\}$, split as $\{2,3\}$ and $\{5,7\}$. The left subsets give log-values for $1,2,3,6$; the right subsets give log-values for $1,5,7,35$.

We then look for the best pair whose combined log does not exceed half of $\log(2\cdot3\cdot5\cdot7)$. The real problem behaves the same way, just on a much larger prime set.

In this toy instance

$$P=2\cdot3\cdot5\cdot7=210,\qquad \sqrt P\approx 14.49.$$

The best valid divisor is therefore

$$14=2\cdot7,$$

whereas the nearby product $15=3\cdot5$ already fails because $15^2=225>210$. This makes the optimization target concrete: we are not looking for the most balanced subset, but for the largest subset product that stays just below the square-root barrier.

7. Modular Output

The final pseudo square root is enormous, so the program prints

$$d \bmod 10^{16}.$$

The value itself is computed exactly in the internal big integer representation, and only the displayed form is reduced modulo $10^{16}$.

How the Code Works

The program only accepts --skip-checkpoints. The function primes_below(190) generates the 42 primes used in the problem. enumerate_subsets computes every subset mask together with its log-sum and sorts the results. The key implementation trick is incremental subset construction: for a nonzero mask, the code removes its lowest set bit, reuses the previous mask sum, and adds one prime log via __builtin_ctz.

solve_for_primes computes the total log target, splits the primes into two halves, and searches the best valid pair by combining the sorted lists. A monotone pointer finds the boundary in the right list, and then only a very small neighborhood around that boundary is rechecked with exact arithmetic. So logarithms guide the search, but exact integer comparison still decides the winner.

Once a promising pair of masks is found, build_from_masks reconstructs the exact product as a BigInt. The helper square_leq checks whether the square stays below the full product. The BigInt type implements fixed-width multiword multiplication and exact division by $10^{16}$ for the final printout.

The checkpoint routine compares the fast solver with a brute-force search for smaller prime limits such as 30, 40, and 50. These checks confirm that the meet-in-the-middle logic and the exact validation agree with exhaustive search on small instances.

Complexity Analysis

If there are $n$ primes, each half has about $2^{n/2}$ subsets. Enumeration, sorting, and the guided neighborhood scan dominate the cost, so the runtime is about

$$O(2^{n/2}\log 2^{n/2})$$

and the memory usage is about

$$O(2^{n/2}).$$

For $n=42$, this is feasible, while the naive $2^{42}$ search is not.

Footnotes and References

Problem page: https://projecteuler.net/problem=266
Meet-in-the-middle technique: Wikipedia - Meet-in-the-middle attack
Subset sum problem: Wikipedia - Subset sum problem
Logarithms and monotone search: Wikipedia - Logarithm
Big integer arithmetic: Wikipedia - Arbitrary-precision arithmetic

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());
    }
}