Square Sum of the Digital Squares

For a positive integer $n$, let $f(n)$ be the sum of the squares of the digits (in base $10$) of $n$, e.g.
\begin{align}
f(3) &= 3^2 = 9,\\
f(25) &= 2^2 + 5^2 = 4 + 25 = 29,\\
f(442) &= 4^2 + 4^2 + 2^2 = 16 + 16 + 4 = 36\\
\end{align}
Find the last nine digits of the sum of all $n$, $0 \lt n \lt 10^{20}$, such that $f(n)$ is a perfect square.

Analysis of this problem hints that we will need to use the generating functions that can hint on the number of ways a number can be formed. Remember that the objective here is to count all numbers from 1 to 10^20 whereby their digit square sums are perfect squares.

The trick is to apply polynomial generating functions let’s say for instance (x^1 + x^4 + x^9 + x^16 + x^25 + x^36 + x^49 + x^64 + x^81)^20

This lengthy polynomial tells you the number of ways you can form a certain number by using 20 digits (where 0 is not counted on the leading digit). Why? because every term’s exponent in the expansion tells you the sum of squares of some 20 digits, and its coefficient tells you the number of such 20-digit numbers (note that it’s less by one dimension because 0 is not allowed as leading number). The exponents run from 20 to 1620, every natural number possible with the square sum of 20 digits.

Ultimately, our mission is to seek terms in the polynomial expansion whose exponents are perfect squares. However, generating such a polynomial manually and subsequently finding its terms with perfect square exponents is practically impossible.

For that reason, we have to be clever. The exponents with perfect squares are sparse. By the Cauchy’s Coefficient Formula, we know that a term in a polynomial with exponent m is generated by the 20th derivative evaluated at 0 and divided by m!. However, this formula applies to exponents less than 20.

Therefore, the construction of the polynomial geared towards filtering perfect square terms will require careful construction of the polynomial to save on unnecessary calculations.

Here’s a pseudo code:
“Sum of n” = 0
Max_A = sqrt(1620), which is the maximum of the all square sums of digits.
For A from 1 to Max_A do
Generating f(x) = (x+ x^4 + x^9 + x^16 + x^25 + x^36 + x^49 + x^64 + x^81)^20
Do 20 times derivation on f(x)
Evaluate f^(20) (0), and divide it by (A^2)!
Multiply the result by (10^20 – 1)//9 (which is the sum of whole numbers up to 10^19)
Add the result to “sum of n”

The terms generated by this procedure are for 20 digit numbers (or 19 digits if you consider the leading digit as non-zero). To account for the rest of the digits, achieve this by dynamic programming. For every digit D and every possible square sum A (where D is from 1 to 19, and A is from 1 to Max_A), Generate g = (x + x^2 + x^3 + … + x^D-1)^A * x^D. Get g’s term with exponent being a perfect square, add all these terms up, multiply by D!, and add it to “sum of n”.

Be cautious when doing this by programming. The coefficients are large and therefore you should use a Big-integer library (or write one yourself). The calculation is still hefty, but it’s absolutely possible to write and run a program for this problem to get the answer.

We see there are only 40 numbers from 1^2 to 40^2. The main computational cost comes from differentiating polynomials. If you first generate a polynomial of degree 1620, then apply 1620 times iterative partial derivatives to find all coefficients of power terms that are perfect squares, it could be computationally heavy.

Final computations would be done mod 10^9 since we only need the last 9 digits. This should greatly speed up the calculation. This works because of the distributive laws which are preserved under modulo arithmetic.
Though this solution needs some knowledge of generating function and could involve quite an amount of coding depending on the language used, it is one of the problems that can be methodologically thought out once the use of polynomial to generate the sequence is realized.

More Answers: