Problem 171: Finding Numbers for Which the Sum of the Squares of the Digits Is a Perfect Square

Project Euler Problem 171 Solution

EulerSolve provides an optimized solution for Project Euler Problem 171, Finding Numbers for Which the Sum of the Squares of the Digits Is a Perfect Square, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $0 \le n < 10^{20}$, write its 20-digit decimal expansion $d_{19}d_{18}\dots d_0$, allowing leading zeros. Define $$f(n)=\sum_{i=0}^{19} d_i^2.$$ The task is to add all numbers $n$ for which $f(n)$ is a perfect square and keep only the last nine digits of that total. Direct enumeration is impossible because there are $10^{20}$ candidates. The useful observation is that the condition depends only on the accumulated square-digit-sum, while the requested output is itself a sum over the actual decimal values. Mathematical Approach Let $L=20$ and $M=10^9$. Since each digit contributes at most $9^2=81$, the largest possible value of $f(n)$ is $20\cdot 81=1620$. That means the entire state space for the square-digit-sum has only 1621 possibilities, which is small enough for a complete dynamic program. Representing the Search Space For each $\ell\in\{0,1,\dots,20\}$ and each $t\in\{0,1,\dots,1620\}$, let $A_{\ell,t}$ be the set of integers $x$ with $0 \le x < 10^\ell$ whose $\ell$-digit expansion, padded with leading zeros if necessary, has square-digit-sum $t$. Thinking in fixed length is exactly what the problem needs: every number below $10^{20}$ is represented once as a 20-digit string, and the all-zero string contributes value 0 anyway....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $f(n)$ be the sum of the squares of the decimal digits of $n$; the statement gives $f(3)=9$, $f(25)=29$ and $f(442)=36$. We must add up all integers $n$ with $0\lt n\lt10^{20}$ for which $f(n)$ is a perfect square and give the last nine digits of this total, that is, its remainder modulo $10^9$. ^[1] For example, $f(34)=9+16=25=5^2$, so $34$ is one of the numbers to be added, while $25$ is not, because $29$ is not a square. The qualifying numbers are the terms of the OEIS sequence A175396 below $10^{20}$, apart from its first term $0$. ^[2]

The whole solution in one sentence: build the zero-padded $20$-digit strings one digit at a time, group them by their digit-square sum, keep for each group the number of strings and the total of their values modulo $10^9$, and finally add the value totals of the groups whose digit-square sum is a positive square.

Mathematical Approach

Testing the numbers one by one is hopeless: there are $10^{20}-1$ of them, and at a billion tests per second the search would take more than $3000$ years. What saves us is that $f(n)$ takes very few different values.

1. Small cases by hand

A one-digit number $n$ has $f(n)=n^2$, which is always a square, so all of $1,2,\dots,9$ qualify: $9$ numbers with sum $45$. A two-digit number $n=10a+b$ with $1\le a\le9$ and $0\le b\le9$ qualifies when $a^2+b^2$ is a square.

If $b=0$, then $f(n)=a^2$ is a square, so $10,20,\dots,90$ qualify, with sum $450$.
If $b\ge1$, we need a Pythagorean triple whose two legs are digits. Checking the $45$ pairs $1\le a\le b\le9$ finds only $3^2+4^2=5^2$ and $6^2+8^2=10^2$, that is, the triple $(3,4,5)$ and its double $(6,8,10)$. ^[3] So $34$, $43$, $68$ and $86$ qualify, with sum $231$.

Below $100$ there are therefore $9+13=22$ qualifying numbers, with sum $45+450+231=726$. A computer search continues the list: below $10^3$ there are $75$ of them with sum $28083$, and below $10^4$ there are $407$ with sum $1719828$.

2. Only 1621 possible sums

A number below $10^{20}$ has at most $20$ digits, and each digit adds at most $9^2=81$, so $0\le f(n)\le20\cdot81=1620$. These numbers therefore have at most $1621$ different values of $f$, and the perfect squares among them are $1^2,2^2,\dots,40^2$, because $40^2=1600\le1620\lt1681=41^2$. Whether $n$ qualifies depends only on $f(n)$, so for each of the $1621$ sums it is enough to know how many numbers have that sum and what they add up to.

3. Strings of fixed length

