Problem 234: Semidivisible Numbers

Mathematical Foundation

Theorem (Interval Decomposition). For consecutive primes $p < q$, every integer $n$ with $p^2 < n < q^2$ satisfies $\mathrm{lps}(n) = p$ and $\mathrm{ups}(n) = q$.

Proof. Since $p$ and $q$ are consecutive primes and $p^2 < n < q^2$, we have $p < \sqrt{n} < q$. As $p$ is the largest prime $\leq p < \sqrt{n}$ and there is no prime between $p$ and $q$, $p$ is the largest prime $\leq \sqrt{n}$, so $\mathrm{lps}(n) = p$. Similarly, $q$ is the smallest prime $\geq \sqrt{n}$, so $\mathrm{ups}(n) = q$. $\square$

Lemma (Semidivisibility Criterion). For $n$ in the interval $(p^2, q^2)$ with consecutive primes $p < q$, $n$ is semidivisible if and only if exactly one of $p \mid n$ or $q \mid n$ holds.

Proof. By the theorem, $\mathrm{lps}(n) = p$ and $\mathrm{ups}(n) = q$. By definition, $n$ is semidivisible when exactly one of these divides $n$. $\square$

Lemma (Arithmetic Sum Formula). The sum of semidivisible numbers in the open interval $(p^2, \min(q^2, L))$ equals

where $\Sigma_d$ denotes the sum of all multiples of $d$ in the interval.

Proof. Let $A$ be the set of multiples of $p$ and $B$ the set of multiples of $q$ in the interval. Semidivisible numbers form the symmetric difference $A \triangle B = (A \setminus B) \cup (B \setminus A)$. Using $\sum_{A \triangle B} = \sum_A + \sum_B - 2\sum_{A \cap B}$ and $A \cap B$ is the set of multiples of $\mathrm{lcm}(p,q) = pq$ (since $p, q$ are distinct primes), the result follows. $\square$

Each $\Sigma_d$ is computed as the sum of an arithmetic progression: for multiples of $d$ in $(a, b)$, the first term is $d \cdot \lceil (a+1)/d \rceil$, the last is $d \cdot \lfloor (b-1)/d \rfloor$, and the sum is $d \cdot (k_1 + k_2)(k_2 - k_1 + 1)/2$.

Editorial

A number n is semidivisible if exactly one of lps(n), ups(n) divides n, where lps(n) = largest prime <= sqrt(n), ups(n) = smallest prime >= sqrt(n). We sum of multiples of d in open interval (lo, hi). Finally, sum of multiples of d in (lo, hi) exclusive.

Pseudocode

L = 999966663333
Sum of multiples of d in open interval (lo, hi)
Sum of multiples of d in (lo, hi) exclusive

Complexity Analysis

• Time: $O(\sqrt{L} \log \log \sqrt{L})$ for the sieve, plus $O(\pi(\sqrt{L}))$ for iterating over consecutive prime pairs. Each pair requires $O(1)$ arithmetic. Total: $O(\sqrt{L} \log \log \sqrt{L})$.
• Space: $O(\sqrt{L})$ for the sieve.

Answer

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

typedef long long ll;
typedef __int128 lll;

// Sum of multiples of d in open interval (lo, hi), i.e., lo < n < hi and d | n
// Returns sum as __int128 to avoid overflow
lll sum_multiples(ll d, ll lo, ll hi) {
    // First multiple of d greater than lo
    ll first_k = lo / d + 1;
    // Last multiple of d less than hi
    ll last_k = (hi - 1) / d;
    if (first_k > last_k) return 0;
    // Sum = d * (first_k + first_k+1 + ... + last_k) = d * (last_k - first_k + 1) * (first_k + last_k) / 2
    lll cnt = last_k - first_k + 1;
    lll s = cnt * (first_k + last_k) / 2;
    return (lll)d * s;
}

int main() {
    const ll LIMIT = 999966663333LL;
    const int SIEVE_SIZE = 1100000; // primes up to ~1.1M covers sqrt(LIMIT) ~ 999983

    // Sieve
    vector<bool> is_prime(SIEVE_SIZE + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; (ll)i * i <= SIEVE_SIZE; i++)
        if (is_prime[i])
            for (int j = i * i; j <= SIEVE_SIZE; j += i)
                is_prime[j] = false;

    vector<ll> primes;
    for (int i = 2; i <= SIEVE_SIZE; i++)
        if (is_prime[i])
            primes.push_back(i);

    lll total = 0;

    for (int i = 0; i + 1 < (int)primes.size(); i++) {
        ll p = primes[i];
        ll q = primes[i + 1];
        ll lo = p * p;    // exclusive lower bound
        ll hi = q * q;    // exclusive upper bound (or LIMIT+1)

        if (lo >= LIMIT) break;
        if (hi > LIMIT + 1) hi = LIMIT + 1;

        // Sum of semidivisible in (lo, hi):
        // divisible by exactly one of p, q
        // = S_p + S_q - 2*S_{pq}
        lll sp = sum_multiples(p, lo, hi);
        lll sq = sum_multiples(q, lo, hi);
        lll spq = sum_multiples(p * q, lo, hi);
        total += sp + sq - 2 * spq;
    }

    // Print __int128
    // Convert to string
    ll hi_part = (ll)(total / 1000000000000000000LL);
    ll lo_part = (ll)(total % 1000000000000000000LL);
    if (hi_part > 0)
        printf("%lld%018lld\n", hi_part, lo_part);
    else
        printf("%lld\n", lo_part);

    return 0;
}

"""
Problem 234: Semidivisible Numbers

Find the sum of all semidivisible numbers not exceeding 999966663333.
A number n is semidivisible if exactly one of lps(n), ups(n) divides n,
where lps(n) = largest prime <= sqrt(n), ups(n) = smallest prime >= sqrt(n).
"""

def sieve_primes(n):
    is_p = bytearray(b'\x01') * (n + 1)
    is_p[0] = is_p[1] = 0
    for i in range(2, int(n**0.5) + 1):
        if is_p[i]:
            is_p[i*i::i] = bytearray(len(is_p[i*i::i]))
    return [i for i in range(2, n + 1) if is_p[i]]

def sum_multiples(d, lo, hi):
    """Sum of multiples of d in open interval (lo, hi)."""
    first_k = lo // d + 1
    last_k = (hi - 1) // d
    if first_k > last_k:
        return 0
    cnt = last_k - first_k + 1
    return d * cnt * (first_k + last_k) // 2

def solve():
    LIMIT = 999966663333

    primes = sieve_primes(1100000)

    total = 0
    for i in range(len(primes) - 1):
        p = primes[i]
        q = primes[i + 1]
        lo = p * p
        hi = q * q

        if lo >= LIMIT:
            break
        if hi > LIMIT + 1:
            hi = LIMIT + 1

        sp = sum_multiples(p, lo, hi)
        sq = sum_multiples(q, lo, hi)
        spq = sum_multiples(p * q, lo, hi)
        total += sp + sq - 2 * spq

    print(total)

if __name__ == "__main__":
    solve()