больше кодовввв

task10/grover.py

@@ -0,0 +1,118 @@
+# pylint: disable=wrong-or-nonexistent-copyright-notice
+"""Demonstrates Grover algorithm.
+
+The Grover algorithm takes a black-box oracle implementing a function
+{f(x) = 1 if x==x', f(x) = 0 if x!= x'} and finds x' within a randomly
+ordered sequence of N items using O(sqrt(N)) operations and O(N log(N)) gates,
+with the probability p >= 2/3.
+
+At the moment, only 2-bit sequences (for which one pass through Grover operator
+is enough) are considered.
+
+=== REFERENCE ===
+Coles, Eidenbenz et al. Quantum Algorithm Implementations for Beginners
+https://arxiv.org/abs/1804.03719
+
+=== EXAMPLE OUTPUT ===
+Secret bit sequence: [1, 0]
+Circuit:
+(0, 0): ───────H───────@───────H───X───────@───────X───H───M───
+                       │                   │               │
+(1, 0): ───────H───X───@───X───H───X───H───X───H───X───H───M───
+                       │
+(2, 0): ───X───H───────X───────────────────────────────────────
+Sampled results:
+Counter({'10': 10})
+Most common bitstring: 10
+Found a match: True
+
+"""
+
+from __future__ import annotations
+
+import random
+
+import cirq
+
+
+def set_io_qubits(qubit_count):
+    """Add the specified number of input and output qubits."""
+    input_qubits = [cirq.GridQubit(i, 0) for i in range(qubit_count)]
+    output_qubit = cirq.GridQubit(qubit_count, 0)
+    return (input_qubits, output_qubit)
+
+
+def make_oracle(input_qubits, output_qubit, x_bits):
+    """Implement function {f(x) = 1 if x==x', f(x) = 0 if x!= x'}."""
+    # Make oracle.
+    # for (1, 1) it's just a Toffoli gate
+    # otherwise negate the zero-bits.
+    yield (cirq.X(q) for (q, bit) in zip(input_qubits, x_bits) if not bit)
+    yield (cirq.TOFFOLI(input_qubits[0], input_qubits[1], output_qubit))
+    yield (cirq.X(q) for (q, bit) in zip(input_qubits, x_bits) if not bit)
+
+
+def make_grover_circuit(input_qubits, output_qubit, oracle):
+    """Find the value recognized by the oracle in sqrt(N) attempts."""
+    # For 2 input qubits, that means using Grover operator only once.
+    c = cirq.Circuit()
+
+    # Initialize qubits.
+    c.append([cirq.X(output_qubit), cirq.H(output_qubit), cirq.H.on_each(*input_qubits)])
+
+    # Query oracle.
+    c.append(oracle)
+
+    # Construct Grover operator.
+    c.append(cirq.H.on_each(*input_qubits))
+    c.append(cirq.X.on_each(*input_qubits))
+    c.append(cirq.H.on(input_qubits[1]))
+    c.append(cirq.CNOT(input_qubits[0], input_qubits[1]))
+    c.append(cirq.H.on(input_qubits[1]))
+    c.append(cirq.X.on_each(*input_qubits))
+    c.append(cirq.H.on_each(*input_qubits))
+
+    # Measure the result.
+    c.append(cirq.measure(*input_qubits, key='result'))
+
+    return c
+
+
+def bitstring(bits):
+    return ''.join(str(int(b)) for b in bits)
+
+
+def main():
+    qubit_count = 2
+    circuit_sample_count = 10
+
+    # Set up input and output qubits.
+    (input_qubits, output_qubit) = set_io_qubits(qubit_count)
+
+    # Choose the x' and make an oracle which can recognize it.
+    x_bits = [random.randint(0, 1) for _ in range(qubit_count)]
+    print(f'Secret bit sequence: {x_bits}')
+
+    # Make oracle (black box)
+    oracle = make_oracle(input_qubits, output_qubit, x_bits)
+
+    # Embed the oracle into a quantum circuit implementing Grover's algorithm.
+    circuit = make_grover_circuit(input_qubits, output_qubit, oracle)
+    print('Circuit:')
+    print(circuit)
+
+    # Sample from the circuit a couple times.
+    simulator = cirq.Simulator()
+    result = simulator.run(circuit, repetitions=circuit_sample_count)
+
+    frequencies = result.histogram(key='result', fold_func=bitstring)
+    print(f'Sampled results:\n{frequencies}')
+
+    # Check if we actually found the secret value.
+    most_common_bitstring = frequencies.most_common(1)[0][0]
+    print(f'Most common bitstring: {most_common_bitstring}')
+    print(f'Found a match: {most_common_bitstring == bitstring(x_bits)}')
+
+
+if __name__ == '__main__':
+    main()

task2/main.py

@@ -0,0 +1,25 @@
+import cirq
+
+q0, q1 = cirq.LineQubit.range(2)
+
+oracles = {
+    "0": [],
+    "1": [cirq.X(q1)],
+    "x": [cirq.CNOT(q0, q1)],
+    "notx": [cirq.CNOT(q0, q1), cirq.X(q1)],
+}
+
+
+def deutsch_algorithm(oracle):
+    yield cirq.X(q1)
+    yield cirq.H(q0), cirq.H(q1)
+    yield oracle
+    yield cirq.H(q0)
+    yield cirq.measure(q0)
+
+
+simulator = cirq.Simulator()
+
+for key, oracle in oracles.items():
+    result = simulator.run(cirq.Circuit(deutsch_algorithm(oracle)), repetitions=10)
+    print("oracle: {:<4} results: {}".format(key, result))

task3/main.py

@@ -0,0 +1,61 @@
+import cirq
+
+# Кубиты: два входных (q0, q1) и целевой вспомогательный (q2)
+q0, q1, q2 = cirq.LineQubit.range(3)
+
+# Оракулы для константных функций: f(x)=0 и f(x)=1
+constant = (
+    [],  # ничего не делаем с целевым
+    [cirq.X(q2)],  # инвертируем целевой
+)
+
+# Оракулы для сбалансированных функций
+balanced = (
+    [cirq.CNOT(q0, q2)],  # f=x0
+    [cirq.CNOT(q1, q2)],  # f=x1
+    [cirq.CNOT(q0, q2), cirq.CNOT(q1, q2)],
+    [cirq.CNOT(q0, q2), cirq.X(q2)],
+    [cirq.CNOT(q1, q2), cirq.X(q2)],
+    [cirq.CNOT(q0, q2), cirq.CNOT(q1, q2), cirq.X(q2)],
+)
+
+
+def deutsch_jozsa_circuit(oracle):
+    """Схема DJ через yield: измеряем только q0 и q1."""
+    # 1) Входы в равномерную суперпозицию
+    yield cirq.H(q0), cirq.H(q1)
+
+    # 2) Подготовка целевого в |-> для фазового трюка
+    yield cirq.X(q2), cirq.H(q2)
+
+    # 3) Вызов оракула U_f
+    yield oracle
+
+    # 4) Интерференция на входах
+    yield cirq.H(q0), cirq.H(q1)
+
+    # 5) Измеряем только входные кубиты
+    yield cirq.measure(q0, q1, key="in")
+
+
+def interpret(bits):
+    # bits — это массив вида [q0, q1]
+    return "constant" if (bits[0] == 0 and bits[1] == 0) else "balanced"
+
+
+# Симулятор
+sim = cirq.Simulator()
+
+print("Результаты для константных функций:")
+for i, oracle in enumerate(constant):
+    result = sim.run(cirq.Circuit(deutsch_jozsa_circuit(oracle)), repetitions=1)
+    bits = result.measurements["in"][0]
+    print(f"  Константный оракул {i}: meas(q0,q1)={bits.tolist()} -> {interpret(bits)}")
+
+print("\nРезультаты для сбалансированных функций:")
+for i, oracle in enumerate(balanced):
+    result = sim.run(cirq.Circuit(deutsch_jozsa_circuit(oracle)), repetitions=1)
+    bits = result.measurements["in"][0]
+    print(
+        f"  Сбалансированный оракул {i}: meas(q0,q1)={bits.tolist()} -> {interpret(bits)}"
+    )

task4/main.py

@@ -0,0 +1,91 @@
+"""Алгоритм Бернштейна–Вазирани в Cirq."""
+
+import random
+
+import cirq
+
+
+def main():
+    """Выполняет алгоритм BV."""
+    # Количество кубитов
+    qubit_count = 8
+    # Количество выборок из схемы
+    circuit_sample_count = 3
+
+    # Выберите кубиты для использования
+    input_qubits = [cirq.GridQubit(i, 0) for i in range(qubit_count)]
+    output_qubit = cirq.GridQubit(qubit_count, 0)
+
+    # Выберите коэффициенты для оракула и создайте схему для его запроса
+    secret_bias_bit = random.randint(0, 1)
+    secret_factor_bits = [random.randint(0, 1) for _ in range(qubit_count)]
+
+    oracle = make_oracle(
+        input_qubits, output_qubit, secret_factor_bits, secret_bias_bit
+    )
+
+    print(
+        "Secret function:\nf(x) = x*<{}> + {} (mod 2)".format(
+            ", ".join(str(e) for e in secret_factor_bits), secret_bias_bit
+        )
+    )
+
+    # Встраиваем оракул в специальную квантовую схему, запрашивая ее ровно один раз
+    circuit = make_bernstein_vazirani_circuit(input_qubits, output_qubit, oracle)
+    print("\nCircuit:")
+    print(circuit)
+
+    # Выберем из схемы пару раз
+    simulator = cirq.Simulator()
+    result = simulator.run(circuit, repetitions=circuit_sample_count)
+    frequencies = result.histogram(key="result", fold_func=bitstring)
+    print("\nSampled results:\n{}".format(frequencies))
+
+    # Проверить, действительно ли мы нашли секретное значение
+    most_common_bitstring = frequencies.most_common(1)[0][0]
+    print(
+        "\nMost common matches secret factors:\n{}".format(
+            most_common_bitstring == bitstring(secret_factor_bits)
+        )
+    )
+
+
+def make_oracle(input_qubits, output_qubit, secret_factor_bits, secret_bias_bit):
+    """Вентили, реализующие функцию f(a) = a*factors + bias (mod 2)."""
+    if secret_bias_bit:
+        yield cirq.X(output_qubit)
+
+    for qubit, bit in zip(input_qubits, secret_factor_bits):
+        if bit:
+            yield cirq.CNOT(qubit, output_qubit)
+
+
+def make_bernstein_vazirani_circuit(input_qubits, output_qubit, oracle):
+    """Находит множители в f (a) = a * factors + bias (mod 2) с помощью одного запроса."""
+    c = cirq.Circuit()
+
+    # Инициализировать кубиты
+    c.append(
+        [
+            cirq.X(output_qubit),
+            cirq.H(output_qubit),
+            cirq.H.on_each(*input_qubits),
+        ]
+    )
+
+    # Запрос оракула
+    c.append(oracle)
+
+    # Измерение по оси X
+    c.append([cirq.H.on_each(*input_qubits), cirq.measure(*input_qubits, key="result")])
+
+    return c
+
+
+def bitstring(bits):
+    """Создает битовую строку из итерации битов."""
+    return "".join(str(int(b)) for b in bits)
+
+
+if __name__ == "__main__":
+    main()

task5/main.py