Write every integer $0\le n\lt10^{L}$ with exactly $L$ digits, padding with leading zeros; for $L=20$ the number $34$ becomes the string of $18$ zeros followed by $3$ and $4$. A string $x=d_{L-1}\dots d_1d_0$ has the value $v(x)=\sum_i d_i\,10^i$, because each digit counts with the power of ten of its place, ^[4] and every integer $0\le n\lt10^L$ is the value of exactly one such string. We write $\sigma(x)=\sum_i d_i^2$ for its digit-square sum. A leading zero adds $0$ to the value and $0^2=0$ to $\sigma$, so padding changes neither $n$ nor $f(n)$. The only string with $\sigma=0$ is the string of zeros, whose value is $0$. Hence the total we want is the sum of $v(x)$ over all strings $x$ of length $20$ for which $\sigma(x)$ is one of $1^2,2^2,\dots,40^2$.

4. Two tables: counts and values

For $\ell\ge0$ and $t\ge0$, let $C_\ell(t)$ be the number of strings $x$ of length $\ell$ with $\sigma(x)=t$, and let $V_\ell(t)$ be the sum of their values. For example, the strings of length $2$ with $\sigma=25$ are $05$, $50$, $34$ and $43$, so $C_2(25)=4$ and $V_2(25)=5+50+34+43=132$. By Section 3, the total asked for is

$$S=\sum_{r=1}^{40}V_{20}\bigl(r^2\bigr).$$

For length $2$ the formula reproduces the sum of Section 1: each square $t=d^2$ with $d\ne5$ is reached only by $0d$ and $d0$, with total $11d$; $t=25$ gives $132$, $t=100$ gives $68+86=154$, and $121$ and $144$ are not reached. Since $1+2+3+4+6+7+8+9=40$, the sum is $11\cdot40+132+154=726$.

5. Adding a leading digit

Every string of length $\ell\ge1$ arises in exactly one way by writing a digit $d$ in front of a string $x$ of length $\ell-1$. Call the result $dx$. The new digit stands at the place with value $10^{\ell-1}$, so

$$\sigma(dx)=d^2+\sigma(x),\qquad v(dx)=d\cdot10^{\ell-1}+v(x).$$

Now fix $t$ and sort the strings of length $\ell$ with $\sigma=t$ by their first digit $d$. Those that begin with $d$ are exactly the strings $dx$ with $\sigma(x)=t-d^2$. There are $C_{\ell-1}(t-d^2)$ of them, and their values add up to $V_{\ell-1}(t-d^2)$ plus $d\cdot10^{\ell-1}$ once for each of them. Adding over the digits gives, for $\ell\ge1$,

$$C_\ell(t)=\sum_{d=0}^{9}C_{\ell-1}\bigl(t-d^2\bigr),$$

$$\begin{aligned}V_\ell(t)&=\sum_{d=0}^{9}V_{\ell-1}\bigl(t-d^2\bigr)\\&\quad+10^{\ell-1}\sum_{d=0}^{9}d\,C_{\ell-1}\bigl(t-d^2\bigr),\end{aligned}$$

where terms with $t-d^2\lt0$ are left out. The second line is the reason for keeping the count table: the place value of the new digit is added once for every string it is written in front of, and only $C_{\ell-1}$ knows how many such strings there are. The sum $\sum_d d\,C_{\ell-1}(t-d^2)$ in that line is the total of the first digits of all strings counted by $C_\ell(t)$. The start is the empty string, of length $0$, value $0$ and sum $0$: $C_0(0)=1$, $V_0(0)=0$, and every other entry of length $0$ is $0$. The argument holds for every $\ell\ge1$, so by induction on $\ell$ the two recurrences produce the true tables for all lengths. Filling a table from the entries of smaller, already solved cases in this way is called dynamic programming. ^[5]

6. Worked example: three digits with sum 9

At length $1$ the tables are simple: $C_1(d^2)=1$ and $V_1(d^2)=d$ for every digit $d$, and all other entries are $0$. One more step gives the length-2 entries we need: $C_2(9)=2$ and $V_2(9)=33$ from $03$ and $30$; $C_2(8)=1$ and $V_2(8)=22$ from $22$; $C_2(5)=2$ and $V_2(5)=33$ from $12$ and $21$; $C_2(0)=1$ and $V_2(0)=0$ from $00$. For length $3$ and $t=9$, only the digits $d=0,1,2,3$ have $d^2\le9$, and the recurrences give

