Problem 387: Harshad Numbers

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

Write $s(n)$ for the sum of the decimal digits of $n$. A positive integer $n$ is a Harshad or Niven number when $s(n)\mid n$, right truncatable when repeatedly deleting its last digit leaves a Harshad number at every step, and strong when $n/s(n)$ is prime. A strong, right truncatable Harshad prime is a prime $p$ whose last digit can be deleted to leave a number with all three properties, and we must sum all of them below $L=10^{14}$. The statement illustrates this with $2011$ and gives the checkpoint that they sum to $90619$ below $10^{4}$. ^[1]

The whole solution in one sentence: grow the admissible prefixes one digit at a time, because every prefix of a right truncatable Harshad number is again one, and test only the four endings $1,3,7,9$ of each prefix that turns out to be strong.

Mathematical Approach

1. The definitions on small numbers

Here $s(18)=9$ and $9\mid 18$, so $18$ is Harshad, while $s(19)=10$ and $10\nmid 19$, so $19$ is not. ^[2] Every one-digit number is Harshad, since $s(n)=n$ there, and $0$ never occurs below, its digit sum being $0$. For $201$ the chain $201\to 20\to 2$ has digit sums $3,2,2$, each dividing its own number, so $201$ is right truncatable; and $201/3=67$ is prime, so $201$ is strong too. As $2011$ is prime, it is one of the numbers we must add.

2. Search the prefix, not the prime

Every integer $p\ge 10$ splits uniquely as $p=10h+d$, with $h=\lfloor p/10\rfloor$ and $d=p\bmod 10$. So $p$ qualifies exactly when $p$ is prime, $p\lt L$, and $h$ is a strong right truncatable Harshad number. Each prime has only this one prefix, so enumerating prefixes and extending them by one digit finds every qualifying prime once and only once. That is the gain: about $3\cdot 10^{12}$ primes below $10^{14}$ against only $12753$ prefixes, counted in Section 7.

3. The prefixes form a tree

Right truncatability is hereditary. For $n\ge 10$ put $n'=\lfloor n/10\rfloor$; deleting the last digit of $n$ and continuing is the same as starting the chain at $n'$, so $n$ is a right truncatable Harshad number precisely when $s(n)\mid n$ and $n'$ is one too, and unwinding this ends at a one-digit number, always Harshad. These numbers are therefore the nodes of a tree: the roots are $1,\dots ,9$, and the children of $n$ are those $10n+d$ with $0\le d\le 9$ that are Harshad. Each node $n\ge 10$ has one parent, $\lfloor n/10\rfloor$, so a single depth-first walk visits every prefix once. ^[3]

A child costs one division only, because appending $d$ to $n\ge 1$ gives the decimal expansion of $10n+d$, so

$$s(10n+d)=s(n)+d\qquad (0\le d\le 9),$$

and the child is Harshad exactly when $10n+d\equiv 0 \pmod{s(n)+d}$. Carrying the digit sum down the recursion makes this test $O(1)$.

4. Only four endings can give a prime

A one-digit prefix is never strong, because $s(h)=h$ forces $h/s(h)=1$, which is not prime. So every usable prefix has $h\ge 10$ and every candidate $p=10h+d$ is at least $100$. A prime that large is divisible by neither $2$ nor $5$, so $d$ is neither even nor $5$: it lies in $\{1,3,7,9\}$. A strong prefix therefore costs at most four primality tests, however long it is.

5. Worked example: the checkpoint $90619$

A two-digit number is right truncatable exactly when it is Harshad, its only proper prefix being a single digit. The $23$ two-digit Harshad numbers, which agree with the published list, ^[4] are

$$\begin{aligned} &10, 12, 18, 20, 21, 24, 27, 30,\\ &36, 40, 42, 45, 48, 50, 54, 60,\\ &63, 70, 72, 80, 81, 84, 90 . \end{aligned}$$

Their quotients $n/s(n)$ are

$$\begin{aligned} &10, 4, 2, 10, 7, 4, 3, 10, 4, 10, 7, 5,\\ &4, 10, 6, 10, 7, 10, 8, 10, 9, 7, 10 \end{aligned}$$

of which exactly seven are prime, so the strong two-digit prefixes are $18, 21, 27, 42, 45, 63, 84$.

Now append the four admissible digits. From $18$ only $181$ is prime, since $183$ and $189$ are multiples of $3$ and $187=11\cdot 17$; $21$ gives $211$, $42$ gives $421$, $63$ gives $631$; $27$ gives two, $271$ and $277$; $45$ gives only $457$, as $451=11\cdot 41$; and $84$ gives none, because $841=29^{2}$, $847=7\cdot 11^{2}$, and $843$ and $849$ are multiples of $3$. So a strong prefix need not produce a prime; conversely the search must grow prefixes that are not strong, since $20$, of quotient $10$, is not strong, yet its child $201$, of quotient $67$, is.

Those seven primes below $1000$ sum to $2449$. One level deeper, the $13$ strong three-digit prefixes

