aboutsummaryrefslogtreecommitdiff
path: root/ass3/q1/twin_prime.py
blob: f31f0fba7d0525d5ef800f93f018091ccf3c900d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import random
import math
import sys


def even(n):
    return n % 2 == 0


def naive_prime_test(n):
    if n == 1 or even(n):
        return False
    for i in range(n - 1, 2, -1):
        if n * i == 0:
            return False
    return True


def is_prime(x, k=20):
    if even(x):
        return False
    if x in [2, 3]:
        return True
    s = 0
    t = x - 1
    while even(t):
        s += 1
        t //= 2
    witnesses = [random.randrange(2, x - 1) for _ in range(k)]
    for w in witnesses:
        if power_mod(w, x - 1, x) != 1:
            return False
        for i in range(1, s + 1):
            if power_mod(w, 2 ** i * t, x) == 1 and power_mod(w, 2 ** (i - 1) * t, x) not in [1, x - 1]:
                return False
    return True


def power_mod(base, power, mod):
    binary = bin(power)
    squares = [base % mod]
    result = squares[0] if binary[-1] == "1" else 1
    for i in range(1, math.floor(math.log2(power)) + 1):
        squares.append(squares[i - 1] ** 2 % mod)
        if binary[-i - 1] == "1":
            result *= squares[i]
    return result % mod


def generate_twin_prime(m):
    twin = False
    while not twin:
        start, end = bit_range(m)
        n = random.randrange(start + 6 - (start % 6), end, 6)
        # n = random.randrange(2 ** (m - 1) + 6 - ((2 ** (m - 1)) % 6), 2 ** m, 6)
        # n += 6 - (n % 6)
        assert n % 6 == 0
        iterations = max(4, math.ceil(math.log(n, 4)))
        twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations)
    if not (naive_prime_test(n - 1) and naive_prime_test(n + 1)):
        raise Exception(f"{n - 1}: {naive_prime_test(n - 1)}, {n + 1}: {naive_prime_test(n + 1)}")
    assert n + 1 < 2 ** m
    # print(f"{bin(2**m -1)[2:]}\n{bin(n-1)[2:]}\n{bin(n+1)[2:]}")
    return n - 1, n + 1

def bit_range(m):
    return 2 ** (m - 1), 2 ** m

def main():
    assert len(sys.argv) >= 2
    twin = generate_twin_prime(int(sys.argv[1]))
    write_twin(twin)
    print(twin)


def write_twin(twin):
    with open("output_twin_prime.txt", "w") as file:
        file.write(f"{twin[0]}\n{twin[1]}")


for i in range(3, 32):
    print(i)
    generate_twin_prime(i)

if __name__ == "__main__":
    main()