Project Euler

Pseudo Square Root

The largest divisor of n that does not exceed sqrt(n) is called the pseudo square root (PSR) of n. Let P = prod_(p < 190, p prime) p be the product of the 42 primes below 190. Find the last 16 digi...

Source sync Apr 19, 2026

Problem #0266

Level Level 10

Solved By 1,981

Languages C++, Python

Answer 1096883702440585

Length 355 words

modular_arithmeticnumber_theorysearch

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 ($\operatorname {PSR}$) of $n$.

It can be seen that $\operatorname {PSR}(3102)=47$.

Let $p$ be the product of the primes below $190$.

Find $\operatorname {PSR}(p) \bmod 10^{16}$.

Problem 266: Pseudo Square Root

Mathematical Foundation

Theorem 1 (Divisor structure of a primorial). Since $P$ is squarefree (a product of distinct primes), every divisor of $P$ has the form $d = \prod_{p \in S} p$ for some subset $S \subseteq \mathcal{P}$ , where $\mathcal{P}$ is the set of primes below 190. The total number of divisors is $2^{42}$ .

Proof. Since $P = \prod_{p \in \mathcal{P}} p$ with all primes distinct, the exponent of each prime in any divisor is either 0 or 1. The divisors are thus in bijection with subsets of $\mathcal{P}$ , of which there are $2^{|\mathcal{P}|} = 2^{42}$ . $\square$

Theorem 2 (Logarithmic reformulation). The PSR problem is equivalent to finding a subset $S \subseteq \mathcal{P}$ maximizing $\sum_{p \in S} \ln p$ subject to $\sum_{p \in S} \ln p \le \frac{1}{2} \sum_{p \in \mathcal{P}} \ln p$ . This is a variant of the 0-1 knapsack problem.

Proof. We seek the largest divisor $d \le \sqrt{P}$ . Taking logarithms, $\ln d = \sum_{p \in S} \ln p$ and the constraint $d \le \sqrt{P}$ becomes $\ln d \le \frac{1}{2} \ln P = \frac{1}{2} \sum_{p \in \mathcal{P}} \ln p$ . Maximizing $d$ is equivalent to maximizing $\ln d$ . $\square$

Theorem 3 (Meet-in-the-middle correctness). Partition $\mathcal{P} = A \sqcup B$ with $|A| = |B| = 21$ . Let $\{(L_i, v_i)\}$ be the list of all subset products from $A$ (with $L_i = \log$ -value and $v_i =$ value $\bmod 10^{16}$ ), sorted by $L_i$ . For each subset product $(L_j, v_j)$ from $B$ , define $L^* = \tfrac{1}{2}\ln P - L_j$ and find the largest $L_i \le L^*$ by binary search. Then the overall maximum of $L_i + L_j$ (subject to $\le \tfrac{1}{2}\ln P$ ) is the logarithm of $\mathrm{PSR}(P)$ , and the corresponding product $v_i \cdot v_j \bmod 10^{16}$ gives the last 16 digits.

Proof. Every divisor $d$ of $P$ factors uniquely as $d = d_A \cdot d_B$ where $d_A$ divides $\prod_{p \in A} p$ and $d_B$ divides $\prod_{p \in B} p$ . The decomposition is exhaustive: all $2^{42}$ divisors are covered. For a fixed $d_B$ , the optimal $d_A$ is the largest one with $\ln d_A \le \frac{1}{2}\ln P - \ln d_B$ , found by binary search on the sorted list. Since both lists have size $2^{21}$ and are fully enumerated, the global optimum is found. $\square$

Editorial

Give the last 16 digits of d. Approach: Meet-in-the-middle on the 42 primes. Split into two halves, enumerate subset products for each, then for each product in the second half, binary search the first half for the largest complementary product that keeps d <= sqrt(P). We enumerate group A. We then iterate over each subset S of A. Finally, iterate over each subset of B, binary search in listA.

Pseudocode

enumerate group A
for each subset S of A
for each subset of B, binary search in listA
for each subset S of B

Complexity Analysis

Time: $O(2^{n/2} \cdot n)$ where $n = 42$ . Enumerating each half takes $O(2^{21} \cdot 21) \approx 4.4 \times 10^7$ operations. Binary search over $2^{21}$ elements costs $O(21)$ per query, with $2^{21}$ queries: total $O(2^{21} \cdot 21) \approx 4.4 \times 10^7$ .
Space: $O(2^{n/2}) = O(2^{21}) \approx 2 \times 10^6$ entries for the sorted list.

Answer

\boxed{1096883702440585}

C++ project_euler/problem_266/solution.cpp

#include <bits/stdc++.h>
using namespace std;

