In [1]:
a = [1] + [0]*249
for i in range(1, 250251):
n = pow(i, i, 250)
b = [0] * 250
for j in range(250):
b[j] = (a[j] + a[j-n]) % 10**16
a = b
print(a[0] - 1)
In [ ]: