Clearly for $x$ coprime to $1000$, the only solutions are the trivial $x = 0$ and $x = 1$, so all nontrivial solutions for $x$ must be divisible by either $2$ or $5$, or both. This leaves us with about $600$ possible candidates. Now it should be clear that if a solution works for $1000$, it will also work for $100$, and for $10$. Let's try and find the solutions for $10$:
$$x^2 \equiv x \pmod{10} \tag{2}$$
We have $5$ potential solutions here, namely $2, 4, 5, 6, 8$ (modulo $10$). Let's check each of them manually, some calculations show that only $5$ and $6$ work. So the solutions must have either $5$ or $6$ as a last digit. Let's now extend our search to $100$, armed with this information (which reduces the set of possible solutions considerably):
$$x^2 \equiv x \pmod{100} \tag{3}$$
We know from the previous step that the potential solutions are $5, 6, 15, 16, \cdots$. A short exhaustive search tells us that only $25$ and $76$ work. Finally, we have narrowed down the search enough to look for the real solutions:
$$x^2 \equiv x \pmod{1000} \tag{4}$$
The only possible solutions are $25, 76, 125, 176, \cdots$. Again, trying those 20 potential solutions shows that only $376$ and $625$ work. Therefore, the solutions are $x \in \{ 0, \, 1, \, 376, \, 625 \}$. QED.
This might look like a non-answer given the amount of trial and error involved, but consider that we've only had to try $45$ potential candidates instead of the original $600$, and the cost of this algorithm is asymptotically logarithmic, for instance, finding the solutions with $m = 100000$ would only take $85$ trials, quite efficient.
Here's an Python 3 script implementing the algorithm outlined above, using recursion, for powers of $10$:
Code:
# Work count
work = 0
def Verify(x, m):
global work
work += 1
return pow(x, 2, 10**m) == x
def Solve(m):
# Solutions
found = []
if m == 1:
for x in [2, 4, 5, 6, 8]: # Precomputed
if Verify(x, m):
found.append(x)
return found
possible = Solve(m - 1)
for p in possible:
for x in range(p, 10**m, 10**(m - 1)):
if Verify(x, m):
found.append(x)
return found
def Enumerate(m):
global work
work = 0
found = Solve(m)
# Add trivial solutions
found.append(0)
found.append(1)
found.sort()
print("Solutions are {" + ', '.join(str(x) for x in found) + "}.")
print("Checked " + str(work) + " candidates.")
Which gives us (spaces added for readability):
Code:
>>> Enumerate(1) # For 10
Solutions are {0, 1, 5, 6}.
Checked 5 candidates.
>>> Enumerate(2) # For 100
Solutions are {0, 1, 25, 76}.
Checked 25 candidates.
>>> Enumerate(3) # For 1000
Solutions are {0, 1, 376, 625}.
Checked 45 candidates.
>>> Enumerate(5) # For 100000
Solutions are {0, 1, 9376, 90625}.
Checked 85 candidates.
>>> Enumerate(8) # For 100000000
Solutions are {0, 1, 12890625, 87109376}.
Checked 145 candidates.
Interestingly, there are only two nontrivial solutions for any power of $10$, and these two solutions add up to $-1$ modulo said power, suggesting that there might exist an analytic solution to this problem, though the prime factorisations of the solutions display no obvious pattern apart from a large number of powers of $2$ and $5$.
So if a better way to do this exists, I'd be interested to see it. It is worth noting that the condition $x^2 \equiv x \pmod{m}$ implies $x^n \equiv x \pmod{m}$ for all $n > 0$, but I could not find how to use this property, sadly.