for q in range(40):
  q = Integer(q)
  if not q.is_prime(): continue
  if q % 4 != 1: continue

  Fq.<alpha> = GF(q)
  def chi(a): return a^((q-1)//2)

  for A in Fq:
    for B in Fq:
      if A*B*(A^2-4*B) == 0: continue
      if (A^2-4*B).is_square(): continue
      for u in Fq:
        if u.is_square(): continue

        def psi(r):
          if r == 0: return (0,0)
          v = -A/(1+u*r^2)
          eps = chi(v^3+A*v^2+B*v)
          x = eps*v-(1-eps)*A/2
          y = -eps*sqrt(x^3+A*x^2+B*x)
          return (x,y)

	b = 0
	while 2^(b+1) <= q: b += 1

	def sigma(r):
	  return Fq(sum([r[i]*2^i for i in range(b)]))

	S = []
	for j in range(2^b):
	  r = [(j>>i)%2 for i in range(b)]
	  if sigma(r) <= (q-1)//2:
	    S.append(r)

	print len(S) == (q+1)//2

	def iota(r): return psi(sigma(r))

	for r in S:
	  x,y = iota(r)
	  print y^2 == x^3 + A*x^2 + B*x

	iotaS = set()
	for r in S:
	  iotaS.add(iota(r))

	print len(iotaS) == len(S)

        psiFq = set()
        for r in Fq:
          psiFq.add(psi(r))

	print iotaS == psiFq