$$\begin{aligned} &201, 207, 209, 247, 402, 407, 423,\\ &481, 603, 801, 803, 804, 846 \end{aligned}$$

yield $18$ primes between $1000$ and $9999$ summing to $88170$. Together $2449+88170=90619$, the value the statement supplies.

6. Why the length must be bounded

Appending the digit $0$ never fails: $s(10n)=s(n)$ by the recurrence above, and $s(n)\mid n$ implies $s(n)\mid 10n$. Every node thus has a child, so the family is infinite, containing $10^{k}$ for every $k$; only the bound on $p$ makes the search finite. Forbid the digit $0$ and the tree dies out after $35$ numbers, the largest $24786$. ^[5]

From $p=10h+d\lt L$ a usable prefix satisfies $h\le H$, where

$$H=\Bigl\lfloor\frac{L-1}{10}\Bigr\rfloor=9999999999999 .$$

The rule used is sharper: stop descending from $n$ once $n\gt\lfloor H/10\rfloor$. Nothing is lost. Then $10n\ge 10\lfloor H/10\rfloor+10\gt H$, so every descendant $h'$ has $h'\ge H+1$; and $H+1\gt (L-1)/10$ gives $10(H+1)\ge L$, so a prime built on such an $h'$ would already reach $L$. The node $n$ is still tested for strength before the descent stops, so its own four candidates are not missed. For $L=10^{14}$ the rule reads $n\gt 999999999999$: prefixes of at most $12$ digits are expanded, those of $13$ digits are leaves.

7. How large the tree is

By decimal length the node counts are

$$\begin{aligned} &9, 23, 54, 94, 175, 271, 451, 646,\\ &1001, 1386, 2019, 2760, 3864 \end{aligned}$$

for lengths $1$ to $13$, so $12753$ prefixes in all, of which $205$ are strong. These grow far more slowly than the tenfold growth of all decimal strings: the last listed ratio is $3864/2760=1.4$, and extending those $3864$ prefixes once more yields $5156$ children, a ratio of $1.33$. Finitely many ratios prove no limiting growth factor, but the mechanism is visible. A child $10n+d$ must be divisible by $s(n)+d$, and that divisor grows as the prefix lengthens — the mean digit sum at length $13$ is $19.14$ — so most of the ten endings fail; only $d=0$, which leaves the digit sum unchanged, is guaranteed to work.

8. Testing primality

Both questions that arise, the primality of $h/s(h)$ and of $10h+d$, concern numbers below $10^{14}$. Trial division would still work, needing divisors only up to $10^{7}$, but repeating that for thousands of candidates is wasteful, and a sieve to $10^{14}$ is out of reach, so a strong probable-prime test is used. Writing $n-1=2^{e}u$ with $u$ odd, $n$ passes for a base $a$ when $a^{u}\equiv 1$, or when $a^{2^{r}u}\equiv -1 \pmod n$ for some $0\le r\lt e$. For $n=13$, where $13-1=4\cdot 3$, base $2$ gives $2^{3}\equiv 8$ and then $8^{2}\equiv -1 \pmod{13}$, so $13$ passes. Every prime passes every base, because $a^{n-1}\equiv 1$ by Fermat's little theorem and $1$ has only the square roots $\pm 1$ modulo a prime, so the repeated squarings must meet $1$ or $-1$. A composite can pass single bases, but the smallest odd composite passing all of $2,3,5,7,11,13,17$ is $341550071728321$. ^[6] Everything tested here is smaller, so passing all seven is a proof, not a probability.

How the Code Works

The C++ program ^[7] defaults limit to $10^{14}$ and lets parse_arguments and parse_u64_after_prefix override it through --limit=, with --skip-checkpoints clearing run_checkpoints. Primality follows Section 8: mul_mod widens its arguments to __uint128_t before reducing, since a product of two values near $10^{14}$ reaches $10^{28}$ ^[8]; pow_mod squares repeatedly, multiplies into the result at each set exponent bit and reduces after every product; and is_prime first divides by the twelve primes up to $37$, discarding most composites, before testing the seven bases.

The search of Sections 3 to 6 is dfs_harshad(n, digit_sum, harshad_limit, prime_limit), returning the subtree sum below $n$. For $n\ge 10$ it applies the strength test and, on success, adds those of $10n+1, 10n+3, 10n+7, 10n+9$ below prime_limit that are prime; then it applies the stopping rule, otherwise looping over $d=0,\dots ,9$ and recursing into every child that passes the divisibility test with digit sum digit_sum + d. solve sets harshad_limit to $\lfloor (L-1)/10\rfloor$ and seeds the nine roots. Two tests there are redundant: next_sum != 0 is always true because every root digit sum is positive, and n % digit_sum == 0 is always true because every visited node is Harshad by construction.

brute_small is the independent check, its lambdas digit_sum and is_right_trunc_harshad recomputing everything for each integer below the limit; run_checkpoints compares solve with $90619$ at $10^{4}$ and with brute_small at $10^{6}$, and main exits nonzero on failure. It shares is_prime with the fast search, so their agreement validates the enumeration but could not expose a fault common to both. Python implements the same recurrence in arbitrary-precision integers, accumulating into a nonlocal total instead of returning subtree sums and shipping its own pow_mod rather than the built-in pow; its sys.setrecursionlimit(100000) is generous, as Section 6 bounds the depth by $13$. Both Python and Java hard-code the limit and run no checkpoints; the options and the comparisons are C++ only. Java mirrors the C++ structure in dfsHarshad and leaves the modular arithmetic to BigInteger.modPow and BigInteger.multiply, so it needs no wide-integer helper. ^[8] A signed long suffices for its accumulator: at most $86$ primes below $10^{14}$ are added, a sum under $10^{16}$.

Complexity and Verification

The work is proportional to the size of the tree, not to $L$. Let $T(L)$ be the number of prefixes the search visits; each tries at most ten children and makes at most five calls to is_prime, so the cost is $O(T(L)\log L)$ operations, counting one step per operation on a bounded integer. For $L=10^{D}$ the visited prefixes are exactly the right truncatable Harshad numbers of at most $D-1$ digits, counted in Section 7; for other bounds a few prefixes above harshad_limit are visited too, harmlessly, since the cand comparison discards them. At $L=10^{14}$ an instrumented copy counts $12753$ visits, $88890$ divisibility tests, $13564$ calls to is_prime and $86$ qualifying primes. Only $719$ of those calls get past the small-prime divisions, together performing $2393$ base exponentiations and $100375$ modular multiplications, a call stopping at the first base that proves compositeness. Auxiliary memory is $O(\log L)$ words: one frame per prefix digit, at most $13$, with no sieve and no table of visited nodes, only fixed-size lists of trial primes, bases and digits. On one machine the search takes about $3$ milliseconds, while a default run takes about $0.11$ seconds, dominated by the checkpoint. There brute_small touches every integer below its bound $B$ and costs about $B\log B$ word operations, since a fixed fraction of integers survives trial division — $148724$ of the $10^{6}$ at the checkpoint — and each survivor then costs $\Theta(\log B)$ steps.

For this article the prefix search and a brute-force scan agree at $10^{4}$, both giving the statement's $90619$, at $10^{5}$, giving $388207$, and at $10^{6}$, giving $1188721$; the search alone gives $10057005$ at $10^{7}$ and $130459097$ at $10^{8}$, the compiled C++ reproduces those and $2449$ at $10^{3}$, and all three programs agree at the full bound. Section 5 was reproduced by code, and the enumeration with $0$ forbidden gave exactly $35$ numbers with largest element $24786$, matching the published zero-free sequence. ^[5] Two caveats tie all of this to the requested bound: the seven bases are proven sufficient only below $341550071728321$, so pushing --limit past roughly $3.4\cdot 10^{14}$ would silently void that guarantee, and the decimal parser has no overflow check, so an argument beyond the unsigned 64-bit range wraps instead of being refused.

Footnotes and References

The references supply the background facts; the tree of prefixes, the pruning rule, the four endings and the operation counts are derived here.

Project Euler 387 — Harshad Numbers. The official statement defines Harshad or Niven numbers, right truncatable and strong Harshad numbers, and the primes whose last-digit truncation is such a number, illustrating them with 201, 20, 2, with 201/3 = 67 and the prime 2011, and gives 90619 as the sum below 10000.
Wikipedia — Harshad number. The source defines a Harshad number by divisibility by its digit sum in the given base and illustrates it with the base-ten examples 18 and 19 used in Section 1.
Wikipedia — Depth-first search. The traversal goes as deep as possible along each branch before backtracking, keeping only the current search path on a stack, so a depth-limited search needs space proportional to the depth limit.
OEIS A005349 — Niven (or Harshad) numbers. The sequence lists the numbers divisible by their digit sum: 1 to 9, 10, 12, 18, 20, 21, 24, 27, 30 and so on; its 23 two-digit entries are those of Section 5.
OEIS A097569 — right-truncatable Harshad numbers (zeros not permitted). The zero-free analogue of the tree of Section 3; the entry records that it is finite, with 24786 as its last term.
Wikipedia — Miller–Rabin primality test defines the strong probable-prime condition of Section 8 and lists the verified small-base bounds, among them that bases 2, 3, 5, 7, 11, 13, 17 suffice below 341550071728321 — the seventh term of OEIS A014233, the smallest odd composite that is a strong pseudoprime to the first n prime bases, citing G. Jaeschke (1993).
C++ — Euler387.cpp. The linked immutable revision contains is_prime, dfs_harshad, solve, brute_small and run_checkpoints; the counts above come from an instrumented re-implementation.
GCC manual — 128-bit integers documents unsigned __int128 as an extension where the target has a wide enough integer mode; it holds the exact product inside mul_mod. Java — BigInteger provides arbitrary-precision integers whose modPow returns a non-negative residue below the modulus.

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());
    }
}