Enigma 1736, or solving $x^3=x$ mod $10^N$

The problem posed by New Scientist Enigma 1736 was to find the largest 4 digit number with all digits different whose cube has the same last 4 digits as the original number.

Generalising the problem to finding an N digit number whose cube has the same last N digits as the original number, i.e. solutions of $x^3=x$ mod $10^N$, (without the constraint of all digits being different) we can rewrite the equation as

$x^3-x = (x-1)x(x+1)=0$ mod $10^N$

For $N$ $\geq$ $3$ there are always fifteen solutions, one for each of the cells in the table. Five (in grey) have an immediate solution, the other ten can be found by solving three linear diophantine equations (red, green and blue).

		$5^N | (x-1)$	$5^N |x $	$5^N | (x+1)$
$2^N | (x-1) $	$2|(x+1)$	$1$	$(r$ mod $2^N)5^N$ $r5^N+s2^N=1$	$2(r$ mod $2^{N-1})5^N-1$ $r5^N+s2^{N-1}=1$
$2^N| x $		$10^N-(r$ mod $2^N)5^N+1$ $r5^N+s2^N=1$	$0$	$(r$ mod $2^N)5^N-1$ $r5^N+s2^N=1$
$2^N| (x+1)$	$2 | (x-1)$	$10^N-2(r$ mod $2^{N-1})5^N+1$ $r5^N+s2^{N-1}=1$	$10^N-(r$ mod $2^N)5^N$ $r5^N+s2^N=1$	$10^N-1$
$2^{N-1}| (x-1)$	$ 2 | (x+1)$	$5^N 2^{N-1}+1$	$(r$ mod $2^{N-1})5^N$ $r5^N+s2^{N-1}=1$	$2(r$ mod $2^{N-2})5^N-1$ $r5^N+s2^{N-2}=1$
$2^{N-1}| (x+1)$	$2 | (x-1)$	$10^N-2(r$ mod $2^{N-2})5^N+1$ $r5^N+s2^{N-2}=1$	$10^N-(r$ mod $2^{N-1})5^N$ $r5^N+s2^{N-1}=1$	$5^N 2^{N-1}-1$

Complementary pairs of solutions, $x$ and $10^N-x$ occur on opposite sides, as illustrated by the example for N=4

		625 | (x−1)	625 | x	625 | (x+1)
16 | (x−1)		1	625	1249
16 | x		9376	0	624
16 | (x+1)		8751	9375	9999
8 | (x−1)		5001	5625	6249
8 | (x+1)		3751	4375	4999

Here is a Python implementation to find solutions:

def egcd(a,b):
  if b == 0:
    return [1,0,a]
  else:
    x,y,g = egcd(b, a%b)
    return [y, x - (a//b)*y, g] 

for N in range(3,101):
  solutions = {0,1,10**N-1, (10**N)//2-1,(10**N)//2+1}
  
  # Blue diophantine 
  r,s,g = egcd(5**N, 2**N)
  x = (r%2**N )*5**N
  solutions|= {x, 10**N-x, x-1, 10**N-x+1}

  # Green diophantine 
  r,s,g = egcd(5**N, 2**(N-1))

  x=(r%2**(N-1))*5**N
  while x < 10**N:
    solutions|= {x, 10**N-x}
    x+= 5**N*2**(N-1)


  x=2*(r%2**(N-1))*5**N - 1
  while x < 10**N:
    solutions|= {x, 10**N-x}
    x+= 5**N*2**(N-1)

  x=(s%5**N)*2**(N-1) - 1
  while x < 10**N:
    solutions|= {x, 10**N-x}
    x+= 5**N*2**(N-1)

  # Red diophantine 
  r,s,g = egcd(5**N, 2**(N-2))

  x = 2*(r%2**(N-2) )*5**N - 1
  while x < 10**N:
    solutions|= {x, 10**N-x}
    x+= 5**N*2**(N-1)

  x = (s%5**N )*2**(N-1) - 1
  while x < 10**N:
    solutions|= {x, 10**N-x}
    x+= 5**N*2**(N-1)

  # sanity check
  if len(solutions)<>15:
    print "ERROR, wrong number of solutions",len(solutions)
  for a in solutions:
    if pow(a,3,10**N) <> a:
      print "ERROR, Invalid value", a
      
  print "N=",N,":", len(solutions),
  print "solutions, highest non-trivial solution =", max( solutions-{10**N-1})