@@ -0,0 +1,104 @@
+from __future__ import annotations
+
+from collections import Counter
+
+import cirq
+import numpy as np
+import scipy as sp
+
+
+def main(qubit_count=3):
+
+    data = []  # we'll store here the results
+
+    # define a secret string:
+    secret_string = np.random.randint(2, size=qubit_count)
+
+    print(f"Secret string = {secret_string}")
+
+    n_samples = 100
+    for _ in range(n_samples):
+        flag = False  # check if we have a linearly independent set of measures
+        while not flag:
+            # Choose qubits to use.
+            input_qubits = [cirq.GridQubit(i, 0) for i in range(qubit_count)]  # input x
+            output_qubits = [
+                cirq.GridQubit(i + qubit_count, 0) for i in range(qubit_count)
+            ]  # output f(x)
+
+            # Pick coefficients for the oracle and create a circuit to query it.
+            oracle = make_oracle(input_qubits, output_qubits, secret_string)
+
+            # Embed oracle into special quantum circuit querying it exactly once
+            circuit = make_simon_circuit(input_qubits, output_qubits, oracle)
+
+            # Sample from the circuit a n-1 times (n = qubit_count).
+            simulator = cirq.Simulator()
+            results = [
+                simulator.run(circuit).measurements["result"][0]
+                for _ in range(qubit_count - 1)
+            ]
+
+            # Classical Post-Processing:
+            flag = post_processing(data, results)
+
+    freqs = Counter(data)
+    print("Circuit:")
+    print(circuit)
+    print(f"Most common answer was : {freqs.most_common(1)[0]}")
+
+
+def make_oracle(input_qubits, output_qubits, secret_string):
+    """Gates implementing the function f(a) = f(b) iff a ⨁ b = s"""
+    # Copy contents to output qubits:
+    for control_qubit, target_qubit in zip(input_qubits, output_qubits):
+        yield cirq.CNOT(control_qubit, target_qubit)
+
+    # Create mapping:
+    if sum(secret_string):  # check if the secret string is non-zero
+        # Find significant bit of secret string (first non-zero bit)
+        significant = list(secret_string).index(1)
+
+        # Add secret string to input according to the significant bit:
+        for j in range(len(secret_string)):
+            if secret_string[j] > 0:
+                yield cirq.CNOT(input_qubits[significant], output_qubits[j])
+
+
+def make_simon_circuit(input_qubits, output_qubits, oracle):
+    """Solves for the secret period s of a 2-to-1 function such that
+    f(x) = f(y) iff x ⨁ y = s
+    """
+
+    c = cirq.Circuit()
+
+    # Initialize qubits.
+    c.append([cirq.H.on_each(*input_qubits)])
+
+    # Query oracle.
+    c.append(oracle)
+
+    # Measure in X basis.
+    c.append([cirq.H.on_each(*input_qubits), cirq.measure(*input_qubits, key="result")])
+
+    return c
+
+
+def post_processing(data, results):
+    """Solves a system of equations with modulo 2 numbers"""
+    sing_values = sp.linalg.svdvals(results)
+    tolerance = 1e-5
+    if (
+        sum(sing_values < tolerance) == 0
+    ):  # check if measurements are linearly dependent
+        flag = True
+        null_space = sp.linalg.null_space(results).T[0]
+        solution = np.around(null_space, 3)  # chop very small values
+        minval = abs(min(solution[np.nonzero(solution)], key=abs))
+        solution = (solution / minval % 2).astype(int)  # renormalize vector mod 2
+        data.append(str(solution))
+        return flag
+
+
+if __name__ == "__main__":
+    main()

task6/main.py

