def pminusone(p,s):
  return Mod(2,p)^s - 1

def ecm(p,s,a):
  x,z = {},{}
  x[1],z[1] = Mod(2,p),Mod(1,p)
  for j in range(2,s+1):
    i = j//2
    if j == 2*i:
      x[j] = (x[i]^2-z[i]^2)^2
      z[j] = 4*x[i]*z[i]*(x[i]^2+a*x[i]*z[i]+z[i]^2)
    else:
      x[j] = 4*(x[i]*x[i+1]-z[i]*z[i+1])^2
      z[j] = 8*(x[i]*z[i+1]-z[i]*x[i+1])^2
  # exercise: do the above in O(log(s)) steps
  return z[s]

s = 60
for p in primes(1,10000):
  if pminusone(p,s) == 0: print('%s divides pminusone(%s)' % (p,s))
  if ecm(p,s,8) == 0: print('%s divides ecm(%s,8)' % (p,s))
  if ecm(p,s,9) == 0: print('%s divides ecm(%s,9)' % (p,s))