diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | ass3/q1/twin_prime.py | 98 |
2 files changed, 100 insertions, 0 deletions
@@ -1,5 +1,7 @@ *.txt .idea +.vscode +venv # ---> Python # Byte-compiled / optimized / DLL files diff --git a/ass3/q1/twin_prime.py b/ass3/q1/twin_prime.py new file mode 100644 index 0000000..aaeeef8 --- /dev/null +++ b/ass3/q1/twin_prime.py @@ -0,0 +1,98 @@ +import random +import math + + +def even(n): + return n % 2 == 0 + + +def is_really_prime(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 pi2(n): + c = 0.661618 + return 2 * c * (n / (math.log(n) ** 2)) + + +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 + assert math.floor(t) == t + 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 + # return all(map(lambda w: is_prime(x, w), [2, 3, 5, 7, 11, 13, 17])) + + +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 + + +# for i in range(3, 100): +# print(i, is_prime(i), is_really_prime(i)) +# if is_prime(i) != is_really_prime(i): +# print("Liar", i, is_prime(i), is_really_prime(i)) + +# is_really_prime(9) + + +# for a in filter(is_really_prime, range(50000, 500000)): +# if not miller_rabins(a, max(4, math.ceil(math.log(a, 4)))): +# raise Exception(f"Didn't work for {a}") + + +def generate_twin_prime(m): + twin = False + while not twin: + n = random.randrange(2 ** (m - 1), 2 ** m) + iterations = max(4, math.ceil(math.log(n, 4))) + twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations) + if not (is_really_prime(n - 1) and is_really_prime(n + 1)): + raise Exception(f"{n-1}: {is_really_prime(n - 1)}, {n+1}: {is_really_prime(n + 1)}") + return n - 1, n + 1 + + +for bits in range(3, 64): + print(f"==== {bits} bits ==== {2 ** (bits - 1), 2 ** bits - 1}") + for _ in range(1000): + generate_twin_prime(bits) + if _ % 50 == 0: + print(".", end="") + print() + +# i = 0 +# while True: +# i += 1 +# base, power, mod = random.randint(1, 1000), random.randint(1, 5000), random.randint(1, 1000) +# mine = power_mod(base, power, mod) +# real = (base ** power) % mod +# pass_test = mine == real +# if i % 10000 == 0: +# print(i) +# if not pass_test: +# # print(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.") +# raise Exception(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.") |