qqz.order_finding

 1from typing import Optional
 2from math import gcd
 3
 4import numpy as np
 5import matplotlib.pyplot as plt
 6from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, Aer, transpile, assemble
 7from qiskit.visualization import plot_histogram
 8from sympy import Rational
 9from sympy.ntheory.continued_fraction import continued_fraction, continued_fraction_convergents
10
11from qft import qft
12from elementary import ax_modM
13from classical_utils import lcm
14
15def order_finding(x: int, N: int, show_hist: Optional[bool] = False) -> int:
16#def order_finding(x: int, N: int, epsilon: Optional[float] = 0.2, show_hist: Optional[bool] = False) -> int:
17    r"""Order-finding algorithm: it finds $r$ of $x^r\equiv 1\pmod N$. It requires 
18
19    Args:
20        x (int): $x$
21        N (int): $N$
22
23    Returns:
24        order $r$
25
26    Examples:
27
28    ```
29    >>> order_finding(x=3, N=5, show_hist=True)
30    4
31    ```
32
33    and get below image in `img` directory:  
34    ![](../img/order_finding_x3_N5.png)  
35    It represents $0/2^6=0$, $2^4/2^6=1/4$, $2^5/2^6=1/2$, and $(2^4+2^5)/2^6=3/4$ from the left.
36    This answer is $r=4$, so $1/2$ looks wrong.
37    However, $\tilde{r}=2$ is a factor of $r$, so we can get correct $r$ by lcm with another $\tilde{r}$.
38    """
39
40    L = int(np.ceil(np.log2(N)))
41    t = 2 * L# + 1 + int(np.ceil(np.log2(3 + 1 / (2 * epsilon)))) # epsilon requires too many qubits to run this program...
42
43    first_register = QuantumRegister(t)
44    second_register = QuantumRegister(2 * t)
45    auxiliary_register_mid = QuantumRegister(t)
46    auxiliary_register_end = QuantumRegister(6 * t - 2)
47    classical_register = ClassicalRegister(len(first_register))
48
49    qc = QuantumCircuit(first_register, auxiliary_register_mid, second_register, auxiliary_register_end, classical_register)
50#    qc.add_register(first_register)
51#    qc.add_register(second_register)
52#    qc.add_register(classical_register)
53
54    qc.h(first_register)
55
56    #import pdb; pdb.set_trace()
57    #qc.append(ax_modM(a=x, M=N, N_len=len(first_register)), [first_register, auxiliary_register_mid, second_register, auxiliary_register_end])
58    qc.append(ax_modM(a=x, M=N, N_len=len(first_register)), qc.qubits[:10 * t - 2])
59
60    qc.append(qft(n=len(first_register)).inverse(), first_register)
61
62    qc.measure(first_register, classical_register)
63
64    backend = Aer.get_backend('aer_simulator_matrix_product_state')#('aer_simulator')
65    qc = transpile(qc, backend)
66    job = backend.run(qc, shots=10000)
67    hist = job.result().get_counts()
68
69    if show_hist:
70        plot_histogram(hist)
71        plt.savefig(f'img/order_finding_x{x}_N{N}.png', bbox_inches='tight')
72        #plt.savefig(f'img/order_finding_x{x}_N{N}_eps{epsilon}.png', bbox_inches='tight')
73
74    all_fractions = []
75    for measured_key, _ in sorted(hist.items(), key=lambda x: x[1], reverse=True):
76        tilde_s_per_r = Rational(int(measured_key[-t:], 2), 2 ** t)
77        if tilde_s_per_r == 0:
78            continue
79
80        fractions = []
81        for fraction in continued_fraction_convergents(continued_fraction(tilde_s_per_r)):
82            if pow(x, fraction.denominator, N) == 1:
83                return fraction.denominator
84            fractions.append(fraction)
85
86            for other_fraction in all_fractions:
87                if math.gcd(fraction.numerator, other_fraction.numerator) == 1:
88                    r_candidate = lcm(fraction.denominator, other_fraction.denominator)
89                    if pow(x, r_candidate, N) == 1:
90                        return r_candidate
91
92    raise Exception('r is NOT found!')
93
94
95if __name__ == '__main__':
96    #print(order_finding(x=5, N=21, show_hist=True))
97    print(order_finding(x=3, N=5, show_hist=True))