$$C_3(9)=C_2(9)+C_2(8)+C_2(5)+C_2(0)=6,$$

$$\begin{aligned}V_3(9)&=(33+22+33+0)\\&\quad+100\,(0\cdot2+1\cdot1+2\cdot2+3\cdot1)\\&=88+800=888.\end{aligned}$$

Listing the strings confirms both numbers. Since $9$ is either $3^2$ alone or $1^2+2^2+2^2$, the strings are $003$, $030$, $300$, $122$, $212$ and $221$, six of them, and $3+30+300+122+212+221=888$. The leading digit $2$ stands in front of both strings counted by $C_2(5)=2$, giving $212$ and $221$, so $200$ is added twice: this is the factor $C_{\ell-1}(t-d^2)$ at work.

7. Only the last nine digits

The exact tables grow quickly. At length $20$ the largest count is $C_{20}(562)\approx3.3\cdot10^{17}$, but the largest value total is about $1.7\cdot10^{37}$, and the exact $S$ has $39$ digits, far beyond the range of $64$-bit integers, which ends below $2^{64}\approx1.8\cdot10^{19}$. Only $S\bmod10^9$ is wanted, and remainders are compatible with the two operations used: if $a\equiv a'$ and $b\equiv b'\pmod m$, then $a+b\equiv a'+b'$ and $ab\equiv a'b'\pmod m$. ^[6] The recurrences use only additions and multiplications, so we may store every entry of both tables, and every power $10^{\ell-1}$, as its remainder modulo $M=10^9$ and reduce after every operation. By induction on $\ell$, each stored entry is then congruent to the true one modulo $M$, and so is the final sum of the $40$ entries; as that sum is reduced to a number between $0$ and $M-1$, it is exactly $S\bmod M$. The programs reduce the count table as well: it enters the value recurrence only through products, so its residues suffice, and small residues keep these products small. Every stored number is then below $10^9\lt2^{30}$, and the largest intermediate result, a reduced product $(d\cdot10^{\ell-1}\bmod M)\cdot C$ plus a stored value, is at most $(M-1)^2+(M-1)\lt10^{18}\lt2^{63}$.

8. A check from symmetry

The value table can also be obtained from the count table alone. Exchanging two positions $i$ and $j$ of a string does not change $\sigma$, so it maps the strings of length $L$ with $\sigma=t$ one-to-one onto themselves, and the digit at position $i$ of a string becomes the digit at position $j$ of its image. Hence the digits at position $i$ of these strings have the same total for every $i$; call it $D_L(t)$. Adding $d_i\,10^i$ over all positions and all strings gives

$$V_L(t)=R_L\,D_L(t),\qquad R_L=\frac{10^L-1}{9},$$

where $R_L$ is the number written with $L$ ones. At the first position $D_L(t)=\sum_d d\,C_{L-1}(t-d^2)$, the second sum of Section 5. In Section 6 the last digits of the six strings add up to $3+0+0+2+2+1=8$, and $111\cdot8=888$; for $t=25$ at length $2$ the last digits $5,0,4,3$ add up to $12$, and $11\cdot12=132$. The published programs do not use this formula, but it gives an independent check of the value table.

9. The count table as a polynomial power

Let $P(x)=\sum_{d=0}^{9}x^{d^2}$, that is, $P(x)=1+x+x^4+\dots+x^{81}$. For example, $x^{25}$ arises in $P(x)^2$ in four ways, as $1\cdot x^{25}$, $x^{25}\cdot1$, $x^9\cdot x^{16}$ and $x^{16}\cdot x^9$, one for each of the strings $05$, $50$, $34$ and $43$ of Section 4. Expanding $P(x)^L$ means choosing one term $x^{d^2}$ from each of the $L$ factors, which is the same as choosing a string of $L$ digits, and the product of the chosen terms is $x^{\sigma}$, with $\sigma$ the digit-square sum of that string. So $C_L(t)$ is the coefficient of $x^t$ in $P(x)^L$: this power is the generating function of the count table, and the recurrence for $C_\ell$ is the multiplication of $P(x)^{\ell-1}$ by $P(x)$, written coefficient by coefficient as a discrete convolution. ^[7] By Section 8, the value table has the generating function $R_L\,Q(x)\,P(x)^{L-1}$, where $Q(x)=\sum_d d\,x^{d^2}$.

How the Code Works

All three programs carry out Sections 4 to 7 directly. The C++ version ^[8] starts in main, which hands the command line to parse_arguments. This accepts --length= and --modulo= followed by decimal digits, read by parse_int_after_prefix, with the defaults $20$ and $10^9$, and the switch --skip-checkpoints; any other argument, or a length or modulus below $1$, ends the program with exit code $1$. Values too large for an int are not detected.

solve(length, modulo) keeps the tables of Section 4 in the vectors count and sum of unsigned 64-bit integers (u64), with one row for each length from $0$ to length and one column for each sum from $0$ to 81 * length. It sets count[0][0] = 1 and stores $10^i\bmod M$ in pow10[i], where $M$ is modulo. Three loops follow: len runs from $1$ to length, digit from $0$ to $9$, and s from the square sq of the digit up to 81 * length. Each step reads cnt_prev and sum_prev from row len - 1 at column s - sq, adds cnt_prev to count[len][s], and adds sum_prev plus the place-value term pow10[len - 1] * digit % modulo * cnt_prev to sum[len][s], reducing modulo $M$ after every addition and product as in Section 7. A step whose two inputs are both $0$ is skipped, which changes nothing, because it would add $0$. Row len is written only from row len - 1, so the order of the two inner loops does not matter. At the end solve adds sum[length][r * r] for $r=1,2,\dots$ as long as r * r <= 81 * length, which is the formula for $S$ in Section 4.

Before the main computation, unless --skip-checkpoints is given, run_checkpoints compares solve(5, 1000000000) and solve(6, 1000000000) with brute_small. That function goes through every $x$ from $1$ to $10^5-1$ or $10^6-1$, extracts the digits with % 10 and / 10, truncates the floating-point square root of the digit-square sum to an integer r and accepts $x$ when r * r == sq_sum. The equality rules out accepting a nonsquare, and every square is accepted, because the square root of a perfect square of at most $6\cdot81=486$ is computed exactly. If a comparison fails, the program reports it on standard error and exits with code $2$; otherwise it prints solve(length, modulo). The options also give small cases: --length=2 prints $726$ and --length=3 prints $28083$, the totals of Section 1.

The Python file defines solve() and calls it at the end; it builds the lists count, sum_val and pow10 and runs the same three loops and the final loop over the squares. The Java program does the same directly in main, with the long arrays count, sumV and pow10. Both fix the length $20$ and the modulus $10^9$ in the code and have no options and no checkpoints. Python integers have unlimited size, so nothing can overflow there; in Java the bound of Section 7 keeps every intermediate result below $10^{18}$, inside the signed 64-bit range.

Complexity and Verification

For each of the $L$ lengths and each digit $d$, the inner loop visits the $81L+1-d^2$ sums from $d^2$ to $81L$. As $0^2+1^2+\dots+9^2=285$, this makes $L(810L-275)$ steps in total, which is $318500$ for $L=20$. Of these, $139660$ perform an update and $178840$ are skipped by the zero test, $164400$ of them because the column read lies above $81(\ell-1)$, the largest sum of a string of length $\ell-1$. The running time is therefore $\Theta(L^2)$. All three programs keep every row of both tables, $2\cdot21\cdot1621=68082$ entries, so the memory also grows as $\Theta(L^2)$; in C++ and Java the entries are $8$-byte integers, $544656$ bytes in total, while Python stores references to integer objects. Since row $\ell$ needs only row $\ell-1$, two rows per table would suffice. In our tests each of the three programs finished in well under a second, the C++ program including its checkpoints.

The C++ checkpoints compare the dynamic program with two brute-force runs, over the numbers below $10^5$ and below $10^6$, $1099998$ tests in total. For this article we also checked the following.

A separate brute force in C, which runs through all numbers below $10^{10}$, found the exact count and sum for every $L=1,\dots,10$; the counts are $9$, $22$, $75$, $407$, $3525$, $35178$, $368816$, $3642462$, $34746543$ and $319925193$. Exact tables built with the recurrences give the same counts and sums, and the C++ program run with --length=L prints their remainders modulo $10^9$.
For the lengths $0$ to $4$, the recurrences agree entry by entry with tables made by listing all strings.
At length $20$, the exact count table equals the list of coefficients of $P(x)^{20}$, computed by repeated squaring, and the exact value table satisfies $V_{20}(t)=R_{20}\,D_{20}(t)$ for all $1621$ sums.
The exact $S$ from these tables, reduced modulo $10^9$, equals the number printed by all three programs.

Footnotes and References

The references below supply the background facts used above. The two recurrences, the reduction modulo $10^9$, the symmetry formula and the operation counts are derived explicitly in this article.

Project Euler 171 — Square Sum of the Digital Squares. The official statement defines f(n) as the sum of the squares of the decimal digits of a positive integer n, gives f(3) = 9, f(25) = 29 and f(442) = 36, and asks for the last nine digits of the sum of all n with 0 < n < 10²⁰ for which f(n) is a perfect square.
OEIS A175396 — Numbers whose sum of squares of digits is a square. The sequence begins 0, 1, 2, …, 9, 10, 20, 30, 34, 40, 43, 50, 60, 68, 70, 80, 86, 90, 100, 122, …, and its examples explain that 34 is a term because 3² + 4² = 25 = 5², and 122 because 1² + 2² + 2² = 9 = 3².
Wikipedia — Pythagorean triple. A Pythagorean triple consists of three positive integers a, b and c with a² + b² = c²; (3, 4, 5) is the best-known example, every multiple (ka, kb, kc) of a triple is again a triple, and (6, 8, 10) is named as a triple that is not primitive.
Wikipedia — Positional notation. In a positional numeral system a digit's value is the digit multiplied by the value of its place, and the place values are powers of the base; in base 10, for example, 5305 = 5·10³ + 3·10² + 0·10¹ + 5·10⁰.
Wikipedia — Dynamic programming. Dynamic programming simplifies a complicated problem by breaking it down into simpler sub-problems, stores the solutions of the sub-problems in a table and builds the solutions of bigger sub-problems from them, so that each sub-problem is solved only once.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition and multiplication: if a₁ ≡ b₁ and a₂ ≡ b₂ (mod m), then a₁ + a₂ ≡ b₁ + b₂ and a₁a₂ ≡ b₁b₂ (mod m), so the results of additions and multiplications may be replaced by their remainders at every step.
Wikipedia — Generating function. The ordinary generating function of a sequence aₙ is the power series whose coefficient of xⁿ is aₙ, and multiplying two ordinary generating functions yields the discrete convolution, or Cauchy product, of their sequences.
C++ — Euler171.cpp. The linked immutable C++ revision contains parse_arguments, solve, brute_small and run_checkpoints. It is the implementation record for this article.

Problem 171 source code

C++

#include <array>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int length = 20;
    int modulo = 1000000000;
    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, "--length=", options.length) ||
            parse_int_after_prefix(arg, "--modulo=", options.modulo)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.length >= 1 && options.modulo >= 1;
}

u64 solve(const int length, const int modulo) {
    const int max_sq_sum = length * 81;
    std::vector<std::vector<u64>> count(static_cast<std::size_t>(length + 1),
                                        std::vector<u64>(static_cast<std::size_t>(max_sq_sum + 1), 0));
    std::vector<std::vector<u64>> sum(static_cast<std::size_t>(length + 1),
                                      std::vector<u64>(static_cast<std::size_t>(max_sq_sum + 1), 0));
    std::vector<u64> pow10(static_cast<std::size_t>(length + 1), 1);
    for (int i = 1; i <= length; ++i) {
        pow10[static_cast<std::size_t>(i)] = (pow10[static_cast<std::size_t>(i - 1)] * 10ULL) %
                                             static_cast<u64>(modulo);
    }

    count[0][0] = 1;
    for (int len = 1; len <= length; ++len) {
        for (int digit = 0; digit <= 9; ++digit) {
            const int sq = digit * digit;
            for (int s = sq; s <= max_sq_sum; ++s) {
                const u64 cnt_prev = count[static_cast<std::size_t>(len - 1)][static_cast<std::size_t>(s - sq)];
                const u64 sum_prev = sum[static_cast<std::size_t>(len - 1)][static_cast<std::size_t>(s - sq)];
                if (cnt_prev == 0 && sum_prev == 0) {
                    continue;
                }
                count[static_cast<std::size_t>(len)][static_cast<std::size_t>(s)] =
                    (count[static_cast<std::size_t>(len)][static_cast<std::size_t>(s)] + cnt_prev) %
                    static_cast<u64>(modulo);

                const u64 add = (sum_prev +
                                 (pow10[static_cast<std::size_t>(len - 1)] *
                                  static_cast<u64>(digit) % static_cast<u64>(modulo) * cnt_prev) %
                                     static_cast<u64>(modulo)) %
                                static_cast<u64>(modulo);
                sum[static_cast<std::size_t>(len)][static_cast<std::size_t>(s)] =
                    (sum[static_cast<std::size_t>(len)][static_cast<std::size_t>(s)] + add) %
                    static_cast<u64>(modulo);
            }
        }
    }

    u64 answer = 0;
    for (int r = 1; r * r <= max_sq_sum; ++r) {
        answer = (answer + sum[static_cast<std::size_t>(length)][static_cast<std::size_t>(r * r)]) %
                 static_cast<u64>(modulo);
    }
    return answer;
}

