144 lines
5.1 KiB
Python
144 lines
5.1 KiB
Python
# https://github.com/dmitrifried/Shors-Algorithm-in-Cirq/blob/master/Shor's%20Algorithm%20in%20Cirq.ipynb
|
||
|
||
from math import gcd,floor,ceil,pow,log2
|
||
from random import randint
|
||
from typing import Optional,Sequence,Union,Callable
|
||
from fractions import Fraction
|
||
|
||
from sympy import isprime
|
||
|
||
from cirq import ArithmeticGate,Circuit,X,H,qft,measure,LineQubit
|
||
from cirq import CircuitDiagramInfoArgs,CircuitDiagramInfo,Result
|
||
from cirq import Simulator
|
||
|
||
def classic_order_finder(x: int, n: int) -> Optional[int]:
|
||
"""
|
||
Поиск делителя классическим способом
|
||
"""
|
||
if x < 2 or n <= x or gcd(x, n) > 1:
|
||
raise ValueError(f'Invalid x={x} for modulus n={n}.')
|
||
r, y = 1, x
|
||
while y != 1:
|
||
y = (x * y) % n
|
||
r += 1
|
||
return r
|
||
|
||
class ModularExp(ArithmeticGate):
|
||
def __init__(self, target: Sequence[int], exponent: Union[int, Sequence[int]], base: int, modulus: int) -> None:
|
||
if len(target) < modulus.bit_length():
|
||
raise ValueError(
|
||
f'Регистр с {len(target)} кубитами маленький для {modulus}'
|
||
)
|
||
self.target = target
|
||
self.exponent = exponent
|
||
self.base = base
|
||
self.modulus = modulus
|
||
|
||
def registers(self) -> Sequence[Union[int, Sequence[int]]]:
|
||
return self.target, self.exponent, self.base, self.modulus
|
||
|
||
def with_registers(self, *new_registers: Union[int, Sequence[int]]) -> 'ModularExp':
|
||
if len(new_registers) != 4:
|
||
raise ValueError(
|
||
f'Expected 4 registers (target, exponent, base, '
|
||
f'modulus), but got {len(new_registers)}'
|
||
)
|
||
target, exponent, base, modulus = new_registers
|
||
if not isinstance(target, Sequence):
|
||
raise ValueError(f'Метка должна быть регистром кубитов {type(target)}')
|
||
if not isinstance(modulus,int):
|
||
raise ValueError(f'модуль это классическая константа {type(modulus)}')
|
||
if not isinstance(base,int):
|
||
raise ValueError(f'база у числа это классическая константа {type(base)}')
|
||
return ModularExp(target, exponent, base, modulus)
|
||
def apply(self, *register_values: int) -> int:
|
||
assert len(register_values) == 4
|
||
target, exponent, base, modulus = register_values
|
||
if target >= modulus:
|
||
return target
|
||
return (target * base**exponent) % modulus
|
||
def _circuit_diagram_info_(self, args: CircuitDiagramInfoArgs) -> CircuitDiagramInfo:
|
||
assert args.known_qubits is not None
|
||
wire_symbols = [f't{i}' for i in range(len(self.target))]
|
||
e_str = str(self.exponent)
|
||
if isinstance(self.exponent, Sequence):
|
||
e_str = 'e'
|
||
wire_symbols += [f'e{i}' for i in range(len(self.exponent))]
|
||
wire_symbols[0] = f'ModularExp(t*{self.base}**{e_str} % {self.modulus})'
|
||
return CircuitDiagramInfo(wire_symbols=tuple(wire_symbols))
|
||
|
||
def make_order_finding_circuit(x: int, n: int) -> Circuit:
|
||
L = n.bit_length()
|
||
target = LineQubit.range(L)
|
||
exponent = LineQubit.range(L, 3 * L + 3)
|
||
return Circuit(
|
||
X(target[L - 1]),
|
||
H.on_each(*exponent),
|
||
ModularExp([2] * len(target), [2] * len(exponent), x, n).on(*target + exponent),
|
||
qft(*exponent, inverse=True),
|
||
measure(*exponent, key='exponent'),
|
||
)
|
||
|
||
def read_eigenphase(result: Result) -> float:
|
||
exponent_as_integer = result.data['exponent'][0]
|
||
exponent_num_bits = result.measurements['exponent'].shape[1]
|
||
return float(exponent_as_integer / 2**exponent_num_bits)
|
||
|
||
def quantum_order_finder(x: int, n: int) -> Optional[int]:
|
||
if x < 2 or n <= x or gcd(x, n) > 1:
|
||
raise ValueError(f'Неправильный x={x} для модуля n={n}.')
|
||
|
||
circuit = make_order_finding_circuit(x, n)
|
||
result = Simulator().run(circuit)
|
||
eigenphase = read_eigenphase(result)
|
||
f = Fraction.from_float(eigenphase).limit_denominator(n)
|
||
if f.numerator == 0:
|
||
return None
|
||
r = f.denominator
|
||
if x**r % n != 1:
|
||
return None
|
||
print(circuit)
|
||
return r
|
||
|
||
def find_factor_of_prime_power(n: int) -> Optional[int]:
|
||
for k in range(2, floor(log2(n)) + 1):
|
||
c = pow(n, 1 / k)
|
||
c1 =floor(c)
|
||
if c1**k == n:
|
||
return c1
|
||
c2 = ceil(c)
|
||
if c2**k == n:
|
||
return c2
|
||
return None
|
||
|
||
def find_factor(
|
||
n: int, order_finder: Callable[[int, int], Optional[int]], max_attempts: int = 30
|
||
) -> Optional[int]:
|
||
if isprime(n):
|
||
return None
|
||
if n % 2 == 0:
|
||
return 2
|
||
c = find_factor_of_prime_power(n)
|
||
if c is not None:
|
||
return c
|
||
for _ in range(max_attempts):
|
||
x = randint(2, n - 1)
|
||
c = gcd(x, n)
|
||
if 1 < c < n:
|
||
return c
|
||
r = order_finder(x, n)
|
||
if r is None:
|
||
continue
|
||
if r % 2 != 0:
|
||
continue
|
||
y = x ** (r // 2) % n
|
||
assert 1 < y < n
|
||
c = gcd(y - 1, n)
|
||
if 1 < c < n:
|
||
return c
|
||
return None
|
||
|
||
n = 33
|
||
d = find_factor(n, quantum_order_finder)
|
||
|