Problem 269: Polynomials with at Least One Integer Root

Project Euler Problem 269 Solution

EulerSolve provides an optimized solution for Project Euler Problem 269, Polynomials with at Least One Integer Root, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write a positive integer $n$ in decimal as $d_{L-1}d_{L-2}\dots d_1d_0$, with $d_{L-1}\neq0$. The problem associates to $n$ the polynomial $$P_n(x)=d_{L-1}x^{L-1}+d_{L-2}x^{L-2}+\cdots+d_1x+d_0.$$ We must count how many integers $1\le n\le10^{16}$ have the property that $P_n(x)$ has at least one integer root. Mathematical Approach 1. Candidate Integer Roots Are Extremely Limited Because every coefficient is a decimal digit, all coefficients are nonnegative and the leading coefficient is positive. Therefore $P_n(r)>0$ for every positive integer $r$. So positive roots are impossible. By the rational-root theorem, any integer root must divide the constant term $d_0$. Since $d_0\in\{0,1,\dots,9\}$, the only possibilities are: 1. $r=0$, which can happen only when $d_0=0$, 2. $r=-t$, where $t$ is a positive divisor of $d_0$, so $1\le t\le 9$. This is the first major reduction: the code never searches an unbounded root range. 2. The Root Condition for $r=-t$ Fix a candidate $t\ge1$. The condition $P_n(-t)=0$ is $$\sum_{i=0}^{L-1} d_i(-t)^i=0.$$ Group digits in pairs from the least significant side: $$e_j=d_{2j},\qquad o_j=d_{2j+1},$$ where a missing top odd digit is treated as $0$. Then $$P_n(-t)=\sum_{j\ge0}(e_j-to_j)t^{2j}.$$ The whole task is now to make those pair contributions cancel exactly. 3....

Detailed mathematical approach

Problem Summary

Write a positive integer $n$ with $L$ decimal digits as $d_{L-1}\dots d_1d_0$, so that $d_0$ is its last digit and $d_{L-1}\ne0$. Its digits become the coefficients of the polynomial

$$P_n(x)=d_{L-1}x^{L-1}+\dots+d_1x+d_0;$$

for example $P_{5703}(x)=5x^3+7x^2+3$, and in general $P_n(0)$ is the last digit of $n$, $P_n(1)$ is its digit sum and $P_n(10)=n$. Let $Z(k)$ be the number of positive integers $n\le k$ for which $P_n$ has at least one integer root, that is, an integer $x$ with $P_n(x)=0$. The statement gives $Z(100\,000)=14696$, and we must find $Z(10^{16})$. ^[1]

The whole solution in one sentence: an integer root can only be $0$, which is a root exactly when $n$ ends in $0$, or $-t$ for a divisor $t$ of the last digit; each such $-t$ can be tested on the digits from right to left with a small carry, as in written addition, so grouping the digit choices by their carries counts the numbers with given roots, and inclusion–exclusion counts every number with at least one root exactly once.

Mathematical Approach

We first look at small numbers, then show that only a few integers can be roots, then test a single candidate with carries, and finally count all numbers of a given length at once.

1. Small cases

A one-digit number gives the constant polynomial $P_n(x)=n\ne0$, which has no root. A two-digit number $n=10a+b$ gives $P_n(x)=ax+b$, whose only root $-b/a$ is an integer exactly when $a$ divides $b$. This holds for the $9$ numbers $10,20,\dots,90$ and for the $23$ numbers $11$ to $19$, $22$, $24$, $26$, $28$, $33$, $36$, $39$, $44$, $48$, $55$, $66$, $77$, $88$ and $99$. Together with $100$, whose polynomial $x^2$ has the root $0$, this gives $Z(100)=9+23+1=33$. Among the three-digit numbers, $P_{121}(x)=(x+1)^2$ has the root $-1$ twice, and $P_{132}(x)=x^2+3x+2$${}=(x+1)(x+2)$ has the two roots $-1$ and $-2$, so adding up the numbers with each root would count $132$ twice.

A computer search that tests every integer $x$ with $|x|\le10$ in every polynomial gives $Z(10^k)$ for $k=1,\dots,8$: $1$, $33$, $172$, $1754$, $14696$, $152960$, $1393262$ and $14435778$, in agreement with the value of the statement for $k=5$. Larger $|x|$ need not be tested: for $|x|\ge10$ the leading term outweighs all the others, since $9\,(|x|^{L-2}+\dots+1)\lt|x|^{L-1}$.

The examples share a pattern. Two integers $u$ and $v$ are congruent modulo $m$, written $u\equiv v\pmod m$, when $u-v$ is a multiple of $m$; congruences survive sums and products, and therefore the evaluation of a polynomial with integer coefficients. ^[3] Since $10\equiv-t\pmod{10+t}$,

$$n=P_n(10)\equiv P_n(-t)\pmod{10+t},$$

so every $n$ with the root $-t$ is a multiple of $10+t$: $132=11\cdot12$ has the roots $-1$ and $-2$, and $19$, with $P_{19}(x)=x+9$, has the root $-9$. For $t=1$, $P_n(-1)$ is the alternating sum of the digits that the divisibility rule for $11$ uses, ^[4] but a root needs this sum to be exactly $0$, not merely a multiple of $11$: the smallest multiple of $11$ without the root $-1$ is $209$, with $P_{209}(-1)=11$.

2. Where an integer root can lie

Claim. The number $0$ is a root of $P_n$ exactly when $d_0=0$, and every other integer root is $-t$ for a positive divisor $t$ of $d_0$.

Proof. The coefficients are digits, so none of them is negative and the leading one is at least $1$; hence $P_n(x)\ge x^{L-1}\gt0$ for $x\gt0$, and there is no positive root. Next, $P_n(0)=d_0$. Finally, if $P_n(-t)=0$ with $t\ge1$, then

$$d_0=-\sum_{i=1}^{L-1}d_i(-t)^i=t\sum_{i=1}^{L-1}(-1)^{i+1}d_i\,t^{i-1},$$

so $t$ divides $d_0$. This is the integer case of the rational root theorem, by which a rational root $p/q$ in lowest terms has $p$ dividing the constant term and $q$ dividing the leading coefficient. ^[2] $\square$

A number that ends in $0$ is therefore counted because of the root $0$, whatever its other roots are. For a last digit $d_0\ne0$ the candidates are $-t$ with $t$ in the set $D(d_0)$ of the divisors of $d_0$: $D(1)=\{1\}$, $D(2)=\{1,2\}$, $D(3)=\{1,3\}$, $D(4)=\{1,2,4\}$, $D(5)=\{1,5\}$, $D(6)=\{1,2,3,6\}$, $D(7)=\{1,7\}$, $D(8)=\{1,2,4,8\}$ and $D(9)=\{1,3,9\}$. For instance, $5703$ ends in $3$, and since $P_{5703}(-1)=5$ and $P_{5703}(-3)=-69$, it has no integer root.

3. Numbers that end in 0

For $L\ge2$ there are $9\cdot10^{L-2}$ numbers with $L$ digits whose last digit is $0$: $9$ choices for the leading digit and $10$ for each of the $L-2$ digits in between. Together with $10^k$ itself, whose polynomial $x^k$ has the root $0$, the numbers up to $10^k$ that end in $0$ contribute

$$1+\sum_{L=2}^{k}9\cdot10^{L-2}=10^{k-1}$$

to $Z(10^k)$, one for each multiple of $10$. From now on the last digit is $d_0\ne0$.

4. Testing one candidate with carries

Fix $t\ge1$. Because $(-t)^{2j}=t^{2j}$ and $(-t)^{2j+1}=-t\cdot t^{2j}$, it pays to read the digits in pairs from the right: let $e_j=d_{2j}$ and $o_j=d_{2j+1}$ for $j=0,\dots,s-1$, where $s=\lceil L/2\rceil$ and $o_{s-1}=0$ when $L$ is odd. Then

$$P_n(-t)=\sum_{j=0}^{s-1}a_j\,t^{2j},\qquad a_j=e_j-t\,o_j.$$

For $t\ge2$ this looks like a number written in base $t^2$, except that its "digits" $a_j$ lie between $-9t$ and $9$. Whether such a sum is $0$ can be checked from the right, as in written addition, where each column passes a carry to the next one. Start with $c_0=0$ and set

$$c_{j+1}=\frac{t\,o_j+c_j-e_j}{t^2}=\frac{c_j-a_j}{t^2},$$

a step that is allowed only when the division is exact.

Claim. $P_n(-t)=0$ exactly when all $s$ divisions are exact and the last carry $c_s$ is $0$.

Proof. If the first $r$ divisions are exact, then $a_j\,t^{2j}=c_j\,t^{2j}-c_{j+1}\,t^{2j+2}$ for $j\lt r$, and in the sum of these terms all carries except the first and the last cancel:

$$\sum_{j=0}^{r-1}a_j\,t^{2j}=c_0-c_r\,t^{2r}=-c_r\,t^{2r}.\qquad(\ast)$$