u64 brute_small(const int length, const int modulo) {
    u64 answer = 0;
    int total_values = 1;
    for (int i = 0; i < length; ++i) {
        total_values *= 10;
    }

    for (int x = 1; x < total_values; ++x) {
        int y = x;
        int sq_sum = 0;
        for (int i = 0; i < length; ++i) {
            const int d = y % 10;
            sq_sum += d * d;
            y /= 10;
        }
        const int r = static_cast<int>(std::sqrt(static_cast<double>(sq_sum)));
        if (r * r == sq_sum) {
            answer = (answer + static_cast<u64>(x)) % static_cast<u64>(modulo);
        }
    }
    return answer;
}

bool run_checkpoints() {
    if (solve(5, 1000000000) != brute_small(5, 1000000000)) {
        std::cerr << "Checkpoint failed for length 5 brute cross-check" << '\n';
        return false;
    }
    if (solve(6, 1000000000) != brute_small(6, 1000000000)) {
        std::cerr << "Checkpoint failed for length 6 brute cross-check" << '\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.length, options.modulo) << '\n';
    return 0;
}

Python

# Problem 171: Finding numbers for which the sum of squares of digits is a perfect square.
# Sum (mod 10^9) of all 20-digit numbers (including leading zeros) with square digit-square-sum.

def solve():
    length = 20
    MOD = 10**9
    max_sq = length * 81
    # DP: count[len][s] and sum_val[len][s]
    count = [[0]*(max_sq+1) for _ in range(length+1)]
    sum_val = [[0]*(max_sq+1) for _ in range(length+1)]
    pow10 = [1]*(length+1)
    for i in range(1, length+1):
        pow10[i] = pow10[i-1]*10 % MOD
    count[0][0] = 1
    for ln in range(1, length+1):
        for d in range(10):
            sq = d*d
            for s in range(sq, max_sq+1):
                cp = count[ln-1][s-sq]
                sp = sum_val[ln-1][s-sq]
                if cp == 0 and sp == 0: continue
                count[ln][s] = (count[ln][s] + cp) % MOD
                add = (sp + pow10[ln-1]*d%MOD*cp) % MOD
                sum_val[ln][s] = (sum_val[ln][s] + add) % MOD
    ans = 0
    r = 1
    while r*r <= max_sq:
        ans = (ans + sum_val[length][r*r]) % MOD
        r += 1
    print(ans)

solve()

Java

public class Euler171 {
    public static void main(String[] args) {
        int length = 20;
        int MOD = 1000000000;
        int maxSq = length * 81;
        long[][] count = new long[length + 1][maxSq + 1];
        long[][] sumV = new long[length + 1][maxSq + 1];
        long[] pow10 = new long[length + 1];
        pow10[0] = 1;
        for (int i = 1; i <= length; i++)
            pow10[i] = pow10[i - 1] * 10 % MOD;
        count[0][0] = 1;
        for (int ln = 1; ln <= length; ln++)
            for (int d = 0; d <= 9; d++) {
                int sq = d * d;
                for (int s = sq; s <= maxSq; s++) {
                    long cp = count[ln - 1][s - sq], sp = sumV[ln - 1][s - sq];
                    if (cp == 0 && sp == 0)
                        continue;
                    count[ln][s] = (count[ln][s] + cp) % MOD;
                    long add = (sp + pow10[ln - 1] * d % MOD * cp) % MOD;
                    sumV[ln][s] = (sumV[ln][s] + add) % MOD;
                }
            }
        long ans = 0;
        for (int r = 1; r * r <= maxSq; r++)
            ans = (ans + sumV[length][r * r]) % MOD;
        System.out.println(ans);
    }
}