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

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()

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])