For $r=s$ this says $P_n(-t)=-c_s\,t^{2s}$, so exact divisions and $c_s=0$ give a root. Conversely, let $P_n(-t)=0$, and suppose that the first $r\lt s$ divisions are exact, which is certainly true for $r=0$. Adding $a_r\,t^{2r}$ to both sides of $(\ast)$ and using that all $s$ terms together give $0$, we get

$$c_r-a_r=\sum_{j=r+1}^{s-1}a_j\,t^{2(j-r)},$$

which is divisible by $t^2$. So the next division is exact too, all $s$ divisions are exact, and $c_s=-P_n(-t)/t^{2s}=0$. $\square$

Formula $(\ast)$ also shows that the carries stay small, for any digits, as long as the divisions are exact. Since $-9t\le a_j\le9$ and $\sum_{j\lt r}t^{2j}=(t^{2r}-1)/(t^2-1)$, it gives for $t\ge2$

$$-\frac{9}{t^2-1}\lt c_r\lt\frac{9t}{t^2-1}.$$

So the carry of $t=2$ is one of $-2,\dots,5$, that of $t=3$ one of $-1,\dots,3$, that of $t=4$ one of $0,1,2$, and for $5\le t\le9$ it is $0$ or $1$. For $t=1$ the carry $c_r=\sum_{j\lt r}(o_j-e_j)$ is the running alternating digit sum with the opposite sign, so $|c_r|\le9r$.

5. Worked example

For $n=132$ the pairs are $(e_0,o_0)=(2,3)$ and $(e_1,o_1)=(1,0)$. For $t=1$ the carries are $c_1=(3+0-2)/1=1$ and $c_2=(0+1-1)/1=0$; for $t=2$ they are $c_1=(6+0-2)/4=1$ and $c_2=(0+1-1)/4=0$. Both last carries are $0$, so $-1$ and $-2$ are roots, as the factorisation $(x+1)(x+2)$ shows. For $n=122$, with the pairs $(2,2)$ and $(1,0)$, the candidate $t=1$ gives $c_1=0$ and $c_2=-1\ne0$, and indeed $P_{122}(-1)=1$; for $t=2$ already the first division, $(4+0-2)/4$, is not exact, and indeed $P_{122}(-2)=2$.

6. A table of carries

Fix a length $L$, a last digit $d_0\ne0$ and a nonempty set $S\subseteq D(d_0)$, and let $N(L,d_0,S)$ be the number of $L$-digit numbers ending in $d_0$ that have every $-t$ with $t\in S$ as a root. By Section 4 all these roots can be tested at once by running the carries of all $t\in S$ side by side: the carry vector $(c^{(t)})_{t\in S}$ starts at zero, every digit pair must make all divisions exact, and the last vector must be zero again. The next vector depends only on the current vector and the current pair of digits, not on the digits before. So all choices of the first pairs that lead to the same vector can be merged: we keep a table that maps every vector to the number of ways to reach it and update the table pair by pair. This is dynamic programming, which counts solutions without visiting them one by one. ^[5] The digits are restricted in three places: $e_0=d_0$ is fixed, the leading digit $d_{L-1}$ runs from $1$ to $9$, and $o_{s-1}=0$ when $L$ is odd; every other digit runs from $0$ to $9$. After the last pair, the entry of the zero vector is $N(L,d_0,S)$.

For three-digit numbers ending in $2$, where $D(2)=\{1,2\}$, the tables can be followed by hand. With $S=\{1\}$, the first pair $(2,o_0)$ leads to the carry $o_0-2$, so the table holds the carries $-2,\dots,7$, one way each. The last pair is $(e_1,0)$ with a leading digit $1\le e_1\le9$, and the last carry $c_1-e_1$ is $0$ exactly when $e_1=c_1$, which is possible for $c_1=1,\dots,7$. So $N(3,2,\{1\})=7$: the numbers $132,242,\dots,792$. With $S=\{2\}$, the division $(2o_0-2)/4$ is exact only for odd $o_0$, with the carries $0,\dots,4$, and the last carry $(c_1-e_1)/4$ is $0$ exactly when $e_1=c_1$, so $N(3,2,\{2\})=4$: the numbers $132$, $252$, $372$ and $492$. With $S=\{1,2\}$, the odd digits $o_0=1,3,5,7,9$ lead to the vectors $(-1,0)$, $(1,1)$, $(3,2)$, $(5,3)$ and $(7,4)$, and only $(1,1)$ with $e_1=1$ ends at zero, so $N(3,2,\{1,2\})=1$: the number $132$.

7. Inclusion–exclusion over the candidates

For a length $L$ and a last digit $d_0\ne0$, let $A_t$ be the set of $L$-digit numbers ending in $d_0$ with the root $-t$. By Section 2 we need the size $U(L,d_0)$ of the union of the sets $A_t$ with $t\in D(d_0)$, and $N(L,d_0,S)$ is the size of the intersection of the sets $A_t$ with $t\in S$. The inclusion–exclusion principle ^[6] gives

$$U(L,d_0)=\sum_{\emptyset\ne S\subseteq D(d_0)}(-1)^{|S|+1}\,N(L,d_0,S).$$

For three-digit numbers ending in $2$ this is $U(3,2)=7+4-1=10$: the number $132$ lies in both sets and is subtracted once. A repeated root causes no trouble; $121$ is simply one element of $A_1$. The sets $D(1),\dots,D(9)$ have $1,2,2,3,2,4,2,4,3$ elements, so a length needs $2^{|D(d_0)|}-1$ intersection counts for each last digit, $57$ in all.

8. Which intersections can be nonempty

Many of these counts are $0$ for a simple reason. By the factor theorem, a root $-t_1$ splits off the factor $x+t_1$, and the long division by $x+t_1$ needs no fractions, because its leading coefficient is $1$. ^[7] A second root $-t_2\ne-t_1$ is then a root of the quotient, since $x+t_1$ does not vanish at $-t_2$, and so on. Distinct roots $-t_1,\dots,-t_m$ therefore give

$$P_n(x)=(x+t_1)\cdots(x+t_m)\,Q(x)$$

with a polynomial $Q$ with integer coefficients, and $x=0$ gives $d_0=t_1\cdots t_m\,Q(0)$: the product $t_1\cdots t_m$ divides $d_0$. For $d_0=6$, only $9$ of the $15$ sets $S$ can have a nonzero count, namely $\{1\}$, $\{2\}$, $\{3\}$, $\{6\}$, $\{1,2\}$, $\{1,3\}$, $\{1,6\}$, $\{2,3\}$ and $\{1,2,3\}$, and the same holds for $d_0=8$ with $3$ and $6$ replaced by $4$ and $8$. The program does not use this shortcut; for the other sets its table simply ends without the zero vector. With $16$ digits, every set whose product divides $d_0$ does have a nonzero count.

9. Assembling the count

One-digit numbers have no root, the lengths $L=1,\dots,k$ cover all numbers below $10^k$, and $10^k$ itself was counted in Section 3. Hence

$$Z(10^k)=10^{k-1}+\sum_{L=2}^{k}\sum_{d_0=1}^{9}U(L,d_0).$$

For $k=3$, the values $U(2,1),\dots,U(2,9)$ are $1,2,2,3,2,4,2,4,3$, together the $23$ numbers of Section 1, and $U(3,1),\dots,U(3,9)$ are $8,10,7,8,4,5,2,3,2$, together $49$. So $Z(1000)=100+23+49=172$, the value that the C++ program checks against a direct search.

How the Code Works

The C++ program ^[9] starts in main, where parse_arguments reads the options --power=$k$, with the default $k=16$, and --skip-checkpoints; any other argument, and $k=0$, ends the program with exit code $1$. The parser reads $k$ into an int without an overflow check, and the 64-bit counts give correct results only up to $k=20$. Unless the checkpoints are skipped, run_checkpoints requires that solve(3) equals brute_count(1000), that solve(4) equals brute_count(10000) and that solve(5) equals $14696$, and ends the program with exit code $2$ otherwise. Then main prints solve of the chosen power, which adds count_for_length for the lengths $1,\dots,k$ and finally $1$ for $10^k$ itself (Section 9).

count_for_length starts with $9\cdot10^{L-2}$ from pow10i when $L\ge2$ (Section 3). For every last digit $d_0=1,\dots,9$ it collects the divisors $t\le9$ of $d_0$ (Section 2), turns every bit mask from $1$ to $2^m-1$, where $m$ is the number of these divisors, into the list roots of the chosen divisors, and adds or subtracts count_subset_for_length of that list according to the parity of __builtin_popcount(mask), in the signed 64-bit union_count (Section 7). count_subset_for_length is the table of Section 6: a std::map from carry vectors std::vector<int> to 64-bit counts, sorted by the vectors, so that a lookup costs a logarithmic number of comparisons. ^[10] It starts with the zero vector and the count $1$. For each of the even_slots $=\lceil L/2\rceil$ pairs it lists the allowed even_digits and odd_digits as in Section 6, without $0$ at the leading position msd_pos. For every vector and digit pair it computes num $=t\,o_j+c_j-e_j$ in long long for each root, drops the pair as soon as num % denom with denom $=t^2$ is not $0$, and otherwise adds the count to the new vector in next_states. If no vector survives a pair, the function returns $0$ at once; at the end it returns the entry of the zero vector, or $0$.

has_integer_root, the direct search of the checkpoints, splits $n$ into digits, tries $0$ when the last digit is $0$ and $-t$ for every divisor $t\le9$ of the last digit, and evaluates the polynomial with Horner's method, value = value * r + d, in long long; ^[8] brute_count counts the numbers that pass. This search relies on the claim of Section 2, so it checks the carries, the tables and the inclusion–exclusion, but not the claim itself.

The Python and Java programs have the same functions (count_subset_for_length, count_for_length, solve_power in Python; countSubsetForLength, countForLength, solve in Java) and compute $Z(10^{16})$ directly, without options or checkpoints. Python keys a dictionary by tuples of carries and counts mask bits with int.bit_count, which needs Python 3.10 or later. Java wraps the carry array in the class StateKey, whose equals and hashCode compare contents via Arrays.equals and Arrays.hashCode, so that it can be a HashMap key, and counts bits with Integer.bitCount. Carries can be negative, and then C++ and Java round an integer quotient toward $0$ while Python rounds it toward $-\infty$; ^[10] ^[11] but every program divides only after checking that the remainder is $0$, so all three get the same carries. For $k=16$ no count exceeds $10^{16}\lt2^{63}$, so the 64-bit integers of C++ and Java suffice; Python's integers are unbounded anyway.

Complexity and Verification

For the bound $10^k$ the program builds one table for each of the $k$ lengths and $57$ sets, with $\lceil L/2\rceil$ pair steps; in each step every vector tries at most $100$ digit pairs, and every pair costs at most $|S|\le4$ divisibility tests. By Section 4 the carry of $t=1$ takes at most $18r+1$ values after $r$ pairs and every other carry at most $8$, so a table never holds more than a constant times $k$ vectors. The work therefore grows at most like $k^3$ table operations, each of which costs a logarithmic number of vector comparisons in the C++ std::map, and since only the current and the next table are kept, the memory grows like $k$. For $k=16$ there are $912$ such runs, which examine about $7.0$ million combinations of a vector with a digit pair and make about $11.8$ million divisibility tests; no table holds more than $716$ vectors, a maximum reached for $L=16$, $d_0=2$ and $S=\{1,2\}$. In our runs the C++ program, checkpoints included, and the Java program finished in about $0.2$ seconds and the Python program in about $1.3$ seconds. The direct search of Section 1 needed about $5$ seconds for $10^8$ numbers; at the same cost per number, $10^{16}$ numbers would take $10^8$ times as long, more than a decade.

All three programs print the same number, and the C++ program passes its checkpoints. The following was verified for this article:

The search of Section 1 over all $n\le10^8$, which tests every integer $|x|\le10$ and so does not rely on Section 2, reproduces solve for every power $k\le8$; it found no positive root and no root with $|x|=10$.
A search over all numbers with at most $6$ digits reproduces all $342$ intersection counts $N(L,d_0,S)$ with $L\le6$.
For every length up to $16$, the sets $S$ whose product does not divide $d_0$ have the count $0$, as Section 8 predicts. The carries stay within the bounds of Section 4: for $t\ge2$ they reach both ends of each range, and for $t=1$ they run from $-72$ to $71$.
A separate program, written for this article, accepts every number ending in $0$ because of the root $0$. For the other numbers it uses single-digit carries, $c_0=0$ and $c_{i+1}=-(d_i+c_i)/t$, for which $P_n(-t)=c_L\,(-t)^L$, and follows all nine candidates $-1,\dots,-9$ in one table, marking a candidate as failed when a division is not exact and accepting a number when a surviving candidate ends with the carry $0$. Without digit pairs, inclusion–exclusion or the divisor sets of Section 2, it gives the same $Z(10^k)$ as the C++ program for every $k\le16$, and the same value as all three programs for $k=16$.

Footnotes and References

The references below supply the background facts used above. The carry test with its proof and bounds, the counting tables, the restriction of the intersections and all counts are derived in this article.

Project Euler 269 — Polynomials with at Least One Integer Root. The official statement defines Pₙ as the polynomial whose coefficients are the digits of n, for example P₅₇₀₃(x) = 5x³ + 7x² + 3, defines Z(k) as the number of positive integers n not exceeding k for which Pₙ has at least one integer root, gives Z(100 000) = 14696 and asks for Z(10¹⁶).
Wikipedia — Rational root theorem. A rational root p/q in lowest terms of a polynomial with integer coefficients has p dividing the constant term and q dividing the leading coefficient. For an integer root q = 1, so the root divides the constant term, as used in Section 2.
Wikipedia — Modular arithmetic. Congruence is compatible with addition and multiplication, and a ≡ b (mod m) implies p(a) ≡ p(b) (mod m) for every polynomial p with integer coefficients; Section 1 applies this with a = 10, b = −t and m = 10 + t.
Wikipedia — Divisibility rule. A number is divisible by 11 exactly when the alternating sum of its digits is divisible by 11; the article explains this rule with 10 ≡ −1 (mod 11).
Wikipedia — Dynamic programming. The method breaks a problem into simpler sub-problems whose solutions are stored and reused; the article notes that it can count solutions without visiting them all, as the tables of Section 6 do.
Wikipedia — Inclusion–exclusion principle. The size of a union of finitely many finite sets is the sum, over all nonempty families of these sets, of the sizes of their intersections, with the sign + for an odd and − for an even number of sets.
Wikipedia — Factor theorem. A polynomial f(x) has the factor x − a exactly when f(a) = 0. This holds over any commutative ring, in particular over the integers, because x − a is monic, so the long division needs no division of coefficients.
Wikipedia — Horner's method. A polynomial of degree n is evaluated in nested form, starting with the leading coefficient and repeatedly multiplying by x and adding the next coefficient, with n multiplications and n additions.
C++ — Euler269.cpp. The linked immutable C++ revision contains count_subset_for_length, count_for_length, solve, has_integer_root, brute_count and run_checkpoints with its three checkpoints. It is the implementation record for this article.
cppreference — std::map: a sorted associative container whose search, removal and insertion operations have logarithmic complexity. cppreference — Arithmetic operators: built-in integer division truncates the quotient towards zero.
Java Language Specification, §15.17.2: integer division rounds toward 0. Python — Numeric types: x // y is always rounded towards minus infinity, integers have unlimited precision, and int.bit_count was added in version 3.10.

Problem 269 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <map>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>

namespace {

using u64 = std::uint64_t;

struct Options {
    int max_power = 16;
    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, "--power=", options.max_power)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_power >= 1;
}

u64 pow10i(const int exp) {
    u64 value = 1;
    for (int i = 0; i < exp; ++i) {
        value *= 10;
    }
    return value;
}

u64 count_subset_for_length(const int length, const int last_digit, const std::vector<int>& roots) {
    if (roots.empty()) {
        return 0;
    }

    const int even_slots = (length + 1) / 2;
    const int msd_pos = length - 1;

    const std::vector<int> zero_state(static_cast<std::size_t>(roots.size()), 0);
    std::map<std::vector<int>, u64> states;
    states[zero_state] = 1;

    for (int j = 0; j < even_slots; ++j) {
        const int even_pos = 2 * j;
        const int odd_pos = 2 * j + 1;

        std::vector<int> even_digits;
        if (j == 0) {
            even_digits.push_back(last_digit);
        } else {
            even_digits.reserve(10);
            for (int d = 0; d <= 9; ++d) {
                even_digits.push_back(d);
            }
        }
        if (even_pos == msd_pos) {
            std::vector<int> filtered;
            filtered.reserve(even_digits.size());
            for (int d : even_digits) {
                if (d != 0) {
                    filtered.push_back(d);
                }
            }
            even_digits.swap(filtered);
        }

        std::vector<int> odd_digits;
        if (odd_pos >= length) {
            odd_digits.push_back(0);
        } else {
            odd_digits.reserve(10);
            for (int d = 0; d <= 9; ++d) {
                odd_digits.push_back(d);
            }
            if (odd_pos == msd_pos) {
                std::vector<int> filtered;
                filtered.reserve(odd_digits.size());
                for (int d : odd_digits) {
                    if (d != 0) {
                        filtered.push_back(d);
                    }
                }
                odd_digits.swap(filtered);
            }
        }

        std::map<std::vector<int>, u64> next_states;
        for (const auto& kv : states) {
            const std::vector<int>& carry = kv.first;
            const u64 ways = kv.second;

            for (int ed : even_digits) {
                for (int od : odd_digits) {
                    std::vector<int> next(carry.size(), 0);
                    bool ok = true;
                    for (std::size_t idx = 0; idx < roots.size(); ++idx) {
                        const int t = roots[idx];
                        const int denom = t * t;
                        const long long num = static_cast<long long>(t) * od + carry[idx] - ed;
                        if (num % denom != 0) {
                            ok = false;
                            break;
                        }
                        next[idx] = static_cast<int>(num / denom);
                    }
                    if (!ok) {
                        continue;
                    }
                    next_states[next] += ways;
                }
            }
        }

        states.swap(next_states);
        if (states.empty()) {
            return 0;
        }
    }

    const auto it = states.find(zero_state);
    return (it == states.end()) ? 0 : it->second;
}

u64 count_for_length(const int length) {
    u64 total = 0;

    // Root x=0: exactly numbers ending in zero.
    if (length >= 2) {
        total += 9ULL * pow10i(length - 2);
    }

    for (int d0 = 1; d0 <= 9; ++d0) {
        std::vector<int> divisors;
        for (int t = 1; t <= 9; ++t) {
            if (d0 % t == 0) {
                divisors.push_back(t);
            }
        }

        const int m = static_cast<int>(divisors.size());
        long long union_count = 0;
        for (int mask = 1; mask < (1 << m); ++mask) {
            std::vector<int> roots;
            roots.reserve(static_cast<std::size_t>(m));
            for (int i = 0; i < m; ++i) {
                if (((mask >> i) & 1) != 0) {
                    roots.push_back(divisors[static_cast<std::size_t>(i)]);
                }
            }

            const u64 cnt = count_subset_for_length(length, d0, roots);
            if ((__builtin_popcount(static_cast<unsigned int>(mask)) & 1) != 0) {
                union_count += static_cast<long long>(cnt);
            } else {
                union_count -= static_cast<long long>(cnt);
            }
        }

        total += static_cast<u64>(union_count);
    }

    return total;
}

u64 solve(const int max_power) {
    u64 total = 0;
    for (int len = 1; len <= max_power; ++len) {
        total += count_for_length(len);
    }

    // Add 10^max_power itself, which always has root x=0.
    total += 1;
    return total;
}

bool has_integer_root(u64 n) {
    std::vector<int> digits;
    while (n > 0) {
        digits.push_back(static_cast<int>(n % 10ULL));
        n /= 10ULL;
    }
    std::reverse(digits.begin(), digits.end());

    const int d0 = digits.back();
    std::vector<int> roots;
    if (d0 == 0) {
        roots.push_back(0);
    }
    for (int t = 1; t <= 9; ++t) {
        if (d0 % t == 0) {
            roots.push_back(-t);
        }
    }

    for (int r : roots) {
        long long value = 0;
        for (int d : digits) {
            value = value * static_cast<long long>(r) + d;
        }
        if (value == 0) {
            return true;
        }
    }
    return false;
}

