Problem 159: Digital Root Sums of Factorisations

Project Euler Problem 159 Solution

EulerSolve provides an optimized solution for Project Euler Problem 159, Digital Root Sums of Factorisations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every positive integer $n$, let $\operatorname{dr}(n)$ be its digital root, obtained by repeatedly summing decimal digits until only one digit remains. If $$n=f_1f_2\cdots f_k,\qquad f_i \ge 2,$$ define the digital-root sum of that multiplicative factorisation by $$\operatorname{drs}(f_1,\dots,f_k)=\sum_{i=1}^{k}\operatorname{dr}(f_i).$$ The quantity of interest is $$\operatorname{mdrs}(n)=\max \operatorname{drs}(f_1,\dots,f_k),$$ where the maximum runs over all multiplicative factorizations of $n$ into integers at least 2, including the trivial one-part factorisation $n$ itself. The problem asks for $$\sum_{n=2}^{10^6-1}\operatorname{mdrs}(n).$$ A brute-force search over all factorizations of every integer below one million would repeat the same subproblems many times. The implementations avoid that by turning $\operatorname{mdrs}$ into a dynamic program on factor pairs. Mathematical Approach The whole solution rests on one observation: once the best value is known for smaller integers, the best value for a larger composite can be assembled from a split $n=ab$. Digital roots and the score of a factorisation For $n \ge 1$, the digital root is given directly by $$\operatorname{dr}(n)=1+((n-1)\bmod 9).$$ So $\operatorname{dr}(n)$ is always one of $1,2,\dots,9$, and multiples of 9 have digital root 9....

Detailed mathematical approach

Problem Summary

A factorisation of an integer $n\gt1$ writes $n$ as a product of integers greater than $1$. The order of the factors does not matter, $n$ itself counts as a factorisation with a single factor, and multiplication by $1$ is not allowed. For example, $24$ has the $7$ factorisations $2\times2\times2\times3$, $2\times3\times4$, $2\times2\times6$, $4\times6$, $3\times8$, $2\times12$ and $24$. The digital root of a number is found by adding its decimal digits and repeating this until a number below $10$ is left: for $467$ we get $4+6+7=17$ and then $1+7=8$. The Digital Root Sum (DRS) of a factorisation is the sum of the digital roots of its factors, and $\operatorname{mdrs}(n)$ is the largest DRS among all factorisations of $n$. The seven factorisations of $24$ have the sums $9$, $9$, $10$, $10$, $11$, $5$ and $6$, so $\operatorname{mdrs}(24)=11$. We must find the sum of $\operatorname{mdrs}(n)$ over all integers $n$ with $1\lt n\lt 1\,000\,000$. ^[1]

The whole solution in one sentence: a best factorisation of $n$ is either $n$ itself or best factorisations of two factors $a$ and $b$ with $ab=n$ put together, so one sweep over all products $ab$ below one million, in increasing order of $a$, turns a table of digital roots into the table of all values $\operatorname{mdrs}(n)$, whose sum is the answer.

Mathematical Approach

We compute digital roots in one step, look at small cases, prove the recurrence that the programs use, and show that one sweep over a table evaluates it.

1. Digital roots in one step

Adding the digits does not change the remainder after division by $9$. Indeed, if $n$ has the digits $d_i$, so that $n=\sum_i d_i10^i$, then $n-\sum_i d_i=\sum_i d_i(10^i-1)$ is a multiple of $9$, because every $10^i-1=99\ldots9$ is one; this is the old rule of casting out nines. ^[3] For $n\ge10$ the digit sum is smaller than $n$, and for $n\ge1$ it is at least $1$, so the repetition stops at a number from $1$ to $9$ with the same remainder as $n$. Exactly one number in that range has this remainder, so for every $n\ge1$ ^[2]

$$\operatorname{dr}(n)=1+\bigl((n-1)\bmod 9\bigr).$$

