|
| 1 | +""" |
| 2 | +Classical simulation of Shor's Algorithm to factor integers. |
| 3 | +
|
| 4 | +Source: https://en.wikipedia.org/wiki/Shor%27s_algorithm |
| 5 | +""" |
| 6 | + |
| 7 | +import random |
| 8 | +import math |
| 9 | +from typing import Tuple, Union |
| 10 | + |
| 11 | + |
| 12 | +def is_prime(n: int) -> bool: |
| 13 | + """ |
| 14 | + Check if a number is prime. |
| 15 | +
|
| 16 | + >>> is_prime(2) |
| 17 | + True |
| 18 | + >>> is_prime(4) |
| 19 | + False |
| 20 | + """ |
| 21 | + if n < 2: |
| 22 | + return False |
| 23 | + if n in (2, 3): |
| 24 | + return True |
| 25 | + if n % 2 == 0: |
| 26 | + return False |
| 27 | + r = int(math.isqrt(n)) |
| 28 | + for i in range(3, r + 1, 2): |
| 29 | + if n % i == 0: |
| 30 | + return False |
| 31 | + return True |
| 32 | + |
| 33 | + |
| 34 | +def modexp(a: int, b: int, m: int) -> int: |
| 35 | + """ |
| 36 | + Modular exponentiation: (a^b) % m |
| 37 | +
|
| 38 | + >>> modexp(2, 5, 13) |
| 39 | + 6 |
| 40 | + """ |
| 41 | + result = 1 |
| 42 | + a = a % m |
| 43 | + while b > 0: |
| 44 | + if b & 1: |
| 45 | + result = (result * a) % m |
| 46 | + a = (a * a) % m |
| 47 | + b >>= 1 |
| 48 | + return result |
| 49 | + |
| 50 | + |
| 51 | +def shor_classical(N: int, max_attempts: int = 10) -> Union[str, Tuple[int, int]]: |
| 52 | + """ |
| 53 | + Classical approximation of Shor's Algorithm to factor a number. |
| 54 | +
|
| 55 | + >>> result = shor_classical(15) |
| 56 | + >>> isinstance(result, tuple) |
| 57 | + True |
| 58 | + >>> sorted(result) == [3, 5] |
| 59 | + True |
| 60 | +
|
| 61 | + >>> shor_classical(13) # Prime |
| 62 | + 'No factors: 13 is prime' |
| 63 | + """ |
| 64 | + if N <= 1: |
| 65 | + return "Failure: input must be > 1" |
| 66 | + if N % 2 == 0: |
| 67 | + return 2, N // 2 |
| 68 | + if is_prime(N): |
| 69 | + return f"No factors: {N} is prime" |
| 70 | + |
| 71 | + for _ in range(max_attempts): |
| 72 | + a = random.randrange(2, N - 1) |
| 73 | + g = math.gcd(a, N) |
| 74 | + if g > 1: |
| 75 | + return g, N // g |
| 76 | + |
| 77 | + r = 1 |
| 78 | + while r < N: |
| 79 | + if modexp(a, r, N) == 1: |
| 80 | + break |
| 81 | + r += 1 |
| 82 | + else: |
| 83 | + continue |
| 84 | + |
| 85 | + if r % 2 != 0: |
| 86 | + continue |
| 87 | + x = modexp(a, r // 2, N) |
| 88 | + if x == N - 1: |
| 89 | + continue |
| 90 | + |
| 91 | + factor1 = math.gcd(x - 1, N) |
| 92 | + factor2 = math.gcd(x + 1, N) |
| 93 | + if factor1 not in (1, N) and factor2 not in (1, N): |
| 94 | + return factor1, factor2 |
| 95 | + |
| 96 | + return "Failure: try more attempts" |
| 97 | + |
| 98 | + |
| 99 | +if __name__ == "__main__": |
| 100 | + import doctest |
| 101 | + |
| 102 | + doctest.testmod() |
0 commit comments