u64 brute_count(const u64 limit) {
    u64 count = 0;
    for (u64 n = 1; n <= limit; ++n) {
        if (has_integer_root(n)) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(3) != brute_count(1000)) {
        std::cerr << "Checkpoint failed for limit 1000" << '\n';
        return false;
    }
    if (solve(4) != brute_count(10000)) {
        std::cerr << "Checkpoint failed for limit 10000" << '\n';
        return false;
    }
    if (solve(5) != 14696ULL) {
        std::cerr << "Checkpoint failed for Z(100000)" << '\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.max_power) << '\n';
    return 0;
}

Python

def pow10i(exp):
    return 10**exp

def count_subset_for_length(length, last_digit, roots):
    if not roots:
        return 0

    even_slots = (length + 1) // 2
    msd_pos = length - 1

    zero_state = tuple(0 for _ in roots)
    states = {zero_state: 1}

    for j in range(even_slots):
        even_pos = 2 * j
        odd_pos = 2 * j + 1

        if j == 0:
            even_digits = [last_digit]
        else:
            even_digits = list(range(10))

        if even_pos == msd_pos:
            even_digits = [d for d in even_digits if d != 0]

        if odd_pos >= length:
            odd_digits = [0]
        else:
            odd_digits = list(range(10))
            if odd_pos == msd_pos:
                odd_digits = [d for d in odd_digits if d != 0]

        next_states = {}
        for carry, ways in states.items():
            for ed in even_digits:
                for od in odd_digits:
                    next_carry = []
                    ok = True
                    for idx, t in enumerate(roots):
                        denom = t * t
                        num = t * od + carry[idx] - ed
                        if num % denom != 0:
                            ok = False
                            break
                        next_carry.append(num // denom)
                    
                    if not ok:
                        continue
                        
                    next_carry_tuple = tuple(next_carry)
                    next_states[next_carry_tuple] = next_states.get(next_carry_tuple, 0) + ways

        states = next_states
        if not states:
            return 0

    return states.get(zero_state, 0)

def count_for_length(length):
    total = 0
    if length >= 2:
        total += 9 * pow10i(length - 2)

    for d0 in range(1, 10):
        divisors = [t for t in range(1, 10) if d0 % t == 0]
        m = len(divisors)
        union_count = 0

        for mask in range(1, 1 << m):
            roots = []
            for i in range(m):
                if (mask >> i) & 1:
                    roots.append(divisors[i])

            cnt = count_subset_for_length(length, d0, roots)
            if mask.bit_count() & 1:
                union_count += cnt
            else:
                union_count -= cnt

        total += union_count

    return total

def solve_power(max_power):
    total = 0
    for length in range(1, max_power + 1):
        total += count_for_length(length)

    total += 1
    return total

def solve():
    ans = solve_power(16)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler269 {

    static long pow10i(int exp) {
        long value = 1;
        for (int i = 0; i < exp; i++)
            value *= 10;
        return value;
    }

    static class StateKey {
        int[] carry;

        StateKey(int[] c) {
            this.carry = c.clone();
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            StateKey stateKey = (StateKey) o;
            return Arrays.equals(carry, stateKey.carry);
        }

        @Override
        public int hashCode() {
            return Arrays.hashCode(carry);
        }
    }

    static long countSubsetForLength(int length, int lastDigit, List<Integer> roots) {
        if (roots.isEmpty())
            return 0;

        int evenSlots = (length + 1) / 2;
        int msdPos = length - 1;

        int[] zeroCarry = new int[roots.size()];
        StateKey zeroState = new StateKey(zeroCarry);
        Map<StateKey, Long> states = new HashMap<>();
        states.put(zeroState, 1L);

        for (int j = 0; j < evenSlots; j++) {
            int evenPos = 2 * j;
            int oddPos = 2 * j + 1;

            List<Integer> evenDigits = new ArrayList<>();
            if (j == 0) {
                evenDigits.add(lastDigit);
            } else {
                for (int d = 0; d <= 9; d++)
                    evenDigits.add(d);
            }
            if (evenPos == msdPos) {
                evenDigits.removeIf(d -> d == 0);
            }

            List<Integer> oddDigits = new ArrayList<>();
            if (oddPos >= length) {
                oddDigits.add(0);
            } else {
                for (int d = 0; d <= 9; d++)
                    oddDigits.add(d);
                if (oddPos == msdPos) {
                    oddDigits.removeIf(d -> d == 0);
                }
            }

            Map<StateKey, Long> nextStates = new HashMap<>();
            for (Map.Entry<StateKey, Long> entry : states.entrySet()) {
                int[] carry = entry.getKey().carry;
                long ways = entry.getValue();

                for (int ed : evenDigits) {
                    for (int od : oddDigits) {
                        int[] nextCarry = new int[roots.size()];
                        boolean ok = true;
                        for (int idx = 0; idx < roots.size(); idx++) {
                            int t = roots.get(idx);
                            int denom = t * t;
                            long num = (long) t * od + carry[idx] - ed;
                            if (num % denom != 0) {
                                ok = false;
                                break;
                            }
                            nextCarry[idx] = (int) (num / denom);
                        }
                        if (!ok)
                            continue;

                        StateKey nextKey = new StateKey(nextCarry);
                        nextStates.put(nextKey, nextStates.getOrDefault(nextKey, 0L) + ways);
                    }
                }
            }

            states = nextStates;
            if (states.isEmpty())
                return 0;
        }

        return states.getOrDefault(zeroState, 0L);
    }

    static long countForLength(int length) {
        long total = 0;
        if (length >= 2) {
            total += 9L * pow10i(length - 2);
        }

        for (int d0 = 1; d0 <= 9; d0++) {
            List<Integer> divisors = new ArrayList<>();
            for (int t = 1; t <= 9; t++) {
                if (d0 % t == 0)
                    divisors.add(t);
            }

            int m = divisors.size();
            long unionCount = 0;
            for (int mask = 1; mask < (1 << m); mask++) {
                List<Integer> roots = new ArrayList<>();
                for (int i = 0; i < m; i++) {
                    if (((mask >> i) & 1) != 0)
                        roots.add(divisors.get(i));
                }

                long cnt = countSubsetForLength(length, d0, roots);
                if ((Integer.bitCount(mask) & 1) != 0) {
                    unionCount += cnt;
                } else {
                    unionCount -= cnt;
                }
            }
            total += unionCount;
        }

        return total;
    }

    static long solve(int maxPower) {
        long total = 0;
        for (int len = 1; len <= maxPower; len++) {
            total += countForLength(len);
        }
        total += 1;
        return total;
    }

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