For $467$ we get $466=9\cdot51+7$, so $\operatorname{dr}(467)=8$, as in the statement, and a multiple of $9$ gets the digital root $9$, not $0$. This also shows why the statement excludes factors $1$: since $\operatorname{dr}(1)=1$, every extra factor $1$ would add $1$ to the sum, and there would be no largest sum.

2. Small cases

For small $n$ all factorisations can be listed by hand. For $n=2,3,\dots,25$ the values of $\operatorname{mdrs}(n)$ are

$2,3,4,5,6,7,8,9,7,2,8,4,$
$9,8,10,8,11,1,9,10,4,5,11,10.$

They show that no simple rule decides how far to split:

A prime $p$ has only the factorisation $p$, so $\operatorname{mdrs}(p)=\operatorname{dr}(p)$: $\operatorname{mdrs}(11)=2$ and $\operatorname{mdrs}(19)=1$.
Splitting can help a lot: $10=2\times5$ gives $2+5=7$, while $\operatorname{dr}(10)=1$.
Splitting as far as possible is not always best: for $24$ the factorisation into primes gives $9$, but $3\times8$ gives $11$; for $16$ the primes give $8$, but $2\times8$ gives $10$.
Not splitting can be best: $26=2\times13$ gives only $2+4=6$, so $\operatorname{mdrs}(26)=\operatorname{dr}(26)=8$. Below $100$, the composite numbers that are strictly better left whole are $6$, $8$, $9$, $26$, $33$, $44$, $62$ and $87$.
Both effects mix: the four factorisations of $52$ have the sums $7$ for $52$, $10$ for $2\times26$, and $8$ for $4\times13$ and for $2\times2\times13$, so the best one splits off the factor $2$ but keeps the composite factor $26$ whole.

So we need a systematic search, and listing every factorisation would do far more work than necessary. The number of factorisations of $n$ is the multiplicative partition function: ^[4] ^[5] it is $7$ for $24$ and $22711$ for $907200$, and the numbers from $2$ to $999999$ have $43\,088\,487$ factorisations together.

3. The recurrence

Claim. For every $n\ge2$,

$$\operatorname{mdrs}(n)=\max\bigl(\operatorname{dr}(n),\,S(n)\bigr),$$

$$S(n)=\max_{ab=n}\bigl(\operatorname{mdrs}(a)+\operatorname{mdrs}(b)\bigr),$$

where $a$ and $b$ run over the integers $a,b\ge2$ with $ab=n$; for a prime $n$ there is no such pair, $S(n)$ is left out, and $\operatorname{mdrs}(n)=\operatorname{dr}(n)$.

Proof. Every value on the right is the DRS of some factorisation of $n$: $\operatorname{dr}(n)$ belongs to the factorisation $n$ itself, and $\operatorname{mdrs}(a)+\operatorname{mdrs}(b)$ to a best factorisation of $a$ written next to a best factorisation of $b$, whose factors multiply to $ab=n$. So the right side is at most $\operatorname{mdrs}(n)$. Conversely, take a best factorisation $n=f_1f_2\cdots f_k$. If $k=1$, its DRS is $\operatorname{dr}(n)$. If $k\ge2$, put $a=f_1$ and $b=f_2\cdots f_k$, both at least $2$. The DRS is $\operatorname{dr}(f_1)$, which is at most $\operatorname{mdrs}(a)$, plus the DRS of the factorisation $f_2\cdots f_k$ of $b$, which is at most $\operatorname{mdrs}(b)$. So $\operatorname{mdrs}(n)$ is at most the right side. $\square$

This is the optimal substructure that dynamic programming needs: a best solution is assembled from best solutions of smaller problems, which are computed once and stored. ^[7] For $24$ the pairs $2\cdot12$, $3\cdot8$ and $4\cdot6$ give $2+8=10$, $3+8=11$ and $4+6=10$ with the values of Section 2, and $\operatorname{dr}(24)=6$, so again $\operatorname{mdrs}(24)=11$. The On-Line Encyclopedia of Integer Sequences lists $\operatorname{mdrs}(n)$ as A347105, with this recurrence. ^[6]

4. How large the values can get

Every factor contributes at most $9$, and the factor $9$ gives the most for its size: for every $f\ge2$,

$$\operatorname{dr}(f)\le9\log_9 f,$$

with equality only for $f=9$. For $f\ge10$ this holds because $\operatorname{dr}(f)\le9\lt9\log_9f$. For $2\le f\le9$ we have $\operatorname{dr}(f)=f$, and $f\le9\log_9f$ is the same as $9^f\le f^9$, which can be checked one value at a time: $81\le512$ for $f=2$, $729\le19683$ for $f=3$, and so on up to $9^8=43046721\le8^9=134217728$ for $f=8$, with equality for $f=9$. The logarithms of the factors of a factorisation add up to $\log_9n$, so adding the inequality over the factors gives

$$\operatorname{mdrs}(n)\le9\log_9n.$$

Equality needs every factor to be $9$, so it holds exactly for the powers of $9$: $\operatorname{mdrs}(9^k)=9k$. For $24$ the bound is about $13.02$. Below $10^6$ it is $9\log_9 999\,999\approx56.59$; the largest value there is $55$, reached only at $944784=2\cdot8\cdot9^5$ with the sum $2+8+45$, while $9^6=531441$ gets $54$. So every entry of the table is small, and the sum asked for is less than $56\cdot10^6\lt2^{31}$.

5. One sweep over all products

The recurrence needs the values of the factors first, and both factors of a split are at most $n/2$. The programs therefore fill a table $t$ with $t[n]=\operatorname{dr}(n)$ and then make one sweep in rounds $a=2,3,\dots,N-1$, with $N=10^6$. Round $a$ runs through $b=2,3,\dots$ as long as $ab\lt N$ and sets

$$t[ab]\leftarrow\max\bigl(t[ab],\,t[a]+t[b]\bigr).$$

So round $a$ walks through the multiples $2a,3a,\dots$ of $a$, and no number is ever factorised. Every ordered pair $(a,b)$ with $a,b\ge2$ and $ab\lt N$ is visited once, so two different factors $x\ne y$ meet twice: in round $x$ as $(x,y)$ and in round $y$ as $(y,x)$.

6. Why one sweep is enough

Two facts make the sweep correct.

No entry is too large. At every moment $t[n]$ is the DRS of some factorisation of $n$: initially of the factorisation $n$ itself, and after an update of the two factorisations behind $t[a]$ and $t[b]$, written side by side. So always $t[n]\le\operatorname{mdrs}(n)$, and entries only grow.
Every entry is finished in time. When round $n$ begins, $t[n]=\operatorname{mdrs}(n)$.

We prove the second fact by induction on $n$. If the factorisation $n$ itself is best, the initial value is already right. Otherwise, by Section 3, $\operatorname{mdrs}(n)=\operatorname{mdrs}(x)+\operatorname{mdrs}(y)$ for some $xy=n$ with $2\le x\le y$. Both factors are at most $n/2$, so round $y$ comes before round $n$, and in round $y$ the pair $(y,x)$ is visited. At that moment $t[y]=\operatorname{mdrs}(y)$ by induction, since round $y$ has begun, and $t[x]=\operatorname{mdrs}(x)$, since round $x$ began no later; by the first fact these entries keep their values. So round $y$ raises $t[n]$ to at least $\operatorname{mdrs}(n)$, and by the first fact to exactly this value. $\square$

The proof uses the later of the two visits, and the earlier one can indeed come too soon. In round $2$, the pair $(2,25)$ finds $t[25]=7$, the digital root of $25$, because $25=5\cdot5$ is improved to $10$ only in round $5$; so it proposes $2+7=9$ for $50$. Round $5$ proposes $5+t[10]=5+7=12$, and round $25$ proposes $10+2=12$ again, which is $\operatorname{mdrs}(50)$.

7. Worked example: all numbers below 25

Take $N=25$. Rounds $2$ to $12$ visit $37$ pairs; from round $13$ on, $2a\ge25$ and nothing happens. Round $2$ raises $t[10]$ from $1$ to $7$, $t[12]$ from $3$ to $8$, $t[14]$ from $5$ to $9$, $t[16]$ from $7$ to $10$, $t[18]$ from $9$ to $11$, $t[20]$ from $2$ to $9$ and $t[24]$ from $6$ to $10$. Two of these already use entries raised earlier in the same round, which is harmless by the first fact of Section 6: $t[20]=2+t[10]=2+7$ and $t[24]=2+t[12]=2+8$. Round $3$ raises $t[15]$ from $6$ to $8$, $t[21]$ from $3$ to $10$, and $t[24]$ from $10$ to $11$ through $3+t[8]=3+8$. Rounds $4$ to $12$ change nothing. The final table for $n=2,\dots,24$ is the list of Section 2 without its last entry, and its sum is $151$. For $N=100$, $1000$ and $10\,000$ the sums are $856$, $10901$ and $124491$.

8. Other ways to organise the search

In the proof of Section 3, the first factor $f_1$ of a best factorisation enters only through $\operatorname{dr}(f_1)$. So the recurrence stays true if $\operatorname{mdrs}(a)$ is replaced by $\operatorname{dr}(a)$, with $a$ running over all divisors $1\lt a\lt n$. This version peels off one factor at a time and gives the same table.
Starting the inner loop at $b=a$ visits every unordered pair only once, in its earlier round, and nearly halves the sweep for $N=10^6$, from $11\,969\,988$ to $5\,985\,493$ steps; the $998$ squares $a\cdot a$ are visited once either way. A computation shows that it gives the same table for every $n\lt10^6$, but the proof of Section 6, which needs the later visit, no longer applies; the published programs keep the full loop.

How the Code Works

The C++ program ^[9] starts in main with the settings of Options: limit $=10^6$ and run_checkpoints switched on. parse_arguments accepts --skip-checkpoints and --limit= followed by digits, which parse_int_after_prefix reads; any other argument is reported as unknown, and both that and a limit below $2$ end the program with exit code $1$. Next, unless --skip-checkpoints was given, run_checkpoints builds the table for the limit $100$ and checks the statement's value $\operatorname{mdrs}(24)=11$; if this fails, the exit code is $2$. Then main prints solve(limit).

solve calls compute_mdrs, which implements Sections 5 and 6. It creates a std::vector<int> of limit zeros, stores digital_root(n), the formula of Section 1, for every $n\ge1$, and runs the double loop, in which b is counted up together with p += a, so that the partner of a needs no division. Entry $1$ is filled but never read, because both factors are at least $2$. solve then adds the entries $2$ to limit - 1 in a std::int64_t, so for a limit $L$ the program computes the sum over $1\lt n\lt L$; with --limit=25 it prints the $151$ of Section 7.

The Python program fixes the limit at $10^6$ inside solve(), defines digital_root there, and runs the same sweep on a list with a while loop, finding the partner as b = p // a; it prints sum(mdrs[2:]). The Java program does the same with an int[] table, the partner p / a and a long sum. Neither takes arguments. For the limit $10^6$, Section 4 shows that even a signed 32-bit sum would fit; Python's integers are unbounded anyway. For large arguments the C++ program has two overflows of int, which are undefined behaviour: parse_int_after_prefix does not check the size of the number it reads, so --limit=3000000000 already overflows while it is read, and a * 2 overflows for limits above $2^{30}$, where the table alone needs 4 GiB.

Complexity and Verification

Filling the table and adding it up take $O(N)$ steps. Round $a$ visits $\lfloor(N-1)/a\rfloor-1$ pairs, so the sweep makes

$$\begin{aligned}&\sum_{a=2}^{N-1}\Bigl(\Bigl\lfloor\frac{N-1}{a}\Bigr\rfloor-1\Bigr)\\&\qquad=D(N-1)-2N+3\end{aligned}$$

steps, where $D(x)=\sum_{k\le x}\lfloor x/k\rfloor$ is the divisor summatory function, the number of ordered pairs of positive integers with product at most $x$. For $x=N-1$, removing the $N-1$ pairs with first factor $1$ and the $N-1$ pairs with second factor $1$, and putting back the pair $(1,1)$, which was removed twice, leaves the pairs that the sweep visits. Dirichlet proved that $D(x)=x\ln x+(2\gamma-1)x+O(\sqrt x)$, where $\gamma\approx0.5772$ is Euler's constant, ^[8] so the sweep makes $N\ln N+(2\gamma-3)N+O(\sqrt N)$ steps, which is $O(N\log N)$. For $N=10^6$ we have $D(999999)=13\,969\,985$, and the sweep makes exactly $11\,969\,988$ steps; the main terms of the formula give about $11\,969\,942$. The memory is $O(N)$, one table of $N$ integers: 4 MB in Java and in C++ with a four-byte int, and on 64-bit CPython about 8 MB of references for the list, plus almost 8 MB more while mdrs[2:] is copied for the sum. In tests on a desktop computer the C++ and Java programs finished within a few hundredths of a second and the Python program in about half a second.

The C++ program checks $\operatorname{mdrs}(24)=11$ before the main run; the Python and Java programs run no check. For this article all three programs were run and print the same number. A separate program lists all factorisations of every $n$ from $2$ to $999999$, takes the largest DRS for each $n$ and agrees with the table of the sweep for every $n$; it needed under two seconds, against a few hundredths for the sweep. The same search reproduces the seven factorisations of $24$ with the sums of the statement, the values $\operatorname{mdrs}(2),\dots,\operatorname{mdrs}(40)$ listed in A347105, ^[6] and the sums $151$, $856$, $10901$ and $124491$ of Section 7, which the C++ program prints with --limit=25, --limit=100, --limit=1000 and --limit=10000. Further computations confirm the counts of Sections 2 and 4 and that the two variants of Section 8 give the same table below $10^6$. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The recurrence, the bound of Section 4, the correctness proof of the sweep and the operation counts are derived explicitly in this article.

Project Euler 159 — Digital Root Sums of Factorisations. The official statement lists the seven factorisations of 24 that do not use multiplication by one, defines the digital root by repeated digit sums, with the digital root of 467 equal to 8, defines the Digital Root Sum of a factorisation and mdrs(n), gives mdrs(24) = 11, and asks for the sum of mdrs(n) over 1 < n < 1 000 000.
Wikipedia — Digital root. For a positive integer n and a base b, the digital root is 1 + ((n − 1) mod (b − 1)), and it is congruent to n modulo b − 1 because b ≡ 1 (mod b − 1). So it is the number from 1 to b − 1 with the same remainder as n; in base 10 it equals the remainder upon division by 9, except that a remainder 0 gives the digital root 9.
Wikipedia — Casting out nines. Replacing a number by the sum of its digits leaves the same remainder after division by 9, because the difference is made of multiples of the numbers 10ⁱ − 1, which are divisible by 9; repeating the digit sum until a single digit is left gives the digital root.
Wikipedia — Multiplicative partition. A multiplicative partition, or unordered factorisation, of n writes n as a product of integers greater than 1, with products that differ only in the order of the factors counted once and n itself counted as one of them; for example, 20 has the four multiplicative partitions 2 × 2 × 5, 2 × 10, 4 × 5 and 20.
OEIS A001055 — The multiplicative partition function. The number of ways of factoring n with all factors greater than 1; its terms include a(20) = 4 and a(24) = 7, the seven factorisations of the statement.
OEIS A347105. The greatest sum of the digital roots of the factorisations of n, with a(24) = 11 worked out from the seven factorisations of 24, and the formula a(n) = max(DR(n), a(b) + a(c)) over all factorisations n = b·c with 2 ≤ b ≤ √n, where DR(x) = 1 + ((x − 1) mod 9). The entry cites Project Euler Problem 159.
Wikipedia — Dynamic programming. The method, developed by Richard Bellman in the 1950s, applies to problems with optimal substructure, where an optimal solution is obtained by combining optimal solutions of subproblems; storing the solutions of the subproblems in a table means that each of them is solved only once.
Wikipedia — Divisor summatory function. D(x) is the sum of the divisor counts d(n) for n ≤ x, which equals the sum of ⌊x/k⌋ for k ≤ x, and D(x) = x log x + x(2γ − 1) + O(√x), first established by Dirichlet in 1849 with the hyperbola method.
C++ — Euler159.cpp. The linked immutable revision, with compute_mdrs, solve and run_checkpoints, is the implementation record for this article.

