"""
Required preamble:  R.<x> = PolynomialRing(QQ)

roots(x^3 + x + 1, 7^2) returns
a list of roots of the given polynomial
in GF(49).

search(x^3 + x + 1, 2, 100, 2)
returns lists of roots for the given
polynomial in GF(p^2) for primes between
2 and 100.
"""

def qq2zz(f):
  # clear denominators of f
  c = f.coeffs()
  d = map( lambda f: f.denom(), c)
  return lcm(d)*f

def roots(f, q):
  # return list of roots of f in finite field of q elements
  K.<T> = GF(q)
  r = [ ]
  g = qq2zz(f).change_ring(K)
  for a in K:
    if g(a) == 0:
      r.append(a)
  return r

def search(f,a,b,k):
  # search for roots of f in GF(p^k) for p in [a,b]
  for p in prime_range(a,b):
    rr = roots(f,p^k)
    if rr != [ ]:
      print p, rr