#graphType = "sextet" graphType = "hexagon" #graphType = "triplet" only_girth = True primes = [] upper_bound = 1000000000 writeToFile = True if writeToFile: if only_girth: filename = "girth" else: filename = "full_sequence" outputfile = open(graphType + str(upper_bound) + "_" + filename + ".csv", "w") if only_girth: outputfile.write("p;|G|;Girth\n") else: outputfile.write("p;|G|;Girth;Diameter;Sequence\n") if graphType == "sextet": p = 3 while p * (p * p - 1) <= 48 * upper_bound: f1 = 2 while f1 * f1 <= p and p % f1 > 0: f1 += 1 if p % f1 > 0: p_res = p % 16 if p_res > 8: p_res = 16 - p_res if p_res == 1: g_size = (p * (p * p - 1)) // 48 else: if p_res == 7: g_size = (p * (p * p - 1)) // 24 else: g_size = (p * p * (p ** 4 - 1)) // 24 if g_size <= upper_bound: primes.append([g_size, p]) p += 2 if graphType == "hexagon": p = 5 while p * (p * p - 1) <= 24 * upper_bound: f1 = 2 while f1 * f1 <= p and p % f1 > 0: f1 += 1 if p % f1 > 0: p_res = p % 24 if p_res == 1 or p_res == 23: g_size = (p * (p * p - 1)) // 24 else: g_size = (p * (p * p - 1)) // 12 if g_size <= upper_bound: primes.append([g_size, p]) p += 2 if graphType == "triplet": p = 3 while p * (p * p - 1) <= 12 * upper_bound: f1 = 2 while f1 * f1 <= p and p % f1 > 0: f1 += 1 if p % f1 > 0: p_res = p % 4 if p_res == 1: g_size = (p * (p * p - 1)) // 12 else: g_size = (p * (p * p - 1)) // 6 if g_size <= upper_bound: primes.append([g_size, p]) p += 2 primes.sort() def quartet(x1, x2, x3, x4): x5 = 0 if x1 == p: x5 = (2 * x2 - x3 - x4) % p else: if x2 == p: x5 = (2 * x1 - x3 - x4) % p else: if x3 == p: x5 = (x1 + x2 - 2 * x4) % p else: if x4 == p: x5 = (x1 + x2 - 2 * x3) % p else: x5 = ((x1 - x3) * (x2 - x4) + (x1 - x4) * (x2 - x3)) % p return(x5) def fourth(x1, x2, x3): if x1 < p and x2 < p and x3 < p and (x1 + x2 - 2 * x3) % p == 0: x4 = p else: if x1 == p: x4 = (2 * x2 - x3) % p else: if x2 == p: x4 = (2 * x1 - x3) % p else: if x3 == p: x4 = ((x1 + x2) * p_inv[2]) % p else: x4 = ((2 * x1 * x2 - (x1 + x2) * x3) * p_inv[(x1 + x2 - 2 * x3) % p]) % p return(x4) def quartet2(x1, y1, x2, y2, x3, y3, x4, y4): x5 = 0 y5 = 0 if y1 == p: x5 = (2 * x2 - x3 - x4) % p y5 = (2 * y2 - y3 - y4) % p else: if y2 == p: x5 = (2 * x1 - x3 - x4) % p y5 = (2 * y1 - y3 - y4) % p else: if y3 == p: x5 = (x1 + x2 - 2 * x4) % p y5 = (y1 + y2 - 2 * y4) % p else: if y4 == p: x5 = (x1 + x2 - 2 * x3) % p y5 = (y1 + y2 - 2 * y3) % p else: x5 = ((x1 - x3) * (x2 - x4) - (y1 - y3) * (y2 - y4) + (x1 - x4) * (x2 - x3) - (y1 - y4) * (y2 - y3)) % p y5 = ((x1 - x3) * (y2 - y4) + (y1 - y3) * (x2 - x4) + (x1 - x4) * (y2 - y3) + (x2 - x3) * (y1 - y4)) % p return(x5 + p * y5) def fourth2(x1, y1, x2, y2, x3, y3): if y1 < p and y2 < p and y3 < p and (x1 + x2 - 2 * x3) % p == 0 and (y1 + y2 - 2 * y3) % p == 0: x4 = 0 y4 = p else: if y1 == p: x4 = (2 * x2 - x3) % p y4 = (2 * y2 - y3) % p else: if y2 == p: x4 = (2 * x1 - x3) % p y4 = (2 * y1 - y3) % p else: if y3 == p: x4 = ((x1 + x2) * p_inv[2]) % p y4 = ((y1 + y2) * p_inv[2]) % p else: c = 2 * (x1 * x2 - y1 * y2) - (x1 + x2) * x3 + (y1 + y2) * y3 d = 2 * (x1 * y2 + y1 * x2) - (x1 + x2) * y3 - (y1 + y2) * x3 if (x1 + x2 - 2 * x3) % p == 0: inv_a = 0 inv_b = (-p_inv[(y1 + y2 - 2 * y3) % p]) % p else: inv_a = (p_inv[(x1 + x2 - 2 * x3 + p_inv[(x1 + x2 - 2 * x3) % p] * (y1 + y2 - 2 * y3) ** 2) % p]) inv_b = (p_inv[(x1 + x2 - 2 * x3) % p] * (-(y1 + y2 - 2 * y3) * inv_a)) % p x4 = (c * inv_a - d * inv_b) % p y4 = (c * inv_b + d * inv_a) % p return(x4 + p * y4) def stab1(x1, y1, x2, y2): if x1 == p: a0 = 1 b0 = (-2 * y1) % p c0 = (y1 * (x2 + y2) - x2 * y2) % p else: if y1 == p: a0 = 1 b0 = (-2 * x1) % p c0 = (x1 * (x2 + y2) - x2 * y2) % p else: if x2 == p: a0 = 1 b0 = (-2 * y2) % p c0 = (y2 * (x1 + y1) - x1 * y1) % p else: if y2 == p: a0 = 1 b0 = (-2 * x2) % p c0 = (x2 * (x1 + y1) - x1 * y1) % p else: a0 = (x1 + y1 - x2 - y2) % p b0 = (-2 * (x1 * y1 - x2 * y2)) % p c0 = (x1 * y1 * (x2 + y2) - x2 * y2 * (x1 + y1)) % p if a0 == 0: z1 = (-c0 * p_inv[b0]) % p z2 = p else: z1 = ((-b0 + p_root[(b0 * b0 - 4 * a0 * c0) % p]) * p_inv[(2 * a0) % p]) % p z2 = ((-b0 - p_root[(b0 * b0 - 4 * a0 * c0) % p]) * p_inv[(2 * a0) % p]) % p if z1 > z2: z3 = z1 z1 = z2 z2 = z3 return(z1, z2) def stab2(x1, y1, x2, y2): x1a = x1 % p x1b = x1 // p y1a = y1 % p y1b = y1 // p x2a = x2 % p x2b = x2 // p y2a = y2 % p y2b = y2 // p if x1 == p * p: a0a = 1 a0b = 0 b0a = (-2 * y1a) % p b0b = (-2 * y1b) % p c0a = (y1a * (x2a + y2a) + alpha * y1b * (x2b + y2b) - x2a * y2a - alpha * x2b * y2b) % p c0b = (y1a * (x2b + y2b) + y1b * (x2a + y2a) - x2a * y2b - x2b * y2a) % p else: if y1 == p * p: a0a = 1 a0b = 0 b0a = (-2 * x1a) % p b0b = (-2 * x1b) % p c0a = (x1a * (x2a + y2a) + alpha * x1b * (x2b + y2b) - x2a * y2a - alpha * x2b * y2b) % p c0b = (x1a * (x2b + y2b) + x1b * (x2a + y2a) - x2a * y2b - x2b * y2a) % p else: if x2 == p * p: a0a = 1 a0b = 0 b0a = (-2 * y2a) % p b0b = (-2 * y2b) % p c0a = (y2a * (x1a + y1a) + alpha * y2b * (x1b + y1b) - x1a * y1a - alpha * x1b * y1b) % p c0b = (y2a * (x1b + y1b) + y2b * (x1a + y1a) - x1a * y1b - x1b * y1a) % p else: if y2 == p * p: a0a = 1 a0b = 0 b0a = (-2 * x2a) % p b0b = (-2 * x2b) % p c0a = (x2a * (x1a + y1a) + alpha * x2b * (x1b + y1b) - x1a * y1a - alpha * x1b * y1b) % p c0b = (x2a * (x1b + y1b) + x2b * (x1a + y1a) - x1a * y1b - x1b * y1a) % p else: a0a = (x1a + y1a - x2a - y2a) % p a0b = (x1b + y1b - x2b - y2b) % p b0a = (-2 * (x1a * y1a + alpha * x1b * y1b - x2a * y2a - alpha * x2b * y2b)) % p b0b = (-2 * (x1a * y1b + x1b * y1a - x2a * y2b - x2b * y2a)) % p d0a = (y1a * (x2a + y2a) + alpha * y1b * (x2b + y2b)) d0b = (y1a * (x2b + y2b) + y1b * (x2a + y2a)) e0a = (y2a * (x1a + y1a) + alpha * y2b * (x1b + y1b)) e0b = (y2a * (x1b + y1b) + y2b * (x1a + y1a)) c0a = (x1a * d0a + alpha * x1b * d0b - x2a * e0a - alpha * x2b * e0b) % p c0b = (x1a * d0b + x1b * d0a - x2a * e0b - x2b * e0a) % p if a0a + p * a0b == 0: if b0a == 0: inv_a = 0 inv_b = p_inv[(alpha * b0b) % p] else: inv_a = p_inv[(b0a - alpha * p_inv[b0a] * b0b * b0b) % p] inv_b = (p_inv[b0a] * (-b0b * inv_a)) % p z1a = (-c0a * inv_a - alpha * c0b * inv_b) % p z1b = (-c0a * inv_b - c0b * inv_a) % p z1 = z1a + p * z1b z2 = p * p else: root_0 = p_root[(b0a * b0a + alpha * b0b * b0b - 4 * a0a * c0a - 4 * alpha * a0b * c0b) % p + p * ((2 * b0a * b0b - 4 * a0a * c0b - 4 * a0b * c0a) % p)] root_a = root_0 % p root_b = root_0 // p if a0a == 0: inv_a = 0 inv_b = p_inv[(2 * alpha * a0b) % p] else: inv_a = p_inv[(2 * a0a - alpha * p_inv[(2 * a0a) % p] * 4 * a0b * a0b) % p] inv_b = (p_inv[(2 * a0a) % p] * (-2 * a0b * inv_a)) % p z1a = ((-b0a + root_a) * inv_a + alpha * (-b0b + root_b) * inv_b) % p z1b = ((-b0a + root_a) * inv_b + (-b0b + root_b) * inv_a) % p z2a = ((-b0a - root_a) * inv_a + alpha * (-b0b - root_b) * inv_b) % p z2b = ((-b0a - root_a) * inv_b + (-b0b - root_b) * inv_a) % p z1 = z1a + p * z1b z2 = z2a + p * z2b if z1 > z2: z3 = z1 z1 = z2 z2 = z3 return(z1, z2) for i0 in range(len(primes)): p = primes[i0][1] print("\np =", p) if writeToFile: outputfile.write(str(p) + ";") alpha = (7 - 3 * (p % 4)) // 2 # -1 for p == 3 (mod 4), +2 for p == 5 (mod 8) p_inv = [0 for x in range(p + 1)] p_inv[0] = p for i in range(1, p): j = 1 while (i * j) % p != 1: j += 1 p_inv[i] = j q = [0] a = 0 b = 1 no_loop = True odd_cycle = False depth = 0 #replace by len(spread) or something? girth = 0 spread = [] if graphType == "sextet": if p % 8 == 1: p_root = [-1 for x in range(p)] for i in range(p): j = (i * i) % p if p_root[j] == -1: p_root[j] = i s = [[[0, p], [1, p - 1], [p_root[p - 1], p - p_root[p - 1]]]] else: p_root = [-1 for x in range(p * p)] for j in range(p): for i in range(p): k = (i * i + alpha * j * j) % p + p * ((2 * i * j) % p) if p_root[k] == -1: p_root[k] = i + p * j s = [[[0, p * p], [1, p - 1], [p_root[p - 1], (-p_root[p - 1]) % p + p * ((-(p_root[p - 1] // p)) % p)]]] if p % 8 == 1: while a < b and (girth == 0 or not only_girth): spread.append(b - a) b0 = b i = a while i < b0 and not (only_girth and odd_cycle): for j in range(3): sextet = [s[i][j]] z1, z2 = stab1(s[i][(j + 1) % 3][0], s[i][(j + 2) % 3][0], s[i][(j + 1) % 3][1], s[i][(j + 2) % 3][1]) sextet.append([z1, z2]) z3, z4 = stab1(s[i][(j + 1) % 3][0], s[i][(j + 2) % 3][1], s[i][(j + 1) % 3][1], s[i][(j + 2) % 3][0]) sextet.append([z3, z4]) sextet.sort() c = -1 d = b e = 0 found = False while not found and d - c >= 2: if d - c == 1: e = c + d - e else: e = (c + d) // 2 if e == -1 or e == b: e = c + d - e f = 0 while f < 6 and sextet[f // 2][f % 2] == s[q[e]][f // 2][f % 2]: f += 1 if f == 6: found = True else: if sextet[f // 2][f % 2] > s[q[e]][f // 2][f % 2]: c = e else: d = e if not found: s.append(sextet) if sextet[f // 2][f % 2] < s[q[e]][f // 2][f % 2]: q.insert(e, b) else: if e == b - 1: q.append(b) else: q.insert(e + 1, b) b += 1 else: if no_loop and q[e] >= a and q[e] < b0: odd_cycle = True i += 1 if no_loop and b - b0 < 2 * (b0 - a): no_loop = False if odd_cycle: girth = 2 * depth + 1 else: girth = 2 * depth + 2 a = b0 depth += 1 else: while a < b and (girth == 0 or not only_girth): spread.append(b - a) b0 = b i = a while i < b0 and not (only_girth and odd_cycle): for j in range(3): sextet = [s[i][j]] z1, z2 = stab2(s[i][(j + 1) % 3][0], s[i][(j + 2) % 3][0], s[i][(j + 1) % 3][1], s[i][(j + 2) % 3][1]) sextet.append([z1, z2]) z3, z4 = stab2(s[i][(j + 1) % 3][0], s[i][(j + 2) % 3][1], s[i][(j + 1) % 3][1], s[i][(j + 2) % 3][0]) sextet.append([z3, z4]) sextet.sort() c = -1 d = b e = 0 found = False while not found and d - c >= 2: if d - c == 1: e = c + d - e else: e = (c + d) // 2 if e == -1 or e == b: e = c + d - e f = 0 while f < 6 and sextet[f // 2][f % 2] == s[q[e]][f // 2][f % 2]: f += 1 if f == 6: found = True else: if sextet[f // 2][f % 2] > s[q[e]][f // 2][f % 2]: c = e else: d = e if not found: s.append(sextet) if sextet[f // 2][f % 2] < s[q[e]][f // 2][f % 2]: q.insert(e, b) else: if e == b - 1: q.append(b) else: q.insert(e + 1, b) b += 1 else: if no_loop and q[e] >= a and q[e] < b0: odd_cycle = True i += 1 if no_loop and b - b0 < 2 * (b0 - a): no_loop = False if odd_cycle: girth = 2 * depth + 1 else: girth = 2 * depth + 2 a = b0 depth += 1 if graphType == "hexagon": if p % 4 == 1: p_root = [-1 for x in range(p)] for i in range(p): j = (i * i) % p if p_root[j] == -1: p_root[j] = i else: p_root = [-1 for x in range(p * p)] for j in range(p): for i in range(p): k = (i * i + alpha * j * j) % p + p * ((2 * i * j) % p) if p_root[k] == -1: p_root[k] = i + p * j if p % 4 == 1: h = [[0, 1, 3, p, p - 3, p - 1]] else: h = [[0, 1, 3, p * p, p - 3, p - 1]] if p % 4 == 1: while a < b and (girth == 0 or not only_girth): spread.append(b - a) b0 = b i = a while i < b0 and not (only_girth and odd_cycle): for j in range(3): z1 = h[i][j] z2 = h[i][j + 3] z3 = h[i][j + 1] z4 = h[i][(j + 5) % 6] z5 = h[i][j + 2] z6 = h[i][(j + 4) % 6] z7, z8 = stab1(z1, z2, z3, z4) z9, z0 = stab1(z1, z2, z5, z6) if quartet(z7, z0, z1, z9) > 0 or quartet(z7, z0, z2, z8) > 0 or quartet(z8, z9, z1, z0) > 0 or quartet(z8, z9, z2, z7) > 0: z_spare = z9 z9 = z0 z0 = z_spare z_min = min(z1, z2, z7, z8, z9, z0) while z1 != z_min: z_spare = z1 z1 = z7 z7 = z9 z9 = z2 z2 = z0 z0 = z8 z8 = z_spare if z7 > z8: z_spare = z7 z7 = z8 z8 = z_spare z_spare = z9 z9 = z0 z0 = z_spare hexagon = [z1, z7, z9, z2, z0, z8] c = -1 d = b e = 0 found = False while not found and d - c >= 2: if d - c == 1: e = c + d - e else: e = (c + d) // 2 if e == -1 or e == b: e = c + d - e f = 0 while f < 6 and hexagon[f] == h[q[e]][f]: f += 1 if f == 6: found = True else: if hexagon[f] > h[q[e]][f]: c = e else: d = e if not found: h.append(hexagon) if hexagon[f] < h[q[e]][f]: q.insert(e, b) else: if e == b - 1: q.append(b) else: q.insert(e + 1, b) b += 1 else: if no_loop and q[e] >= a and q[e] < b0: odd_cycle = True i += 1 if no_loop and b - b0 < 2 * (b0 - a): no_loop = False if odd_cycle: girth = 2 * depth + 1 else: girth = 2 * depth + 2 a = b0 depth += 1 else: while a < b and (girth == 0 or not only_girth): spread.append(b - a) b0 = b i = a while i < b0 and not (only_girth and odd_cycle): for j in range(3): z1 = h[i][j] z2 = h[i][j + 3] z3 = h[i][j + 1] z4 = h[i][(j + 5) % 6] z5 = h[i][j + 2] z6 = h[i][(j + 4) % 6] z7, z8 = stab2(z1, z2, z3, z4) z9, z0 = stab2(z1, z2, z5, z6) if quartet2(z7 % p, z7 // p, z0 % p, z0 // p, z1 % p, z1 // p, z9 % p, z9 // p) > 0 or quartet2(z7 % p, z7 // p, z0 % p, z0 // p, z2 % p, z2 // p, z8 % p, z8 // p) > 0 or quartet2(z8 % p, z8 // p, z9 % p, z9 // p, z1 % p, z1 // p, z0 % p, z0 // p) > 0 or quartet2(z8 % p, z8 // p, z9 % p, z9 // p, z2 % p, z2 // p, z7 % p, z7 // p) > 0: z_spare = z9 z9 = z0 z0 = z_spare z_min = min(z1, z2, z7, z8, z9, z0) while z1 != z_min: z_spare = z1 z1 = z7 z7 = z9 z9 = z2 z2 = z0 z0 = z8 z8 = z_spare if z7 > z8: z_spare = z7 z7 = z8 z8 = z_spare z_spare = z9 z9 = z0 z0 = z_spare hexagon = [z1, z7, z9, z2, z0, z8] c = -1 d = b e = 0 found = False while not found and d - c >= 2: if d - c == 1: e = c + d - e else: e = (c + d) // 2 if e == -1 or e == b: e = c + d - e f = 0 while f < 6 and hexagon[f] == h[q[e]][f]: f += 1 if f == 6: found = True else: if hexagon[f] > h[q[e]][f]: c = e else: d = e if not found: h.append(hexagon) if hexagon[f] < h[q[e]][f]: q.insert(e, b) else: if e == b - 1: q.append(b) else: q.insert(e + 1, b) b += 1 else: if no_loop and q[e] >= a and q[e] < b0: odd_cycle = True i += 1 if no_loop and b - b0 < 2 * (b0 - a): no_loop = False if odd_cycle: girth = 2 * depth + 1 else: girth = 2 * depth + 2 a = b0 depth += 1 if graphType == "triplet": t = [[0, 1, 2]] while a < b and (girth == 0 or not only_girth): spread.append(b - a) b0 = b i = a while i < b0 and not (only_girth and odd_cycle): for j in range(3): z1 = t[i][(j + 1) % 3] z2 = t[i][(j + 2) % 3] z3 = t[i][j] z4 = fourth(z1, z2, z3) triplet = [z1, z2, z4] triplet.sort() c = -1 d = b e = 0 found = False while not found and d - c >= 2: if d - c == 1: e = c + d - e else: e = (c + d) // 2 if e == -1 or e == b: e = c + d - e f = 0 while f < 3 and triplet[f] == t[q[e]][f]: f += 1 if f == 3: found = True else: if triplet[f] > t[q[e]][f]: c = e else: d = e if not found: t.append(triplet) if triplet[f] < t[q[e]][f]: q.insert(e, b) else: if e == b - 1: q.append(b) else: q.insert(e + 1, b) b += 1 else: if no_loop and q[e] >= a and q[e] < b0: odd_cycle = True i += 1 if no_loop and b - b0 < 2 * (b0 - a): no_loop = False if odd_cycle: girth = 2 * depth + 1 else: girth = 2 * depth + 2 a = b0 depth += 1 if only_girth: print("Order =", primes[i0][0]) if writeToFile: outputfile.write(str(primes[i0][0])) else: print("Order =", len(q)) if writeToFile: outputfile.write(str(len(q))) print("Girth =", girth) if writeToFile: outputfile.write(";" + str(girth)) if not only_girth: print("Diameter =", depth - 1) if writeToFile: outputfile.write(";" + str(depth - 1)) for j in range(depth): print(spread[j], end = " ") if writeToFile: outputfile.write(";" + str(spread[j])) print("") if writeToFile: outputfile.write("\n") if writeToFile: outputfile.close()