Problem 387: Harshad Numbers

View on Project Euler

Project Euler Problem 387 Solution

EulerSolve provides an optimized solution for Project Euler Problem 387, Harshad Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want the sum of all primes \(p \lt L\) such that deleting the last decimal digit of \(p\) leaves a strong right-truncatable Harshad number . In the local solution files the default bound is \(L = 10^{14}\). Let \(s(n)\) denote the sum of the decimal digits of \(n\). A positive integer \(n\) is a Harshad number if $$s(n)\mid n.$$ It is strong if $$\frac{n}{s(n)}$$ is prime. It is right-truncatable if every non-empty decimal prefix is Harshad. Writing \(n=\overline{d_1d_2\cdots d_k}\), each prefix \(\overline{d_1d_2\cdots d_i}\) for \(1 \le i \le k\) must satisfy the Harshad divisibility condition. Every candidate answer prime has a unique decomposition $$p = 10h + d,\qquad d = p \bmod 10,\qquad h=\left\lfloor\frac{p}{10}\right\rfloor.$$ The problem therefore reduces to generating exactly those prefixes \(h\) that are strong and right-truncatable Harshad numbers, then checking whether one more digit turns them into a prime below \(L\). Mathematical Approach Step 1: Search the prefix, not the final prime A brute-force approach would iterate over primes \(p \lt L\), remove the last digit, and test the whole prefix chain. The repository code does the opposite: it constructs only valid Harshad prefixes and never explores arbitrary integers. All one-digit numbers \(1,\dots,9\) are Harshad because for a single decimal digit we have \(s(n)=n\), hence \(n/s(n)=1\)....

Detailed mathematical approach

Problem Summary

We want the sum of all primes \(p \lt L\) such that deleting the last decimal digit of \(p\) leaves a strong right-truncatable Harshad number. In the local solution files the default bound is \(L = 10^{14}\).

Let \(s(n)\) denote the sum of the decimal digits of \(n\). A positive integer \(n\) is a Harshad number if

$$s(n)\mid n.$$

It is strong if

$$\frac{n}{s(n)}$$

is prime. It is right-truncatable if every non-empty decimal prefix is Harshad. Writing \(n=\overline{d_1d_2\cdots d_k}\), each prefix \(\overline{d_1d_2\cdots d_i}\) for \(1 \le i \le k\) must satisfy the Harshad divisibility condition.

Every candidate answer prime has a unique decomposition

$$p = 10h + d,\qquad d = p \bmod 10,\qquad h=\left\lfloor\frac{p}{10}\right\rfloor.$$

The problem therefore reduces to generating exactly those prefixes \(h\) that are strong and right-truncatable Harshad numbers, then checking whether one more digit turns them into a prime below \(L\).

Mathematical Approach

Step 1: Search the prefix, not the final prime

A brute-force approach would iterate over primes \(p \lt L\), remove the last digit, and test the whole prefix chain. The repository code does the opposite: it constructs only valid Harshad prefixes and never explores arbitrary integers.

All one-digit numbers \(1,\dots,9\) are Harshad because for a single decimal digit we have \(s(n)=n\), hence \(n/s(n)=1\). These values are the roots of the search tree.

Step 2: Extend prefixes with a digit-sum recurrence

Suppose the current node is a Harshad number \(n\) with digit sum \(s=s(n)\). Appending a digit \(d\in\{0,1,\dots,9\}\) gives

$$n' = 10n + d,\qquad s(n') = s + d.$$

The child is useful only if it is Harshad, namely

$$10n + d \equiv 0 \pmod{s+d}.$$

This test is exactly what the DFS applies before recursing. Because every child is created from an already valid Harshad parent, right-truncatability is automatic by induction: each visited node has a Harshad parent chain all the way back to a one-digit root.

Step 3: Identify strong Harshad prefixes immediately

Whenever a visited prefix \(n \ge 10\) is Harshad, the code computes

$$q=\frac{n}{s(n)}.$$

If \(q\) is prime, then \(n\) is strong. This is the decisive filter, because a final answer prime must come from a strong prefix. One-digit prefixes are excluded from this stage in the code since they would only produce \(q=1\), which is not prime.

Step 4: Only four last digits can produce a prime

For any prime \(p \gt 10\), the last digit cannot be even and cannot be \(5\). Therefore every strong prefix \(h\) produces at most four relevant prime candidates:

$$p = 10h + d,\qquad d \in \{1,3,7,9\}.$$

The implementations test each such \(p\) for the two remaining conditions: \(p \lt L\) and \(p\) is prime.

If both hold, \(p\) is added to the total. This characterization is complete and duplicate-free, because decimal truncation is unique: every valid prime determines exactly one prefix \(h=\lfloor p/10 \rfloor\).

Step 5: Bound the recursion by the largest useful prefix

If \(p=10h+d \lt L\), then necessarily

$$h \le \left\lfloor\frac{L-1}{10}\right\rfloor.$$

The C++ solution stores this value in harshad_limit. Any Harshad number larger than this can never be the prefix of a valid answer, so exploring beyond it is pointless.

The recursive stop test in the code is slightly sharper. If

$$n \gt \left\lfloor\frac{\left\lfloor(L-1)/10\right\rfloor}{10}\right\rfloor,$$

then every child \(10n+d\) already exceeds harshad_limit. The current node may still produce final primes \(10n+d\), but it cannot produce deeper Harshad prefixes that will later matter.

Worked Example

Take \(n=18\). Its digit sum is \(s(18)=9\), and \(18/9=2\), so \(18\) is a strong Harshad number. Since \(1\) and \(18\) are both Harshad, it is also right-truncatable.

The only possible prime endings are

$$181,\ 183,\ 187,\ 189.$$

Among these, \(181\) is prime, \(183\) and \(189\) are divisible by \(3\), and \(187=11\cdot 17\). Hence \(181\) is one of the required primes.

The local C++ program also validates the sample checkpoint

$$\sum_{p \lt 10^4} p = 90619,$$

and then cross-checks the optimized solver against a brute-force routine at \(L=10^6\). Those checks are consistent with the derivation above.

How the Code Works

The three local implementations use the same core state: the current prefix \(n\) and its digit sum \(s(n)\). Carrying the digit sum forward via \(s(n')=s(n)+d\) makes each tree transition constant time instead of recomputing all decimal digits from scratch.

The C++ file exposes solve(limit), optional argument parsing with --limit=..., and a --skip-checkpoints flag. Its DFS function dfs_harshad returns the sum contributed by one subtree. For primality it first removes small prime divisors up to \(37\), then runs deterministic Miller-Rabin with bases

$$\{2,3,5,7,11,13,17\}.$$

To avoid overflow in modular multiplication, C++ uses __uint128_t.

The Python file implements the same recurrence and the same Miller-Rabin base set, but accumulates the answer in a nonlocal variable inside dfs. The Java file mirrors the same DFS structure in dfsHarshad and delegates modular exponentiation to BigInteger.modPow. Mathematically all three versions are identical.

Complexity Analysis

Let \(H(L)\) be the number of right-truncatable Harshad prefixes actually visited, and let \(S(L)\) be the number of those prefixes that are strong. The DFS work per visited node is \(O(1)\): one digit-sum update and one divisibility test for each appended digit.

Each strong prefix triggers at most four primality tests, so a precise description is

$$O\bigl(H(L) + 4S(L)\,T_{\mathrm{prime}}\bigr),$$

where \(T_{\mathrm{prime}}\) is the cost of one Miller-Rabin test in the 64-bit range used here. Since \(S(L)\le H(L)\), this is often summarized as

$$O\bigl(H(L)\,T_{\mathrm{prime}}\bigr).$$

The recursion depth is the number of decimal digits of the largest explored prefix, so the auxiliary stack space is \(O(\log_{10} L)\). This is far smaller than scanning all integers or all primes below \(L\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=387
  2. Harshad numbers: https://en.wikipedia.org/wiki/Harshad_number
  3. Miller-Rabin primality test: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
  4. Deterministic 64-bit primality testing reference: https://cp-algorithms.com/algebra/primality_tests.html

Problem 387 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 100'000'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0ULL;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(ch - '0');
    }
    value = parsed;
    return 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;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 100ULL;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    u64 result = 1ULL;
    base %= mod;
    while (exp > 0ULL) {
        if (exp & 1ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

bool is_prime(const u64 n) {
    if (n < 2ULL) {
        return false;
    }
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0ULL) {
            return false;
        }
    }

    u64 d = n - 1ULL;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1ULL;
        ++s;
    }

    for (u64 a : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL}) {
        if (a >= n) {
            continue;
        }
        u64 x = pow_mod(a, d, n);
        if (x == 1ULL || x == n - 1ULL) {
            continue;
        }
        bool witness = true;
        for (int r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1ULL) {
                witness = false;
                break;
            }
        }
        if (witness) {
            return false;
        }
    }
    return true;
}

u64 dfs_harshad(const u64 n, const int digit_sum, const u64 harshad_limit, const u64 prime_limit) {
    u64 acc = 0ULL;

    if (n >= 10ULL && n % static_cast<u64>(digit_sum) == 0ULL) {
        const u64 q = n / static_cast<u64>(digit_sum);
        if (is_prime(q)) {
            for (int d : {1, 3, 7, 9}) {
                const u64 cand = n * 10ULL + static_cast<u64>(d);
                if (cand < prime_limit && is_prime(cand)) {
                    acc += cand;
                }
            }
        }
    }

    if (n > harshad_limit / 10ULL) {
        return acc;
    }

    for (int d = 0; d <= 9; ++d) {
        const u64 next = n * 10ULL + static_cast<u64>(d);
        const int next_sum = digit_sum + d;
        if (next_sum != 0 && next % static_cast<u64>(next_sum) == 0ULL) {
            acc += dfs_harshad(next, next_sum, harshad_limit, prime_limit);
        }
    }
    return acc;
}

