Problem 110: Diophantine Reciprocals II

Mathematical Foundation

Theorem 1. (Algebraic Transformation — same as Problem 108) The equation $\frac{1}{x} + \frac{1}{y} = \frac{1}{n}$ with positive integers $x, y$ is equivalent to $(x-n)(y-n) = n^2$, and the number of solutions with $x \leq y$ is $f(n) = \lceil d(n^2)/2 \rceil$.

Proof. See Problem 108, Theorems 1 and 2. $\square$

Theorem 2. (Divisor Count of a Square) If $n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$, then $d(n^2) = \prod_{i=1}^{k}(2a_i + 1)$.

Proof. See Problem 108, Theorem 3. $\square$

Theorem 3. (Required Bound) We need $f(n) > 4{,}000{,}000$, i.e., $\lceil d(n^2)/2 \rceil > 4{,}000{,}000$.

Since $d(n^2) = \prod(2a_i + 1)$ is always odd, $\lceil d(n^2)/2 \rceil = (d(n^2) + 1)/2$. The condition $(d(n^2) + 1)/2 > 4{,}000{,}000$ gives $d(n^2) > 7{,}999{,}999$, so $d(n^2) \geq 8{,}000{,}001$ (since $d(n^2)$ is odd).

Proof. $d(n^2)$ is a product of odd numbers $(2a_i + 1)$, hence odd. For odd $m$, $\lceil m/2 \rceil = (m+1)/2$. The inequality $(m+1)/2 > 4{,}000{,}000$ gives $m > 7{,}999{,}999$, and since $m$ is odd, $m \geq 8{,}000{,}001$. $\square$

Theorem 4. (Optimization as a Discrete Minimization Problem) Minimize $n = \prod p_i^{a_i}$ (equivalently, minimize $\log n = \sum a_i \log p_i$) subject to $\prod(2a_i + 1) \geq 8{,}000{,}001$ and $a_1 \geq a_2 \geq \cdots \geq a_k \geq 1$, where $p_i$ is the $i$-th prime.

Proof. The constraint $a_1 \geq a_2 \geq \cdots \geq a_k$ is necessary for optimality by the same argument as Problem 108 Lemma 1: swapping a larger exponent to a smaller prime strictly decreases $n$ while preserving $d(n^2)$. $\square$

Lemma 1. (Search Space Bounds) The number of primes $k$ satisfies $k \leq \lfloor \log_3(8{,}000{,}001) \rfloor = 14$, and the maximum exponent satisfies $a_1 \leq \lfloor (8{,}000{,}001 - 1)/2 \rfloor$ but in practice $a_1 \leq 14$ suffices (since $3^{14} = 4{,}782{,}969 < 8{,}000{,}001 < 3^{15}$, using more primes with exponent 1 is more efficient than increasing one exponent).

Proof. With $k$ primes all having exponent $\geq 1$, $d(n^2) \geq 3^k$. For $3^k \geq 8{,}000{,}001$, we need $k \geq 15$. However, some primes may have exponent $> 1$, so $k$ can be smaller. The bound $k \leq 14$ follows from $\log_3(8{,}000{,}001) \approx 14.5$. For the maximum exponent: if $a_1 = m$ and all other exponents are 1, then $d(n^2) = (2m+1) \cdot 3^{k-1}$, and it suffices to have $m$ small with $k$ large. The practical bound $a_1 \leq 14$ is sufficient. $\square$

Theorem 5. (Solution) The minimum $n$ is:

with exponent sequence $(3, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1)$ giving:

and $f(n) = \lceil 8{,}037{,}225/2 \rceil = 4{,}018{,}613 > 4{,}000{,}000$.

Proof. The exponent sequence $(3, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1)$ is non-increasing with 12 primes. The divisor product is $7^2 \cdot 5^2 \cdot 3^8 = 8{,}037{,}225 \geq 8{,}000{,}001$. Assigning exponents to the first 12 primes in decreasing order: $2^3 \cdot 3^3 \cdot 5^2 \cdot 7^2 \cdot 11 \cdot 13 \cdot 17 \cdot 19 \cdot 23 \cdot 29 \cdot 31 \cdot 37 = 9{,}350{,}130{,}049{,}860{,}600$.

To verify minimality: any alternative exponent sequence with $\prod(2a_i+1) \geq 8{,}000{,}001$ yields $n \geq 9{,}350{,}130{,}049{,}860{,}600$. This is confirmed by exhaustive search over all feasible non-increasing exponent sequences (bounded by $k \leq 15$ primes and $a_1 \leq 14$), computing $n$ for each and retaining the minimum. $\square$

Editorial

This is the same divisor-count minimization as Problem 108, but the target is now so large that comparing candidate integers directly during the search becomes awkward. The structure of the optimal solution is unchanged: we still assign a non-increasing exponent sequence to successive primes, because moving a larger exponent onto a smaller prime always makes $n$ smaller.

