1  def random_mersenne_prime():
 2      while True:
 3          p = random_prime(100, 200)
 4          m = 2^p - 1
 5          if is_prime(m):
 6              return m
 7  
 8  #p = random_mersenne_prime()
 9  p = 8191
10  
11  # -1 should not be a quadratic residue modulo p
12  assert legendre_symbol(-1, p) == -1
13  assert p % 4 == 3
14  
15  F = GF(p)
16  R.<x> = F[]
17  K.<i> = F.extension(x^2 + 1)
18  
19  def get_point(t):
20      x = (1 - t^2)/(1 + t^2)
21      y = 2*t / (1 + t^2)
22      return (x, y)
23  
24  (p1_x, p1_y) = get_point(K(3))
25  assert p1_x^2 + p1_y^2 == 1
26  z = p1_x + i*p1_y
27  
28  (p2_x, p2_y) = get_point(K(7))
29  assert p2_x^2 + p2_y^2 == 1
30  w = p2_x + i*p2_y
31  
32  def abs(z):
33      return z[0]^2 + z[1]^2
34  
35  # z and w are now elements of F_p(i) which is a field
36  assert abs(z) == 1
37  assert abs(w) == 1
38  assert abs(z * w) == 1
39  
40  # We can also construct the inverse
41  def conjugate(z):
42      return z[0] - i*z[1]
43  # Remember that (x + iy)(x - iy) = x^2 - i^2 y^2 = x^2 + y^2
44  assert z * conjugate(z) in F
45  # Now we find the multiplicative inverse
46  z_inv = conjugate(z) / (z * conjugate(z))
47  assert z * z_inv == 1
48  
49  # Size of K is p^2
50  assert len(K) == p^2
51  
52  # Because in sage we cannot construct a homomorphism to GF(p^2) directly
53  # we instead construct the isomorphic field extension, and use that instead.
54  Fp2.<a> = GF(p^2)
55  conway_fp2 = a.minimal_polynomial()
56  L.<j> = F.extension(conway_fp2(x=x))
57  
58  y = L.multiplicative_generator()
59  phi = K.hom([y^(y.multiplicative_order()/i.multiplicative_order())])
60  # K ≌ GF(p^2)
61  assert phi.is_injective() and phi.is_surjective()
62  
63  assert a.multiplicative_order() == p^2 - 1
64  g_K = phi.inverse()(y)
65  g1 = g_K^int((p^2 - 1)/(p + 1))
66  assert abs(g1) == 1
67  
68  g2 = K.multiplicative_generator()
69  g2 = g2^(p - 1)
70  assert abs(g2) == 1
71  
72  # bug in sage where .unit_group is missing
73  # https://ask.sagemath.org/question/62822/make-morphism-from-gfp2s-multiplicative-group-to-gfps-multiplicative-group/
74  C = AbelianGroup([p + 1])
75  g, = C.gens()
76  while True:
77      print(f"Group of order = {g.order()}")
78      for C2 in C.subgroups():
79          print(f"  {C2}")
80      print()
81  
82      if g.order() == 1:
83          break
84  
85      # For some annoying reason, I cannot iterate on subgroups
86      # C = C.subgroup([g^2])
87      # Trying to get the subgroup of this will give me some error.
88      # Lets just construct it manually.
89      C = AbelianGroup([(g^2).order()])
90      g, = C.gens()
91