From 7c58e681b231da8fce81560bf82763ecb988d621 Mon Sep 17 00:00:00 2001 From: Arity-T Date: Wed, 19 Nov 2025 11:27:58 +0300 Subject: [PATCH] =?UTF-8?q?=D0=B1=D0=BE=D0=BB=D1=8C=D1=88=D0=B5=20=D0=BA?= =?UTF-8?q?=D0=BE=D0=B4=D0=BE=D0=B2=D0=B2=D0=B2=D0=B2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- task10/grover.py | 118 ++++++++++++++++ task8/naive.py | 89 ++++++++++++ task8/qft.py | 69 +++++++++ task8/qpe.py | 106 ++++++++++++++ task9/shor_cq.py | 357 +++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 739 insertions(+) create mode 100644 task10/grover.py create mode 100644 task8/naive.py create mode 100644 task8/qft.py create mode 100644 task8/qpe.py create mode 100644 task9/shor_cq.py diff --git a/task10/grover.py b/task10/grover.py new file mode 100644 index 0000000..53dbcf3 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/task8/naive.py b/task8/naive.py new file mode 100644 index 0000000..b7737c5 --- /dev/null +++ b/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) diff --git a/task8/qft.py b/task8/qft.py new file mode 100644 index 0000000..dc39822 --- /dev/null +++ b/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() diff --git a/task8/qpe.py b/task8/qpe.py new file mode 100644 index 0000000..23bdea5 --- /dev/null +++ b/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() diff --git a/task9/shor_cq.py b/task9/shor_cq.py new file mode 100644 index 0000000..65af721 --- /dev/null +++ b/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])