The implementation therefore performs the same recursive search over exponent patterns, but it compares candidates in log-space. Tracking $\log n = \sum a_i \log p_i$ preserves the ordering of candidate values while avoiding overflow during recursion. Once the best exponent pattern has been identified, the actual integer $n$ is reconstructed from those stored exponents.

Pseudocode

Store the primes and their logarithms, and begin with no exponents chosen.
Recursively assign an exponent to each successive prime, never exceeding the previous exponent.

At each recursive call:
    If the current divisor product has reached 8,000,001, compare the current log(n) with the best one seen so far and keep the better exponent pattern.
    Otherwise try the next exponent values in order, update both the divisor product and log(n), and abandon any branch that is already no better than the current best.
    After returning from a branch, remove the exponent and continue with the next choice.

After the search finishes, rebuild n from the stored best exponents and return it.

Complexity Analysis

• Time: The recursive search explores non-increasing exponent sequences with aggressive pruning ($n < \text{best\_n}$ and $\prod(2a_i+1) \geq 8{,}000{,}001$). The number of feasible sequences is bounded by the number of partitions of $\lceil\log_3(8{,}000{,}001)\rceil \approx 15$ into at most 15 parts, weighted by decreasing exponent constraints. In practice, the search visits on the order of $10^4$ nodes and completes in milliseconds.
• Space: $O(k)$ for the recursion stack, where $k \leq 15$.

Answer

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

/*
 * Problem 110: Diophantine Reciprocals II
 *
 * Same as Problem 108 but threshold is 4,000,000.
 * Need d(n^2) >= 8,000,001.
 *
 * We use logarithms to avoid overflow when comparing n values,
 * since n can be very large. We track log(n) and reconstruct
 * the actual n at the end using __int128 or careful multiplication.
 */

const int NPRIMES = 20;
int primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71};
double log_primes[NPRIMES];
double best_log_n;
vector<int> best_exponents;

void dfs(int idx, int max_exp, double log_n, long long div_count, vector<int>& exponents) {
    if (div_count >= 8000001LL) {
        if (log_n < best_log_n) {
            best_log_n = log_n;
            best_exponents = exponents;
        }
        return;
    }
    if (idx >= NPRIMES) return;
    if (log_n >= best_log_n) return;

    for (int e = 1; e <= max_exp; e++) {
        double new_log = log_n + e * log_primes[idx];
        if (new_log >= best_log_n) break;
        exponents.push_back(e);
        dfs(idx + 1, e, new_log, div_count * (2 * e + 1), exponents);
        exponents.pop_back();
    }
}

int main() {
    for (int i = 0; i < NPRIMES; i++)
        log_primes[i] = log((double)primes[i]);

    best_log_n = 1e18;
    vector<int> exponents;
    dfs(0, 40, 0.0, 1, exponents);

    // Reconstruct the actual n using the best exponents
    // Use __int128 or careful computation
    // Since answer fits in long long (9350130049860600), we can use long long
    long long n = 1;
    for (int i = 0; i < (int)best_exponents.size(); i++) {
        for (int j = 0; j < best_exponents[i]; j++) {
            n *= primes[i];
        }
    }

    cout << n << endl;
    return 0;
}

"""
Problem 110: Diophantine Reciprocals II

Find the least n such that the number of distinct solutions of
1/x + 1/y = 1/n exceeds 4,000,000.

Same approach as Problem 108 but with threshold 8,000,001 for d(n^2).
Uses log-space comparison to handle large n values.
"""

import math

PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]
LOG_PRIMES = [math.log(p) for p in PRIMES]
TARGET = 8_000_001

best_log_n = float('inf')
best_exponents = []

def search(idx, max_exp, log_n, div_count, exponents):
    global best_log_n, best_exponents

    if div_count >= TARGET:
        if log_n < best_log_n:
            best_log_n = log_n
            best_exponents = exponents[:]
        return

    if idx >= len(PRIMES):
        return

    if log_n >= best_log_n:
        return

    for e in range(1, max_exp + 1):
        new_log = log_n + e * LOG_PRIMES[idx]
        if new_log >= best_log_n:
            break
        exponents.append(e)
        search(idx + 1, e, new_log, div_count * (2 * e + 1), exponents)
        exponents.pop()

search(0, 40, 0.0, 1, [])

# Reconstruct n
n = 1
for i, e in enumerate(best_exponents):
    n *= PRIMES[i] ** e

print(n)

# Verification
if __name__ == "__main__":
    d_n2 = 1
    for e in best_exponents:
        d_n2 *= (2 * e + 1)
    solutions = (d_n2 + 1) // 2

    factors = {PRIMES[i]: best_exponents[i] for i in range(len(best_exponents))}
    print(f"n = {n}")
    print(f"Factorization: {factors}")
    print(f"d(n^2) = {d_n2}")
    print(f"Solutions = {solutions}")
    print(f"Exceeds 4,000,000: {solutions > 4_000_000}")