@@ -0,0 +1,335 @@
+import cirq
+import numpy as np
+
+
+def bitstring(bits):
+    """Преобразует список битов в строку, состоящую из символов '0' и '1'.
+
+    Аргументы:
+    bits -- список или массив битов (0 или 1)
+
+    Возвращает:
+    строку, например, '0101101'
+    """
+    return "".join(str(int(b)) for b in bits)
+
+
+class BB84(object):
+
+    def __init__(self, repetitions=1, random_seed=None):
+        """
+        Инициализация протокола BB84.
+
+        repetitions -- число повторений измерений (обычно 1, для статистики можно больше)
+        random_seed -- зерно генератора случайных чисел для воспроизводимости (если None, сгенерируется случайно)
+        """
+        super().__init__()
+        if random_seed is None:
+            random_seed = np.random.randint(0, 2**31 - 1)  # Генерируем случайное зерно
+        np.random.seed(random_seed)  # Настраиваем генератор случайных чисел
+        self.repetitions = repetitions  # Количество повторов симуляции схемы
+        self.circuits = {}  # Словарь для хранения схем с разными названиями
+
+    def setup_demo(self, num_qubits):
+        """
+        Генерация случайных параметров для одной сессии BB84.
+
+        Выбираются:
+        - Базисы Алисы (0 - вычислительный, 1 - базис Адамара)
+        - Состояния Алисы для каждого кубита (0 или 1)
+        - Базисы Боба (случайные 0 или 1)
+        - Базисы Евы (случайные 0 или 1)
+        Также вычисляется ожидаемый секретный ключ для совпадающих базисов Алисы и Боба
+        """
+        self.alice_basis = [
+            np.random.randint(0, 2) for _ in range(num_qubits)
+        ]  # 0 или 1
+        self.alice_state = [
+            np.random.randint(0, 2) for _ in range(num_qubits)
+        ]  # 0 или 1
+        self.bob_basis = [np.random.randint(0, 2) for _ in range(num_qubits)]
+        self.eve_basis = [np.random.randint(0, 2) for _ in range(num_qubits)]
+
+        # Формируем ожидаемый ключ - отбираем биты Алисы там, где базисы совпали с Бобом
+        self.expected_key = bitstring(
+            [
+                self.alice_state[i]
+                for i in range(num_qubits)
+                if self.alice_basis[i] == self.bob_basis[i]
+            ]
+        )
+
+    def make_bb84_circuit(
+        self, num_qubits, sender_basis, receiver_basis, sender_state, name="1"
+    ):
+        """
+        Построение квантовой схемы для протокола BB84.
+
+        Параметры:
+        num_qubits -- число кубитов
+        sender_basis -- базисы отправителя (0 - вычислительный, 1 - Адамара)
+        receiver_basis -- базисы получателя (0 или 1)
+        sender_state -- биты состояния отправителя (0 или 1)
+        name -- название схемы (ключ в словаре circuits)
+
+        Алгоритм:
+        1) Отправитель подготавливает кубиты:
+           - Применяет X к кубитам с состоянием 1 (переключение |0> -> |1>)
+           - Применяет H к кубитам, если базис Адамара (1)
+        2) Получатель измеряет кубиты:
+           - Применяет H перед измерением, если базис Адамара (1)
+        3) Все кубиты измеряются
+        """
+        self.qubits = cirq.LineQubit.range(num_qubits)  # Создание линии кубитов
+        self.circuits[name] = cirq.Circuit()  # Инициализация пустой схемы
+
+        sender_flips = []  # Операции X для изменения состояния с |0> на |1>
+        sender_rotations = []  # Операции H для переключения базиса
+
+        # Цикл по кубитам на стороне отправителя
+        for index in range(num_qubits):
+            if sender_state[index] == 1:
+                sender_flips.append(
+                    cirq.X(self.qubits[index])
+                )  # Применяем X, если бит 1
+            if sender_basis[index] == 1:
+                sender_rotations.append(
+                    cirq.H(self.qubits[index])
+                )  # Применяем H, если выбран базис Адамара
+
+        # Добавляем операции отправителя в схему как отдельные "слои" времени (Moments),
+        # что соответствует выполнению этих операций одновременно независимо для разных кубитов
+        self.circuits[name].append(cirq.Moment(sender_flips))
+        self.circuits[name].append(cirq.Moment(sender_rotations))
+
+        # Операции получателя - подготовка базиса измерения
+        recv_basis_choice = []
+        for index, rbasis in enumerate(receiver_basis):
+            if rbasis == 1:
+                recv_basis_choice.append(
+                    cirq.H(self.qubits[index])
+                )  # Применяем H для базиса Адамара перед измерением
+
+        self.circuits[name].append(
+            cirq.Moment(recv_basis_choice)
+        )  # Добавляем перед измерением
+        self.circuits[name].append(
+            cirq.Moment(cirq.measure_each(*self.qubits))
+        )  # Измеряем все кубиты
+
+    def simulate_base_protocol(self, num_qubits=8, print_legend=None):
+        """
+        Симуляция базового протокола BB84 без присутствия перехватчика (Евы).
+
+        - Генерация параметров
+        - Создание схемы для Алисы и Боба
+        - Симуляция измерений
+        - Выделение ключа на основе совпадения базисов
+        - Печать схемы и результатов
+        """
+        self.setup_demo(num_qubits)
+
+        # Строим схему взаимодействия Алисы и Боба
+        self.make_bb84_circuit(
+            num_qubits, self.alice_basis, self.bob_basis, self.alice_state, name="base"
+        )
+
+        # Запуск симуляции схемы c заданным числом повторений
+        result = cirq.Simulator().run(
+            program=self.circuits["base"], repetitions=self.repetitions
+        )
+
+        # Получение бит результата измерений: для каждого кубита берётся первый результат первого повторения (.item())
+        result_bitstring = bitstring(
+            [int(result.measurements[str(q)][0].item()) for q in self.qubits]
+        )
+
+        # Формируем ключ Боба только для тех позиций, где базисы совпали с Алисой
+        obtained_key = "".join(
+            [
+                result_bitstring[i]
+                for i in range(num_qubits)
+                if self.alice_basis[i] == self.bob_basis[i]
+            ]
+        )
+
+        # Выводим схему и результаты
+        print("########## Печать схемы ##########")
+        print(self.circuits["base"])
+        print("########## Печать результатов ##########")
+        self.print_results(
+            self.alice_basis,
+            self.bob_basis,
+            self.alice_state,
+            self.expected_key,
+            obtained_key,
+        )
+        print(
+            "Общий секретный ключ успешно распределен и может использоваться для шифрования."
+        )
+
+    def simulate_eavesdropped_protocol(self, num_qubits=8, print_legend=None):
+        """
+        Симуляция протокола BB84 с перехватом.
+
+        - Генерация параметров
+        - Симуляция перехвата Евы между Алисой и Бобом
+        - Симуляция передачи от Евы к Бобу
+        - Формирование ключей по разным базисам
+        - Печать схемы и результатов
+        """
+        self.setup_demo(num_qubits)
+
+        # Создаем схему передачи Алиса -> Ева (с их базисами)
+        self.make_bb84_circuit(
+            num_qubits,
+            self.alice_basis,
+            self.eve_basis,
+            self.alice_state,
+            name="alice_eve",
+        )
+
+        # Симуляция измерений Евы
+        result = cirq.Simulator().run(
+            program=self.circuits["alice_eve"], repetitions=self.repetitions
+        )
+        # Ева сохраняет измеренные биты
+        eve_state = [int(result.measurements[str(q)][0].item()) for q in self.qubits]
+
+        # Ева формирует новую схему для передачи Бобу, используя свои измерения и выбранный базис
+        self.make_bb84_circuit(
+            num_qubits, self.eve_basis, self.bob_basis, eve_state, name="eve_bob"
+        )
+
+        # Симуляция измерений Боба
+        result = cirq.Simulator().run(
+            program=self.circuits["eve_bob"], repetitions=self.repetitions
+        )
+        result_bitstring = bitstring(
+            [int(result.measurements[str(q)][0].item()) for q in self.qubits]
+        )
+
+        # Формируем ключи для анализа совпадения (Алиса-Ева и Алиса-Боб)
+        intercepted_key = "".join(
+            [
+                result_bitstring[i]
+                for i in range(num_qubits)
+                if self.alice_basis[i] == self.eve_basis[i]
+            ]
+        )
+        obtained_key = "".join(
+            [
+                result_bitstring[i]
+                for i in range(num_qubits)
+                if self.alice_basis[i] == self.bob_basis[i]
+            ]
+        )
+
+        # Объединяем цепочки схем Алиса->Ева и Ева->Боб для вывода
+        circuit = self.circuits["alice_eve"] + self.circuits["eve_bob"]
+        print("########## Печать схемы ##########")
+        print(circuit)
+        print("########## Печать результатов ##########")
+        self.print_eavesdropped_results(
+            self.alice_basis,
+            self.bob_basis,
+            self.eve_basis,
+            self.alice_state,
+            eve_state,
+            self.expected_key,
+            obtained_key,
+        )
+
+        print("Ева попыталась подслушать и повторно отправить кубиты.")
+        print("Это приводит к ошибкам в ключе, которые обнаруживает Боб.")
+
+    def print_results(
+        self, alice_basis, bob_basis, alice_state, expected_key, obtained_key
+    ):
+        """
+        Вывод результатов базового протокола без перехвата.
+
+        Печатаются:
+        - исходное сообщение Алисы,
+        - базисы Алисы и Боба,
+        - совпадения базисов (в виде X),
+        - ключи Алисы и Боба.
+        """
+        num_qubits = len(alice_basis)
+        basis_match = "".join(
+            ["X" if alice_basis[i] == bob_basis[i] else "_" for i in range(num_qubits)]
+        )
+        alice_basis_str = "".join(
+            ["C" if alice_basis[i] == 0 else "H" for i in range(num_qubits)]
+        )
+        bob_basis_str = "".join(
+            ["C" if bob_basis[i] == 0 else "H" for i in range(num_qubits)]
+        )
+
+        print(f"Только Алиса знает сообщение:\t{bitstring(alice_state)}")
+        print("Базисы Алисы и Боба (C - вычислительный, H - Адамара):")
+        print(f"Алиса: {alice_basis_str}")
+        print(f"Боб:   {bob_basis_str}")
+        print(f"Совпадение базисов:   {basis_match}")
+        print(f"Ключ Алисы: {expected_key}")
+        print(f"Ключ Боба: {obtained_key}")
+
+    def print_eavesdropped_results(
+        self,
+        alice_basis,
+        bob_basis,
+        eve_basis,
+        alice_state,
+        eve_state,
+        expected_key,
+        obtained_key,
+    ):
+        """
+        Вывод результатов протокола с перехватом.
+
+        Печатаются:
+        - исходное сообщение Алисы,
+        - предположение Евы (ее измеренные биты),
+        - базисы Алисы, Евы и Боба,
+        - совпадения базисов,
+        - ключи Алисы, Евы (ставится равным ключу Алисы, тк она пытается угадать),
+          и ключ Боба.
+        """
+        num_qubits = len(alice_basis)
+        eve_basis_match = "".join(
+            ["X" if alice_basis[i] == eve_basis[i] else "_" for i in range(num_qubits)]
+        )
+        bob_basis_match = "".join(
+            ["X" if alice_basis[i] == bob_basis[i] else "_" for i in range(num_qubits)]
+        )
+        alice_basis_str = "".join(
+            ["C" if alice_basis[i] == 0 else "H" for i in range(num_qubits)]
+        )
+        eve_basis_str = "".join(
+            ["C" if eve_basis[i] == 0 else "H" for i in range(num_qubits)]
+        )
+        bob_basis_str = "".join(
+            ["C" if bob_basis[i] == 0 else "H" for i in range(num_qubits)]
+        )
+
+        print(f"Только Алиса знает сообщение:\t{bitstring(alice_state)}")
+        print(f"Предположение Евы о cообщении:\t{bitstring(eve_state)}")
+        print("Базисы (C - вычислительный, H - Адамара):")
+        print(f"Алиса: {alice_basis_str}")
+        print(f"Ева:   {eve_basis_str}")
+        print(f"Боб:   {bob_basis_str}")
+        print(f"Совпадение базисов Алисы и Евы: {eve_basis_match}")
+        print(f"Совпадение базисов Алисы и Боба: {bob_basis_match}")
+        print(f"Ключ Алисы: {expected_key}")
+        print(f"Ключ Евы: {expected_key}")  # По идее, Ева пытается угадать ключ Алисы
+        print(f"Ключ Боба: {obtained_key}")
+
+
+if __name__ == "__main__":
+    bb84 = BB84()
+
+    print("=== БАЗОВЫЙ ПРОТОКОЛ BB84 ===")
+    bb84.simulate_base_protocol(num_qubits=8, print_legend=True)
+
+    print("\n\n=== ПРОТОКОЛ BB84 С ПЕРЕХВАТОМ ===")
+    bb84.simulate_eavesdropped_protocol(num_qubits=8, print_legend=False)

