Schoof algorithm (without Elkies/Atkins improvements) to compute number of points on a curve

Schoof.sage

from sage.all_cmdline import *

def curve_fx(e: EllipticCurve):
  a4 = e.a4()
  a6 = e.a6()
  X,Y=e.multiplication_by_m(1)
  return X**3 + a4*X + a6  

def trace_mod_2(e: EllipticCurve):
  (X,Y) = e.multiplication_by_m(1)
  q=e.base_field().order()
  if gcd( X**q - X, curve_fx(e)).degree() > 0:
    return 0
  else:
    return 1

def normalize(polyList, fx):
  const_term  = 0
  y_term = 0

  for i, v in enumerate(polyList):
      xterm = polyList[i]
      if i % 2 == 0:
          expo = int(i/2)
          const_term = const_term + xterm * (fx ** expo) 
      else:
          expo = int(int(i-1)/2)
          y_term = y_term + xterm * (fx ** expo)
  return (y_term,const_term)


def mul_iso_mod_l(iso_1, iso_2, div_poly, curve : EllipticCurve):
  x1,y1 = iso_1
  x2,y2 = iso_2
  x3 = x1(x2,y2) % div_poly
  y3 = y1(x2,y2) % div_poly
  return (x3,y3)

def add_iso_mod_l(iso_1, iso_2, div_poly, curve : EllipticCurve):
  X,Y = curve.multiplication_by_m(1)
  (a1,b1) = iso_1
  (a2,b2) = iso_2 

  if a1 == None and b1 == None:
    return iso_2

  if a2 == None and b2 == None:
    return iso_1

  if a1 == a2 and b1 == -b2:
    return (None, None)

  b1 = b1 // Y
  b2 = b2 // Y

  r = None
  fx = curve_fx(curve)

  if b1 == b2:
    r_den = 2*b1*fx
    r_num = 3*a1*a1 + curve.a4()
  else:
    r_num = b1 - b2
    r_den = a1 - a2

  poly_factor = gcd(r_den, div_poly)
  
  if poly_factor.degree() > 0:
    div_poly  = div_poly // poly_factor
    return add_iso_mod_l((a1 % div_poly, b1 % div_poly),
                         (a2 % div_poly, b2 % div_poly),
                         div_poly, curve)  

  r_den_inverse = r_den.univariate_polynomial().inverse_mod(div_poly.univariate_polynomial())

  r = r_num * r_den_inverse 
  a3 = (r*r*fx - a1 - a2) 
  b3 = r*(a1 - a3 ) - b1 

  return (a3 % div_poly, (b3 % div_poly)*Y)


def scalar_mul_iso_mod_l(iso, n : int, div_poly, curve : EllipticCurve):
  X,Y=curve.multiplication_by_m(1)
  n = Integer(n)
  result = (None, None)

  for bit in n.bits()[::-1]:
    result = add_iso_mod_l(result, result, div_poly, curve)
    # Nothing about this is constant time
    if bit == 1:
      result = add_iso_mod_l(result, iso, div_poly, curve)

  return result


def trace_mod_l(e : EllipticCurve, l : int):
  q = e.base_field().characteristic()
  assert gcd(q, l) == 1
  assert l.is_prime()

  if l == 2:
    return trace_mod_2(e)

  curve_x = curve_fx(e)
  X,Y=e.multiplication_by_m(1)
  l_torsion = e.division_polynomial(l,two_torsion_multiplicity=1)
  l_iso = ((X**q) % l_torsion, (((curve_x)**((q-1)//2)) % l_torsion)*Y)
  l_sq_iso = mul_iso_mod_l(l_iso, l_iso, l_torsion, e)
  q_iso = scalar_mul_iso_mod_l((X,Y), q, l_torsion, e)
  target_sum = add_iso_mod_l(q_iso, l_sq_iso, l_torsion, e);

  if target_sum[0] == None and target_sum[1] == None:
    return 0
  
  for i in range(1,l):
    n_time_frob = scalar_mul_iso_mod_l(l_iso, i, l_torsion, e)
    if n_time_frob[0] == target_sum[0] and n_time_frob[1] == target_sum[1]:
      return i

  return 0


def compute_trace(curve : EllipticCurve):
  prime_prods = Integer(1)
  l = Integer(1)
  t = 0
  q = curve.base_field().characteristic()
  search_bound = 4 * int(sqrt(float(q)))

  while prime_prods < search_bound:
    l = next_prime(l)
    if l % q == 0:
      continue

    tl = trace_mod_l(curve, l)
    prime_prods_inv = Integer(prime_prods.inverse_mod(l))
    l_inv           = Integer(l.inverse_mod(prime_prods))

    t = (prime_prods * prime_prods_inv * tl + l * l_inv * t ) % (l*prime_prods)
    prime_prods     = l * prime_prods

  if t > (prime_prods // 2):
    return prime_prods - t 
  
  return t

if __name__=='__main__':
  E=EllipticCurve(GF(11),[1,9]);
  print(trace_mod_2(E)); print(E.trace_of_frobenius() % 2)
  trace = E.trace_of_frobenius()
  trace2 = compute_trace(E)
  print("Trace-1: {}, Trace-2: {}".format(trace, trace2))