Problem 159 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;

struct Options {
    int limit = 1000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 2;
}

int digital_root(const int n) {
    return (n == 0) ? 0 : 1 + (n - 1) % 9;
}

std::vector<int> compute_mdrs(const int limit) {
    std::vector<int> mdrs(static_cast<std::size_t>(limit), 0);
    for (int n = 1; n < limit; ++n) {
        mdrs[static_cast<std::size_t>(n)] = digital_root(n);
    }

    for (int a = 2; a < limit; ++a) {
        for (int b = 2, p = a * 2; p < limit; ++b, p += a) {
            const int cand = mdrs[static_cast<std::size_t>(a)] + mdrs[static_cast<std::size_t>(b)];
            if (cand > mdrs[static_cast<std::size_t>(p)]) {
                mdrs[static_cast<std::size_t>(p)] = cand;
            }
        }
    }

    return mdrs;
}

i64 solve(const int limit) {
    const std::vector<int> mdrs = compute_mdrs(limit);

    i64 sum = 0;
    for (int n = 2; n < limit; ++n) {
        sum += mdrs[static_cast<std::size_t>(n)];
    }
    return sum;
}

bool run_checkpoints() {
    const std::vector<int> mdrs = compute_mdrs(100);
    if (mdrs[24] != 11) {
        std::cerr << "Checkpoint failed for mdrs(24)" << '\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

# Problem 159: Digital root clocks / Maximum Digital Root Sum (MDRS)

def solve():
    limit = 1000000
    
    def digital_root(n):
        if n == 0: return 0
        return 1 + (n - 1) % 9
    
    mdrs = [0] * limit
    for n in range(1, limit):
        mdrs[n] = digital_root(n)
    
    for a in range(2, limit):
        p = a * 2
        while p < limit:
            b = p // a
            cand = mdrs[a] + mdrs[b]
            if cand > mdrs[p]:
                mdrs[p] = cand
            p += a
    
    print(sum(mdrs[2:]))

solve()

Java

public class Euler159 {
    static int digitalRoot(int n) {
        return n == 0 ? 0 : 1 + (n - 1) % 9;
    }

    public static void main(String[] args) {
        int limit = 1000000;
        int[] mdrs = new int[limit];
        for (int n = 1; n < limit; n++)
            mdrs[n] = digitalRoot(n);
        for (int a = 2; a < limit; a++)
            for (int p = a * 2; p < limit; p += a) {
                int cand = mdrs[a] + mdrs[p / a];
                if (cand > mdrs[p])
                    mdrs[p] = cand;
            }
        long sum = 0;
        for (int n = 2; n < limit; n++)
            sum += mdrs[n];
        System.out.println(sum);
    }
}