task7/main.py

@@ -0,0 +1,81 @@
+from __future__ import annotations
+
+import random
+
+import cirq
+import numpy as np
+
+
+def make_quantum_teleportation_circuit(ranX, ranY):
+    circuit = cirq.Circuit()
+    msg, alice, bob = cirq.LineQubit.range(3)
+
+    # Создаём запутанную пару кубитов между Алисой и Бобом
+    circuit.append([cirq.H(alice), cirq.CNOT(alice, bob)])
+
+    # Создаём случайное состояние для сообщения
+    circuit.append([cirq.X(msg) ** ranX, cirq.Y(msg) ** ranY])
+
+    # Запутываем сообщение и кубит Алисы
+    circuit.append([cirq.CNOT(msg, alice), cirq.H(msg)])
+
+    # Измеряем сообщение и кубит Алисы
+    circuit.append(cirq.measure(msg, key="msg"))
+    circuit.append(cirq.measure(alice, key="alice"))
+
+    # Используем два классических бита
+    # для восстановления исходного квантового сообщения на кубите Боба
+    circuit.append(cirq.X(bob).with_classical_controls("alice"))
+    circuit.append(cirq.Z(bob).with_classical_controls("msg"))
+
+    return circuit
+
+
+def main(seed=None):
+    random.seed(seed)
+
+    ranX = random.random()
+    ranY = random.random()
+    circuit = make_quantum_teleportation_circuit(ranX, ranY)
+
+    print("Circuit:")
+    print(circuit)
+
+    sim = cirq.Simulator(seed=seed)
+
+    # Запускаем простую симуляцию, которая применяет случайные X и Y-гейты,
+    # создающие наше сообщение
+    q0 = cirq.LineQubit(0)
+    message = sim.simulate(cirq.Circuit([cirq.X(q0) ** ranX, cirq.Y(q0) ** ranY]))
+
+    print("\nBloch Sphere of Message After Random X and Y Gates:")
+    expected = cirq.bloch_vector_from_state_vector(message.final_state_vector, 0)
+    print(
+        "x: ",
+        np.around(expected[0], 4),
+        "y: ",
+        np.around(expected[1], 4),
+        "z: ",
+        np.around(expected[2], 4),
+    )
+
+    final_results = sim.simulate(circuit)
+
+    print("\nBloch Sphere of Qubit 2 at Final State:")
+    teleported = cirq.bloch_vector_from_state_vector(
+        final_results.final_state_vector, 2
+    )
+    print(
+        "x: ",
+        np.around(teleported[0], 4),
+        "y: ",
+        np.around(teleported[1], 4),
+        "z: ",
+        np.around(teleported[2], 4),
+    )
+
+    return expected, teleported
+
+
+if __name__ == "__main__":
+    main()

task8/naive.py

@@ -0,0 +1,89 @@
+from math import acos, cos, pi, sqrt
+
+import cirq
+import numpy as np
+
+
+def naive_kitaev_single(phi_real: float, repetitions: int = 2000):
+    """
+    Наивный вариант алгоритма Китаева
+    ---------------------------------------------------
+    Реализует схему:
+         H --●-- H -- измерение
+               |
+               U
+    где U|1> = e^{i2πφ}|1>.
+
+    После схемы:
+        P(0) = cos²(πφ)
+        P(1) = sin²(πφ)
+
+    Измеряя P(0), можно получить множество кандидатов для φ из уравнения:
+        φ = (±arccos(√P(0)) + mπ) / π
+        где m ∈ ℤ, φ ∈ [0,1).
+    """
+
+    # Симулятор и кубиты
+    sim = cirq.Simulator()
+    a, u = cirq.LineQubit.range(2)  # a — фазовый кубит, u — системный
+
+    # Создаём унитар U с фазой φ
+    U = cirq.MatrixGate(np.diag([1.0, np.exp(1j * 2 * np.pi * phi_real)]))
+
+    # Собираем схему
+    circuit = cirq.Circuit(
+        cirq.X(u),  # подготавливаем |ψ⟩ = |1⟩
+        cirq.H(a),  # суперпозиция на фазовом кубите
+        cirq.ControlledGate(U).on(a, u),
+        cirq.H(a),  # интерференция
+        cirq.measure(a, key="m"),  # измеряем фазовый кубит
+    )
+
+    # Запускаем симуляцию
+    result = sim.run(circuit, repetitions=repetitions)
+    m = result.measurements["m"][:, 0]
+
+    # Эмпирическая вероятность
+    p0_emp = np.mean(m == 0)
+    p0_theory = cos(pi * phi_real) ** 2
+
+    # Решаем уравнение cos²(πφ) = p0_emp
+    c = sqrt(max(0.0, min(1.0, p0_emp)))  # защита от ошибок округления
+    alpha = acos(c)  # угол α ∈ [0, π/2]
+
+    # Все кандидаты φ = (±α + mπ) / π в [0,1)
+    candidates = []
+    for m in range(2):  # m=0 и m=1 дают все уникальные решения в [0,1)
+        for sign in (+1, -1):
+            phi_c = (sign * alpha + m * pi) / pi
+            phi_c_mod = phi_c % 1.0
+            candidates.append((m, sign, phi_c_mod))
+
+    # Удалим дубли и отсортируем
+    uniq = {}
+    for m, sign, val in candidates:
+        key = round(val, 12)
+        if key not in uniq:
+            uniq[key] = (m, sign, val)
+    candidates_unique = sorted(uniq.values(), key=lambda t: t[2])
+
+    print(f"φ(real)      = {phi_real:.6f}")
+    print(f"P(0)_emp     = {p0_emp:.6f}")
+    print(f"P(0)_theory  = {p0_theory:.6f}")
+    print(f"c = sqrt(P0) = {c:.6f},  α = arccos(c) = {alpha:.6f} рад")
+    print()
+    print("Все кандидаты φ (m, sign, φ ∈ [0,1)):")
+    for m, sign, val in candidates_unique:
+        s = "+" if sign > 0 else "-"
+        print(f"  m={m:1d}, sign={s}  →  φ = {val:.12f}")
+    print()
+    print(f"repetitions = {repetitions}")
+
+    return candidates_unique
+
+
+# === Пример запуска ===
+if __name__ == "__main__":
+    phi_real = 0.228  # истинная фаза
+    repetitions = 2000
+    naive_kitaev_single(phi_real, repetitions)