int main() {
    // Find all primes below 190
    vector<int> primes;
    for (int i = 2; i < 190; i++) {
        bool is_prime = true;
        for (int j = 2; j * j <= i; j++) {
            if (i % j == 0) { is_prime = false; break; }
        }
        if (is_prime) primes.push_back(i);
    }

    int n = primes.size(); // 42
    int half = n / 2;      // 21
    int other = n - half;   // 21

    const long long MOD = 10000000000000000LL; // 10^16

    // Compute log(sqrt(P)) using long double for precision
    long double log_sqrtP = 0;
    for (int p : primes) log_sqrtP += logl((long double)p);
    log_sqrtP /= 2.0L;

    // Modular multiplication using __int128
    auto mulmod = [&](long long a, long long b) -> long long {
        return (long long)((__int128)a * b % MOD);
    };

    // Generate all subset products for first half
    int sz1 = 1 << half;
    vector<long double> logA(sz1);
    vector<long long> modA(sz1);

    for (int mask = 0; mask < sz1; mask++) {
        long double lg = 0;
        long long val = 1;
        for (int i = 0; i < half; i++) {
            if (mask & (1 << i)) {
                lg += logl((long double)primes[i]);
                val = mulmod(val, primes[i]);
            }
        }
        logA[mask] = lg;
        modA[mask] = val;
    }

    // Sort by log value
    vector<int> idxA(sz1);
    iota(idxA.begin(), idxA.end(), 0);
    sort(idxA.begin(), idxA.end(), [&](int a, int b) { return logA[a] < logA[b]; });

    vector<long double> sortedLogA(sz1);
    vector<long long> sortedModA(sz1);
    for (int i = 0; i < sz1; i++) {
        sortedLogA[i] = logA[idxA[i]];
        sortedModA[i] = modA[idxA[i]];
    }

    // For each subset of second half, find best match
    int sz2 = 1 << other;
    long double best_log = -1;
    long long best_mod = 0;

    for (int mask = 0; mask < sz2; mask++) {
        long double lg = 0;
        long long val = 1;
        for (int i = 0; i < other; i++) {
            if (mask & (1 << i)) {
                lg += logl((long double)primes[half + i]);
                val = mulmod(val, primes[half + i]);
            }
        }

        long double target = log_sqrtP - lg;
        if (target < -1e-15) continue;

        // Binary search: find largest index where sortedLogA[idx] <= target
        int lo = 0, hi = sz1 - 1, pos = -1;
        while (lo <= hi) {
            int mid = (lo + hi) / 2;
            if (sortedLogA[mid] <= target + 1e-15) {
                pos = mid;
                lo = mid + 1;
            } else {
                hi = mid - 1;
            }
        }

        if (pos < 0) continue;

        // Check a few candidates near pos for floating point safety
        for (int j = max(0, pos - 3); j <= min(sz1 - 1, pos + 3); j++) {
            if (sortedLogA[j] > target + 1e-15) continue;
            long double total_log = sortedLogA[j] + lg;
            if (total_log > best_log + 1e-15) {
                best_log = total_log;
                best_mod = mulmod(sortedModA[j], val);
            } else if (fabsl(total_log - best_log) < 1e-12) {
                // Tie: compare mod values... but we actually want the LARGER divisor.
                // When logs are equal within precision, we need true comparison.
                // Use the fact that if log(a*b) > log(c*d), then a*b > c*d.
                // If they're truly equal in log space within precision, we can't distinguish.
                // Try to pick the larger mod value (heuristic, since actual value matters mod 10^16).
                // Actually for the correct answer, the maximum is unique (no exact ties).
                // So one of the candidates will have strictly larger log.
                long long candidate = mulmod(sortedModA[j], val);
                if (candidate > best_mod) best_mod = candidate;
            }
        }
    }

    cout << best_mod << endl;
    return 0;
}

Python project_euler/problem_266/solution.py

"""
Problem 266: Pseudo Square Root

Find the largest divisor d of the product P of all primes below 190
such that d <= sqrt(P). Give the last 16 digits of d.

Approach: Meet-in-the-middle on the 42 primes.
Split into two halves, enumerate subset products for each,
then for each product in the second half, binary search the first half
for the largest complementary product that keeps d <= sqrt(P).
"""

from math import log
from bisect import bisect_right

def solve():
    # Sieve primes below 190
    primes = []
    for i in range(2, 190):
        if all(i % j != 0 for j in range(2, int(i**0.5) + 1)):
            primes.append(i)

    n = len(primes)  # 42
    half = n // 2     # 21
    MOD = 10**16

    # log(sqrt(P))
    log_sqrtP = sum(log(p) for p in primes) / 2.0

    # Generate subset products for first half
    sz1 = 1 << half
    groupA = []
    for mask in range(sz1):
        lg = 0.0
        val = 1
        for i in range(half):
            if mask & (1 << i):
                lg += log(primes[i])
                val = (val * primes[i]) % MOD
        groupA.append((lg, val))

    # Sort by log
    groupA.sort()
    sorted_logs = [x[0] for x in groupA]
    sorted_mods = [x[1] for x in groupA]

    best_log = -1.0
    best_mod = 0

    # For each subset of second half
    other = n - half
    sz2 = 1 << other
    for mask in range(sz2):
        lg = 0.0
        val = 1
        for i in range(other):
            if mask & (1 << i):
                lg += log(primes[half + i])
                val = (val * primes[half + i]) % MOD

        target = log_sqrtP - lg
        if target < -1e-15:
            continue

        # Binary search: largest index where sorted_logs[idx] <= target
        idx = bisect_right(sorted_logs, target + 1e-15) - 1
        if idx < 0:
            continue

        # Check a few neighbors for float safety
        for j in range(max(0, idx - 2), min(sz1, idx + 3)):
            if sorted_logs[j] > target + 1e-15:
                continue
            total_log = sorted_logs[j] + lg
            if total_log > best_log + 1e-15:
                best_log = total_log
                best_mod = (sorted_mods[j] * val) % MOD
            elif abs(total_log - best_log) < 1e-12:
                candidate = (sorted_mods[j] * val) % MOD
                if candidate > best_mod:
                    best_mod = candidate

    print(best_mod)

solve()