u64 solve(const u64 limit) {
    const u64 harshad_limit = (limit - 1ULL) / 10ULL;
    u64 total = 0ULL;
    for (int d = 1; d <= 9; ++d) {
        total += dfs_harshad(static_cast<u64>(d), d, harshad_limit, limit);
    }
    return total;
}

u64 brute_small(const u64 limit) {
    auto digit_sum = [](u64 n) {
        int s = 0;
        while (n > 0ULL) {
            s += static_cast<int>(n % 10ULL);
            n /= 10ULL;
        }
        return s;
    };
    auto is_right_trunc_harshad = [&](u64 n) {
        while (n > 0ULL) {
            const int s = digit_sum(n);
            if (s == 0 || n % static_cast<u64>(s) != 0ULL) {
                return false;
            }
            n /= 10ULL;
        }
        return true;
    };

    u64 sum = 0ULL;
    for (u64 p = 2ULL; p < limit; ++p) {
        if (!is_prime(p)) {
            continue;
        }
        const u64 h = p / 10ULL;
        if (h < 10ULL) {
            continue;
        }
        if (!is_right_trunc_harshad(h)) {
            continue;
        }
        const int s = digit_sum(h);
        if (h % static_cast<u64>(s) == 0ULL && is_prime(h / static_cast<u64>(s))) {
            sum += p;
        }
    }
    return sum;
}

bool run_checkpoints() {
    if (solve(10'000ULL) != 90'619ULL) {
        std::cerr << "Checkpoint failed: limit=10000 sample" << '\n';
        return false;
    }
    if (solve(1'000'000ULL) != brute_small(1'000'000ULL)) {
        std::cerr << "Checkpoint failed: brute-force cross-check for limit=1e6" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

import sys
sys.setrecursionlimit(100000)

def solve():
    LIMIT = 100_000_000_000_000

    def pow_mod(base, exp, mod):
        result = 1
        base %= mod
        while exp > 0:
            if exp & 1:
                result = result * base % mod
            base = base * base % mod
            exp >>= 1
        return result

    def is_prime(n):
        if n < 2:
            return False
        for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
            if n == p:
                return True
            if n % p == 0:
                return False
        d = n - 1
        s = 0
        while d % 2 == 0:
            d //= 2
            s += 1
        for a in [2, 3, 5, 7, 11, 13, 17]:
            if a >= n:
                continue
            x = pow_mod(a, d, n)
            if x == 1 or x == n - 1:
                continue
            witness = True
            for _ in range(1, s):
                x = x * x % n
                if x == n - 1:
                    witness = False
                    break
            if witness:
                return False
        return True

    harshad_limit = (LIMIT - 1) // 10
    total = 0

    def dfs(n, digit_sum):
        nonlocal total
        if n >= 10 and n % digit_sum == 0:
            q = n // digit_sum
            if is_prime(q):
                for d in [1, 3, 7, 9]:
                    cand = n * 10 + d
                    if cand < LIMIT and is_prime(cand):
                        total += cand
        if n > harshad_limit // 10:
            return
        for d in range(10):
            nxt = n * 10 + d
            ns = digit_sum + d
            if ns != 0 and nxt % ns == 0:
                dfs(nxt, ns)

    for d in range(1, 10):
        dfs(d, d)

    return str(total)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;

public class Euler387 {
    private static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }
        BigInteger bn = BigInteger.valueOf(n);
        BigInteger bnMinus1 = bn.subtract(BigInteger.ONE);
        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }
        long[] bases = { 2, 3, 5, 7, 11, 13, 17 };
        for (long a : bases) {
            if (a >= n)
                continue;
            BigInteger ba = BigInteger.valueOf(a);
            BigInteger x = ba.modPow(BigInteger.valueOf(d), bn);
            if (x.equals(BigInteger.ONE) || x.equals(bnMinus1))
                continue;
            boolean witness = true;
            for (int r = 1; r < s; r++) {
                x = x.multiply(x).mod(bn);
                if (x.equals(bnMinus1)) {
                    witness = false;
                    break;
                }
            }
            if (witness)
                return false;
        }
        return true;
    }

    private static long dfsHarshad(long n, int digitSum, long harshadLimit, long primeLimit) {
        long acc = 0;
        if (n >= 10 && n % digitSum == 0) {
            long q = n / digitSum;
            if (isPrime(q)) {
                int[] ds = { 1, 3, 7, 9 };
                for (int d : ds) {
                    long cand = n * 10 + d;
                    if (cand < primeLimit && isPrime(cand)) {
                        acc += cand;
                    }
                }
            }
        }

        if (n > harshadLimit / 10)
            return acc;

        for (int d = 0; d <= 9; d++) {
            long next = n * 10 + d;
            int nextSum = digitSum + d;
            if (nextSum != 0 && next % nextSum == 0) {
                acc += dfsHarshad(next, nextSum, harshadLimit, primeLimit);
            }
        }
        return acc;
    }

    public static String solve() {
        long limit = 100000000000000L;
        long harshadLimit = (limit - 1) / 10;
        long total = 0;
        for (int d = 1; d <= 9; d++) {
            total += dfsHarshad(d, d, harshadLimit, limit);
        }
        return String.valueOf(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}