Problem 387: Harshad Numbers

Mathematical Analysis

The key observation is that right truncatable Harshad numbers (RTHN) can be built incrementally. If $n$ is an RTHN with digit sum $s$, then $10n + d$ (for digit $d \in \{0,\dots,9\}$) is an RTHN if and only if $(s+d) \mid (10n+d)$.

The algorithm proceeds by breadth-first generation of all RTHN up to $10^{13}$ (since appending one more digit to get the prime can reach up to $10^{14} - 1$). At each RTHN, we check whether it is strong (i.e., $n/s(n)$ is prime). For each strong RTHN, we append digits $1, 3, 7, 9$ (the only digits that can make the result prime) and test primality.

Derivation / Algorithm

1. Seed: Single-digit numbers $1, 2, \dots, 9$ are all RTHN (each is trivially divisible by itself).
2. Expand: For each RTHN $n$ with digit sum $s$ and $n < 10^{13}$, try appending each digit $d \in \{0,\dots,9\}$. If $(s + d) \mid (10n + d)$, then $10n + d$ is a new RTHN with digit sum $s + d$.
3. Check strong: For each RTHN $n$, if $n / s(n)$ is prime, mark $n$ as a strong RTHN.
4. Collect primes: For each strong RTHN $n$ (with $n < 10^{13}$), try $p = 10n + d$ for $d \in \{1, 3, 7, 9\}$. If $p < 10^{14}$ and $p$ is prime, add $p$ to the answer.

Proof of Correctness

• Completeness: Every RTHN of $k$ digits has an RTHN prefix of $k-1$ digits, so the BFS from single-digit seeds generates all RTHN.
• Primality of strong RTHN primes: We only count a prime $p$ if its prefix $\lfloor p/10 \rfloor$ is simultaneously right truncatable Harshad and strong Harshad, matching the problem definition exactly.
• Verification: The algorithm yields a sum of 90619 for the limit $10^4$, matching the given test case.

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

• Time: The number of RTHN up to $10^{13}$ is relatively small (on the order of tens of thousands). For each, we perform a constant number of divisibility checks and at most 4 primality tests. Primality testing of numbers up to $10^{14}$ via deterministic Miller-Rabin with appropriate witnesses runs in $O(\log^2 n)$. Overall: effectively $O(R \cdot \log^2 N)$ where $R$ is the count of RTHN and $N = 10^{14}$.
• Space: $O(R)$ to store the current frontier of RTHN.

Answer

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

typedef long long ll;
typedef unsigned long long ull;

// Modular exponentiation: a^e mod m (handles large intermediates via __int128)
ll mod_pow(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1)
            result = (__int128)result * base % mod;
        base = (__int128)base * base % mod;
        exp >>= 1;
    }
    return result;
}

// Deterministic Miller-Rabin for n < 3.317e14
bool is_prime(ll n) {
    if (n < 2) return false;
    if (n < 4) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    ll d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    for (ll a : {2, 3, 5, 7, 11, 13, 17}) {
        if (a >= n) continue;
        ll x = mod_pow(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int i = 0; i < r - 1; i++) {
            x = (__int128)x * x % n;
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const ll LIMIT = 100000000000000LL; // 10^14
    const ll PREFIX_LIMIT = LIMIT / 10;

    ll total = 0;
    // (number, digit_sum)
    vector<pair<ll, int>> frontier;
    for (int d = 1; d <= 9; d++)
        frontier.push_back({d, d});

    while (!frontier.empty()) {
        vector<pair<ll, int>> next_frontier;
        for (auto& [n, s] : frontier) {
            // Check if strong Harshad: n/s is prime
            if (n % s == 0 && is_prime(n / s)) {
                for (int d : {1, 3, 7, 9}) {
                    ll candidate = 10 * n + d;
                    if (candidate < LIMIT && is_prime(candidate))
                        total += candidate;
                }
            }
            // Expand: append digits 0..9
            for (int d = 0; d <= 9; d++) {
                ll new_n = 10 * n + d;
                int new_s = s + d;
                if (new_n < PREFIX_LIMIT && new_s > 0 && new_n % new_s == 0)
                    next_frontier.push_back({new_n, new_s});
            }
        }
        frontier = move(next_frontier);
    }

    cout << total << endl;
    return 0;
}

"""
Problem 387: Harshad Numbers
Find the sum of all strong, right truncatable Harshad primes less than 10^14.
"""

def is_prime(n: int) -> bool:
    """Deterministic Miller-Rabin primality test for n < 3.317e14."""
    if n < 2:
        return False
    if n < 4:
        return True
    if n % 2 == 0 or n % 3 == 0:
        return False
    # Small primes check
    if n < 25:
        return n in {2, 3, 5, 7, 11, 13, 17, 19, 23}
    # Write n-1 = d * 2^r
    d, r = n - 1, 0
    while d % 2 == 0:
        d //= 2
        r += 1
    # Witnesses sufficient for n < 3.317e14
    for a in [2, 3, 5, 7, 11, 13, 17]:
        if a >= n:
            continue
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

def solve(limit: int = 10**14):
    """
    BFS over right truncatable Harshad numbers (RTHN).
    For each strong RTHN, try appending a digit to form a prime.
    """
    total = 0
    # frontier: list of (number, digit_sum)
    # Seeds: single-digit Harshad numbers 1..9
    frontier = [(d, d) for d in range(1, 10)]
    prefix_limit = limit // 10  # RTHN must be < limit/10

    while frontier:
        next_frontier = []
        for n, s in frontier:
            # Check if n is a strong Harshad number
            if n % s == 0 and is_prime(n // s):
                # Try appending digits that could make a prime
                for d in [1, 3, 7, 9]:
                    candidate = 10 * n + d
                    if candidate < limit and is_prime(candidate):
                        total += candidate

            # Expand: try appending each digit 0..9
            for d in range(10):
                new_n = 10 * n + d
                new_s = s + d
                if new_n < prefix_limit and new_s > 0 and new_n % new_s == 0:
                    next_frontier.append((new_n, new_s))

        frontier = next_frontier

    return total

# Compute and verify
answer_small = solve(10_000)
assert answer_small == 90619, f"Small test failed: {answer_small}"

answer = solve(10**14)
assert answer == 696067597313468

print(answer)