task8/qft.py

@@ -0,0 +1,69 @@
+"""
+Реализация Квантового Преобразования Фурье (QFT) для произвольного числа кубитов n.
+"""
+
+import cirq
+import numpy as np
+
+
+# === Функция построения схемы QFT ===
+def qft_circuit(n: int) -> cirq.Circuit:
+    """Создаёт схему квантового преобразования Фурье (QFT) для n кубитов."""
+    qubits = [cirq.LineQubit(i) for i in range(n)]
+    circuit = cirq.Circuit()
+
+    # Для каждого кубита применяем Hadamard и контролируемые фазовые повороты
+    for i in range(n):
+        circuit.append(cirq.H(qubits[i]))
+        for j in range(i + 1, n):
+            angle = 1 / (2 ** (j - i))
+            circuit.append(cirq.CZ(qubits[j], qubits[i]) ** angle)
+
+    # После всех поворотов меняем порядок кубитов на обратный
+    for i in range(n // 2):
+        circuit.append(cirq.SWAP(qubits[i], qubits[n - i - 1]))
+
+    return circuit
+
+
+# === Главная функция ===
+def main():
+    n = 4  # число кубитов
+    input_state_str = "0111"
+
+    print(f"=== Квантовое преобразование Фурье (QFT) для {n} кубитов ===")
+    print(f"Исходное состояние: |{input_state_str}>")
+
+    # --- Создаём кубиты ---
+    qubits = [cirq.LineQubit(i) for i in range(n)]
+    circuit = cirq.Circuit()
+
+    # Переводим строку в установку нужных кубитов в |1>
+    for i, bit in enumerate(
+        reversed(input_state_str)
+    ):  # reversed, чтобы младший бит был справа
+        if bit == "1":
+            circuit.append(cirq.X(qubits[i]))  # применяем оператор X (NOT)
+
+    # --- Добавляем саму QFT ---
+    circuit += qft_circuit(n)
+
+    print("\nСхема:")
+    print(circuit)
+
+    # --- Симуляция ---
+    simulator = cirq.Simulator()
+    result = simulator.simulate(circuit)
+    state = result.final_state_vector
+
+    # --- Вычисляем  фазы ---
+    phases = np.angle(state)
+
+    print("\n=== Результирующее состояние (амплитуды и фазы) ===")
+    for i, (amp, phi) in enumerate(zip(state, phases)):
+        if abs(amp) > 1e-9:  # пропускаем элементы с амплитудой близкой к 0
+            print(f"|{i:0{n}b}>: amp={abs(amp):.3f}, phase={phi:.3f} rad")
+
+
+if __name__ == "__main__":
+    main()

task8/qpe.py

@@ -0,0 +1,106 @@
+"""
+Адаптированный пример алгоритма квантовой оценки фазы (QPE)
+
+Для каждого значения целевой фазы φ (от 0.0 до 0.9) проверяется,
+насколько точно QPE оценивает её с использованием заданного числа кубитов.
+
+Схема QPE реализует следующую структуру:
+
+   ---H---@------------------|QFT⁻¹|---M---
+          |
+   ---H---@------------------|QFT⁻¹|---M---
+          |
+   ---H---@------------------|QFT⁻¹|---M---
+          |
+   ---|U|--|U²|--|U⁴|--|U⁸|------------
+
+Измеренные биты M дают приближение фазы φ.
+"""
+
+import cirq
+import numpy as np
+
+
+def run_phase_estimation(u_gate, n_qubits: int, repetitions: int):
+    """
+    Строит и симулирует схему QPE (Quantum Phase Estimation).
+
+    Параметры:
+      u_gate:   унитарный оператор U, чья фаза e^{2πiφ} оценивается.
+      n_qubits:       количество кубитов в регистре фазы (точность ~1/2^n).
+      repetitions:    количество повторов измерений (для статистики).
+    """
+
+    # Кубит, к которому применяется оператор U.
+    target = cirq.LineQubit(-1)
+
+    # Регистр для измерения фазы — n_qubits управляющих кубитов.
+    phase_qubits = cirq.LineQubit.range(n_qubits)
+
+    # Контролируемые операторы U^(2^i)
+    controlled_powers = [
+        (u_gate ** (2**i)).on(target).controlled_by(phase_qubits[i])
+        for i in range(n_qubits)
+    ]
+
+    # Схема QPE:
+    circuit = cirq.Circuit(
+        cirq.H.on_each(*phase_qubits),  # 1. Адамар на все управляющие кубиты
+        controlled_powers,  # 2. Контролируемые степени U
+        cirq.qft(*phase_qubits, without_reverse=True)
+        ** -1,  # 3. Обратное квантовое преобразование Фурье
+        cirq.measure(*phase_qubits, key="phase"),  # 4. Измерение регистра фазы
+    )
+
+    simulator = cirq.Simulator()
+    result = simulator.run(circuit, repetitions=repetitions)
+    return result
+
+
+def example_gate(phi: float):
+    """
+    Пример 1-кубитного унитарного оператора U:
+      U|0⟩ = e^{2πiφ}|0⟩,
+      U|1⟩ = |1⟩.
+    """
+    matrix = np.array([[np.exp(2j * np.pi * phi), 0], [0, 1]])
+    return cirq.MatrixGate(matrix)
+
+
+def experiment(n_qubits: int, repetitions=100):
+    """
+    Запускает серию экспериментов для разных значений фазы φ.
+    """
+    print(f"\nТест с {n_qubits} кубитами")
+    errors = []
+
+    for target_phi in np.arange(0, 1.0, 0.1):
+        gate = example_gate(target_phi)
+        result = run_phase_estimation(gate, n_qubits, repetitions)
+
+        # Наиболее часто встречающееся значение результата измерения
+        mode = result.data["phase"].mode()[0]
+
+        # Переводим измеренный результат в оценку фазы
+        phi_est = mode / (2**n_qubits)
+
+        print(
+            f"target={target_phi:0.4f}, estimate={phi_est:0.4f} = {mode}/{2**n_qubits}"
+        )
+        errors.append((target_phi - phi_est) ** 2)
+
+    rms = np.sqrt(np.mean(errors))
+    print(f"Среднеквадратичная ошибка (RMS): {rms:.4f}")
+
+
+def main():
+    """
+    Запускает эксперименты для разных количеств кубитов:
+    2, 4 и 8 — для демонстрации роста точности.
+    """
+    for n in (2, 4, 8):
+        experiment(n, repetitions=200)
+
+
+if __name__ == "__main__":
+    main()

task9/shor_cq.py

@@ -0,0 +1,357 @@
+# pylint: disable=wrong-or-nonexistent-copyright-notice
+"""Demonstrates Shor's algorithm.
+
+Shor's algorithm [1] is a quantum algorithm for integer factorization. Given
+a composite integer n, it finds its non-trivial factor d, i.e. a factor other
+than 1 or n.
+
+The algorithm consists of two parts: quantum order-finding subroutine and
+classical probabilistic reduction of the factoring problem to the order-
+finding problem. Given two positive integers x and n, the order-finding
+problem asks for the smallest positive integer r such that x**r mod n == 1.
+
+The classical reduction algorithm first handles two corner cases which do
+not rely on quantum computation: when n is even or a prime power. For other
+n, the algorithm first draws a random x uniformly from 2..n-1 and then uses
+the quantum order-finding subroutine to compute the order r of x modulo n,
+i.e. it finds the smallest positive integer r such that x**r == 1 mod n. Now,
+if r is even, then y = x**(r/2) is a solution to the equation
+
+    y**2 == 1 mod n.                                                   (*)
+
+It's easy to see that in this case gcd(y - 1, n) or gcd(y + 1, n) divides n.
+If in addition y is a non-trivial solution, i.e. if it is not equal to -1,
+then gcd(y - 1, n) or gcd(y + 1, n) is a non-trivial factor of n (note that
+y cannot be 1). If r is odd or if y is a trivial solution of (*), then the
+algorithm is repeated for a different random x.
+
+It turns out [1] that the probability of r being even and y = x**(r/2) being
+a non-trivial solution of equation (*) is at least 1 - 1/2**(k - 1) where k
+is the number of distinct prime factors of n. Since the case k = 1 has been
+handled by the classical part, we have k >= 2 and the success probability of
+a single attempt is at least 1/2.
+
+The subroutine for finding the order r of a number x modulo n consists of two
+steps. In the first step, Quantum Phase Estimation is applied to a unitary such
+as
+
+    U|y⟩ = |xy mod n⟩  0 <= y < n
+    U|y⟩ = |y⟩         n =< y
+
+whose eigenvalues are s/r for s = 0, 1, ..., r - 1. In the second step, the
+classical continued fractions algorithm is used to recover r from s/r. Note
+that when gcd(s, r) > 1 then an incorrect r is found. This can be detected
+by verifying that r is indeed the order of x. If it is not, Quantum Phase
+Estimation algorithm is retried.
+
+[1]: https://arxiv.org/abs/quant-ph/9508027
+"""
+
+from __future__ import annotations
+
+import argparse
+import fractions
+import math
+import random
+from typing import Callable, Sequence
+
+import cirq
+import sympy
+
+parser = argparse.ArgumentParser(description="Factorization demo.")
+parser.add_argument("n", type=int, help="composite integer to factor")
+parser.add_argument(
+    "--order_finder",
+    type=str,
+    choices=("naive", "quantum"),
+    default="naive",
+    help=(
+        'order finder to use; must be either "naive" '
+        'for a naive classical algorithm or "quantum" '
+        "for a quantum circuit; note that in practice "
+        '"quantum" is substantially slower since it '
+        "incurs the overhead of classical simulation."
+    ),
+)
+
+
+def naive_order_finder(x: int, n: int) -> int | None:
+    """Computes smallest positive r such that x**r mod n == 1.
+
+    Args:
+        x: integer whose order is to be computed, must be greater than one
+           and belong to the multiplicative group of integers modulo n (which
+           consists of positive integers relatively prime to n),
+        n: modulus of the multiplicative group.
+
+    Returns:
+        Smallest positive integer r such that x**r == 1 mod n.
+        Always succeeds (and hence never returns None).
+
+    Raises:
+        ValueError: When x is 1 or not an element of the multiplicative
+            group of integers modulo n.
+    """
+    if x < 2 or n <= x or math.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(cirq.ArithmeticGate):
+    """Quantum modular exponentiation.
+
+    This class represents the unitary which multiplies base raised to exponent
+    into the target modulo the given modulus. More precisely, it represents the
+    unitary V which computes modular exponentiation x**e mod n:
+
+        V|y⟩|e⟩ = |y * x**e mod n⟩ |e⟩     0 <= y < n
+        V|y⟩|e⟩ = |y⟩ |e⟩                  n <= y
+
+    where y is the target register, e is the exponent register, x is the base
+    and n is the modulus. Consequently,
+
+        V|y⟩|e⟩ = (U**e|r⟩)|e⟩
+
+    where U is the unitary defined as
+
+        U|y⟩ = |y * x mod n⟩      0 <= y < n
+        U|y⟩ = |y⟩                n <= y
+
+    in the header of this file.
+
+    Quantum order finding algorithm (which is the quantum part of the Shor's
+    algorithm) uses quantum modular exponentiation together with the Quantum
+    Phase Estimation to compute the order of x modulo n.
+    """
+
+    def __init__(
+        self,
+        target: Sequence[int],
+        exponent: int | Sequence[int],
+        base: int,
+        modulus: int,
+    ) -> None:
+        if len(target) < modulus.bit_length():
+            raise ValueError(
+                f"Register with {len(target)} qubits is too small for modulus {modulus}"
+            )
+        self.target = target
+        self.exponent = exponent
+        self.base = base
+        self.modulus = modulus
+
+    def registers(self) -> Sequence[int | Sequence[int]]:
+        return self.target, self.exponent, self.base, self.modulus
+
+    def with_registers(self, *new_registers: 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"Target must be a qubit register, got {type(target)}")
+        if not isinstance(base, int):
+            raise ValueError(f"Base must be a classical constant, got {type(base)}")
+        if not isinstance(modulus, int):
+            raise ValueError(
+                f"Modulus must be a classical constant, got {type(modulus)}"
+            )
+        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: cirq.CircuitDiagramInfoArgs
+    ) -> cirq.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 cirq.CircuitDiagramInfo(wire_symbols=tuple(wire_symbols))
+
+
+def make_order_finding_circuit(x: int, n: int) -> cirq.Circuit:
+    """Returns quantum circuit which computes the order of x modulo n.
+
+    The circuit uses Quantum Phase Estimation to compute an eigenvalue of
+    the unitary
+
+        U|y⟩ = |y * x mod n⟩      0 <= y < n
+        U|y⟩ = |y⟩                n <= y
+
+    discussed in the header of this file. The circuit uses two registers:
+    the target register which is acted on by U and the exponent register
+    from which an eigenvalue is read out after measurement at the end. The
+    circuit consists of three steps:
+
+    1. Initialization of the target register to |0..01⟩ and the exponent
+       register to a superposition state.
+    2. Multiple controlled-U**2**j operations implemented efficiently using
+       modular exponentiation.
+    3. Inverse Quantum Fourier Transform to kick an eigenvalue to the
+       exponent register.
+
+    Args:
+        x: positive integer whose order modulo n is to be found
+        n: modulus relative to which the order of x is to be found
+
+    Returns:
+        Quantum circuit for finding the order of x modulo n
+    """
+    L = n.bit_length()
+    target = cirq.LineQubit.range(L)
+    exponent = cirq.LineQubit.range(L, 3 * L + 3)
+    return cirq.Circuit(
+        cirq.X(target[L - 1]),
+        cirq.H.on_each(*exponent),
+        ModularExp([2] * len(target), [2] * len(exponent), x, n).on(*target + exponent),
+        cirq.qft(*exponent, inverse=True),
+        cirq.measure(*exponent, key="exponent"),
+    )
+
+
+def read_eigenphase(result: cirq.Result) -> float:
+    """Interprets the output of the order finding circuit.
+
+    Specifically, it returns s/r such that exp(2πis/r) is an eigenvalue
+    of the unitary
+
+        U|y⟩ = |xy mod n⟩  0 <= y < n
+        U|y⟩ = |y⟩         n <= y
+
+    described in the header of this file.
+
+    Args:
+        result: trial result obtained by sampling the output of the
+            circuit built by make_order_finding_circuit
+
+    Returns:
+        s/r where r is the order of x modulo n and s is in [0..r-1].
+    """
+    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) -> int | None:
+    """Computes smallest positive r such that x**r mod n == 1.
+
+    Args:
+        x: integer whose order is to be computed, must be greater than one
+           and belong to the multiplicative group of integers modulo n (which
+           consists of positive integers relatively prime to n),
+        n: modulus of the multiplicative group.
+
+    Returns:
+        Smallest positive integer r such that x**r == 1 mod n or None if the
+        algorithm failed. The algorithm fails when the result of the Quantum
+        Phase Estimation is inaccurate, zero or a reducible fraction.
+
+    Raises:
+        ValueError: When x is 1 or not an element of the multiplicative
+            group of integers modulo n.
+    """
+    if x < 2 or n <= x or math.gcd(x, n) > 1:
+        raise ValueError(f"Invalid x={x} for modulus n={n}.")
+
+    circuit = make_order_finding_circuit(x, n)
+    result = cirq.sample(circuit)
+    eigenphase = read_eigenphase(result)
+    f = fractions.Fraction.from_float(eigenphase).limit_denominator(n)
+    if f.numerator == 0:
+        return None  # pragma: no cover
+    r = f.denominator
+    if x**r % n != 1:
+        return None  # pragma: no cover
+    return r
+
+
+def find_factor_of_prime_power(n: int) -> int | None:
+    """Returns non-trivial factor of n if n is a prime power, else None."""
+    for k in range(2, math.floor(math.log2(n)) + 1):
+        c = math.pow(n, 1 / k)
+        c1 = math.floor(c)
+        if c1**k == n:
+            return c1
+        c2 = math.ceil(c)
+        if c2**k == n:
+            return c2
+    return None
+
+
+def find_factor(
+    n: int, order_finder: Callable[[int, int], int | None], max_attempts: int = 30
+) -> int | None:
+    """Returns a non-trivial factor of composite integer n.
+
+    Args:
+        n: integer to factorize,
+        order_finder: function for finding the order of elements of the
+            multiplicative group of integers modulo n,
+        max_attempts: number of random x's to try, also an upper limit
+            on the number of order_finder invocations.
+
+    Returns:
+        Non-trivial factor of n or None if no such factor was found.
+        Factor k of n is trivial if it is 1 or n.
+    """
+    if sympy.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 = random.randint(2, n - 1)
+        c = math.gcd(x, n)
+        if 1 < c < n:
+            return c  # pragma: no cover
+        r = order_finder(x, n)
+        if r is None:
+            continue  # pragma: no cover
+        if r % 2 != 0:
+            continue  # pragma: no cover
+        y = x ** (r // 2) % n
+        assert 1 < y < n
+        c = math.gcd(y - 1, n)
+        if 1 < c < n:
+            return c
+    return None  # pragma: no cover
+
+
+def main(n: int, order_finder: Callable[[int, int], int | None] = naive_order_finder):
+    if n < 2:
+        raise ValueError(
+            f"Invalid input {n}, expected positive integer greater than one."
+        )
+
+    d = find_factor(n, order_finder)
+
+    if d is None:
+        print(f"No non-trivial factor of {n} found. It is probably a prime.")
+    else:
+        print(f"{d} is a non-trivial factor of {n}")
+
+        assert 1 < d < n
+        assert n % d == 0
+
+
+if __name__ == "__main__":  # pragma: no cover
+    ORDER_FINDERS = {"naive": naive_order_finder, "quantum": quantum_order_finder}
+    args = parser.parse_args()
+    main(n=args.n, order_finder=ORDER_FINDERS[args.order_finder])