Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
The Discrete Logarithm Problem [1] was shown by Shor [2] to be solved in a polynomial time using quantum computers, while the fastest classical algorithms take a superpolynomial time. The problem is at least as hard as the factoring problem. In fact, the hardness of the problem is the basis for the Diffie-Hellman [3] protocol for key exchange. The algorithm is a specific instance of the Abelian Hidden Subgroup Problem [4]. The algorithm treats the following problem:
  • Input: A cyclic group G=⟨g⟩G = \langle g \rangle with a generator g∈Gg\in G, an element x∈Gx\in G and the order of the group is known r=∣G∣r=|G|. A quantum oracle UfU_f applying the transformation
Uf∣α⟩∣β⟩∣0⟩=∣α⟩∣β⟩∣f(α,β)⟩  ,U_f|\alpha\rangle |\beta\rangle |0\rangle = |\alpha\rangle |\beta\rangle |f(\alpha, \beta)\rangle~~, where f(α,β)=xαgβf(\alpha, \beta) = x^\alpha g^\beta with α,β∈Zr\alpha, \beta \in \mathbb{Z}_r.
  • Promise: There is a positive integer s∈Ns\in \mathbb{N} such that gs=xg^s = x.
  • Output: The least positive integer s=log⁡gxs = \log_g x, i.e, the discrete logarithm.
  • Complexity: A single application of the algorithm succeeds with probability O(log⁡log⁡r)O(\log \log r) and requires O((log⁡t)2)O((\log t)^2) time, where t=⌈log⁡r+log⁡(1/ϵ)⌉t = \lceil \log r +\log(1/\epsilon) \rceil and 1−ϵ1-\epsilon is success probability. Hence, boosting the success probability to a constant O(1)O(1) via repetition yields a total expected runtime of
O((log⁡t)2log⁡log⁡r) .O((\log t)^2 \log \log r)~.
Keywords: Abelain Hidden subgroup problem, quantum Fourier transform, period finding.

​
Algorithm Steps

We first consider the case where r=2mr=2^m and later generalize. We begin by introducing an input state ∣ψ0⟩=∣0m⟩∣0m⟩∣0R⟩  ,|\psi_0\rangle = |0^m\rangle|0^m\rangle|0^R\rangle~~, where RR is the number of bits required to represent ff‘s image. In addition, we note that f(α,β)=xαgβ=gαlog⁡gx+βf(\alpha,\beta) = x^\alpha g^\beta = g^{\alpha\log_g x + \beta} is constant on lines satisfying {(α,β)∈Zr2,αlog⁡gx+β=λmod  r} .\{(\alpha,\beta)\in \mathbb{Z}_r^2, \alpha\log_g x + \beta = \lambda \mod r\}~. The algorithm includes four main steps, and has a similar structure to the other quantum algorithms for the hidden subgroup problem.
  1. We first prepare a uniform superposition over the input states of the first two registers:
∣0m⟩∣0m⟩∣0m⟩→H⊗m1r∑α,β∈Zr∣α⟩∣β⟩∣0m⟩  .|0^m\rangle|0^m\rangle|0^m\rangle \xrightarrow{H^{\otimes m}} \frac{1}{r}\sum_{\alpha,\beta \in \mathbb{Z}_r}|\alpha \rangle |\beta \rangle |0^m\rangle ~~.
  1. Perform a query to the oracle
→Uf1r∑α,β∈Zr∣α⟩∣β⟩∣f(α,β)⟩  .\xrightarrow{U_f} \frac{1}{r}\sum_{\alpha, \beta \in \mathbb{Z}_r}|\alpha \rangle |\beta \rangle |f(\alpha, \beta)\rangle~~. Utilizing the periodicity of ff, we can express the state as =1r∑α,λ∈Zr∣α⟩∣λ−αlog⁡gx⟩∣gλ⟩  .= \frac{1}{r}\sum_{\alpha, \lambda \in \mathbb{Z}_r}|\alpha \rangle |\lambda -\alpha \log_g x \rangle |g^{\lambda}\rangle~~.
  1. Perform a double inverse Fourier transform over the group Zr\mathbb{Z}_r
→QFTZr†×QFTZr†1r2∑λ,μ,ν∈Zr(∑αe−i2πα(μ−νlog⁡gx)/r)e−i2πλν/r∣μ⟩∣ν⟩∣gλ⟩\xrightarrow{\text{QFT}_{\mathbb{Z}_r}^\dagger \times \text{QFT}_{\mathbb{Z}_r}^\dagger } \frac{1}{r^2}\sum_{\lambda, \mu, \nu \in \mathbb{Z}_r}\left(\sum_\alpha e^{-i 2\pi \alpha(\mu - \nu \log_g x)/r}\right)e^{-i 2\pi \lambda \nu/r}|\mu \rangle |\nu \rangle |g^{\lambda}\rangle =1r∑λ,ν∈Zrei2πλν/r∣νlog⁡gx⟩∣ν⟩∣gλ⟩  ,= \frac{1}{r}\sum_{\lambda, \nu \in \mathbb{Z}_r}e^{i 2\pi \lambda \nu/r}|\nu \log_g x \rangle |\nu \rangle |g^{\lambda}\rangle~~, where the Fourier transform over Zr\mathbb{Z}_r is QFTZr∣α⟩=1r∑μ∈Zrei2παμ/r∣μ⟩\text{QFT}_{\mathbb{Z}_r} |\alpha \rangle = \frac{1}{\sqrt{r}}\sum_{\mu\in \mathbb{Z}_r} e^{i 2\pi \alpha \mu/r}|\mu \rangle and in the second line we utilized the identity ∑αei2παη=rδ0η\sum_\alpha e^{i 2\pi \alpha \eta} = r \delta_{0\eta} (easily derived by use of a sum over a geometric series). 4. Measure the three quantum variables. The measurement of the third variable collapses the quantum state to ∣gλ⟩|g^\lambda \rangle with a uniform distribution. After the collapse, the resulting state is independent of λ\lambda, therefore, we can discard this measurement outcome. From the measurement of the second quantum variable, we obtain with a uniform distribution over ν\nu the outcome νlog⁡gx\nu \log_g x. With a probablility of order O(log⁡log⁡r)O(\log \log r), the obtained result ν\nu is co-prime to rr, and there exists a modular rr multiplicative inverse ν−1\nu^{-1}. Under this condition, we can multiply the outcome of the second register to obtain the discrete log s=log⁡gxs=\log_g x. Hence, repeating the experiment O(log⁡log⁡r)O(\log \log r) times leads to a success probability of order O(1)O(1).

​
Building the Algorithm with Classiq

We begin by importing software packages 
import matplotlib.pyplot as plt
import numpy as np

from classiq import *
from classiq.qmod.symbolic import ceiling, log
We consider two exemplary cases. First, we study the group Z5×\mathbb{Z}_5^\times, where the order is a power of r=4=22r=4=2^2. For this case, quantum variables with only log⁡r=m\log r = m qubits are required. Following, the group Z13×\mathbb{Z}_{13}^\times exemplifies the alternative case where the order r≠2m r\neq 2^m for an integer mm. For these examples, we denote the modulus by NN. Since, these are cyclic groups, we have N=r+1N = r+1, therefore R=⌈log⁡N⌉R = \lceil \log N\rceil. In the following examples, the third register, containing the function f(α,β)f(\alpha, \beta ) after the second step, is represented by ⌈log⁡N⌉\lceil \log N\rceil. The heart of the algorithm’s logic is the implementation of the function ∣α⟩∣β⟩∣0⟩→∣α⟩∣β⟩∣xαgβ⟩ .|\alpha \rangle|\beta \rangle|0\rangle \rightarrow |\alpha \rangle|\beta \rangle|x^{\alpha } g^{\beta}\rangle~. This is done using two applications of the modular exponentiation function, described in detail in the Shor’s Factoring Algorithm notebook. So here we import it from the Classiq library. The modular_exponentiation function defined below accepts these arguments:
  • N: CInt - modulo number
  • a: CInt - base of the exponentiation
  • x: QArray[QBit] - unsigned integer to multiply by the exponentiation
  • pw: QArray[QBit]- power of the exponentiation
So the function implements ∣pw⟩∣x⟩→∣pw⟩∣x⋅apwmod  N⟩|pw\rangle|x\rangle \rightarrow |pw\rangle|x \cdot a ^ {pw}\mod N\rangle. 
from classiq import *


@qfunc
def modular_exponentiation(N: CInt, a: CInt, x: QArray, pw: QArray):
    repeat(
        count=pw.len,
        iteration=lambda index: control(
            pw[index],
            # lambda: inplace_modular_multiply(N, (a ** (2**index)) % N, x),  #modular_multiply_constant_inplace
            lambda: modular_multiply_constant_inplace(N, a ** (2**index), x),
        ),
    )


@qfunc
def discrete_log_oracle(
    g_generator: CInt,
    x_element: CInt,
    N_modulus: CInt,
    alpha: QArray,
    beta: QArray,
    func_res: Output[QNum],
) -> None:
    allocate(ceiling(log(N_modulus, 2)), func_res)

    func_res ^= 1
    modular_exponentiation(N_modulus, x_element, func_res, alpha)
    modular_exponentiation(N_modulus, g_generator, func_res, beta)

​

​
Full Algorithm

  1. Prepare a uniform superposition over the first two quantum variables alpha, beta.
Each variable has size ⌈log⁡r⌉+log⁡(1/ϵ)\lceil \log r\rceil + \log({1/{\epsilon}}). In the special case where rr is a power of 2, log⁡r\log r is enough.
  1. Compute discrete_log_oracle on the func_res variable. func_res is of size ⌈log⁡N⌉\lceil \log N\rceil.
  2. Apply the inverse Fourier transform alpha, beta.
  3. Measure.
@qfunc
def discrete_log(
    g: CInt,
    x: CInt,
    N: CInt,
    order: CInt,
    alpha: Output[QArray],
    beta: Output[QArray],
    func_res: Output[QArray],
) -> None:
    reg_len = ceiling(log(order, 2))
    allocate(reg_len, alpha)
    allocate(reg_len, beta)

    hadamard_transform(alpha)
    hadamard_transform(beta)

    discrete_log_oracle(g, x, N, alpha, beta, func_res)

    invert(lambda: qft(alpha))
    invert(lambda: qft(beta))
After the inverse QFTs (under the assumption of r=2mr=2^m for some mm, therefore t=2rt=2^r), the variables become ∣ψ⟩=1r∑λ,ν∈Zrei2πλν/r∣νlog⁡gx⟩α∣ν⟩β∣gλ⟩res|\psi\rangle = \frac{1}{r}\sum_{\lambda, \nu \in \mathbb{Z}_r}e^{i 2\pi \lambda \nu/r}|\nu \log_g x \rangle_{\alpha} |\nu \rangle_{\beta} |g^{\lambda}\rangle_{\text{res}} where we added a subscript to the quantum variables for clarity, and ∣⟩res| \rangle_{\text{res}} designates the func_res variable. If ν∈Zr\nu\in \mathbb{Z}_r has an inverse, ν−1∈Zr\nu^{-1}\in\mathbb{Z}_r we can extract s=log⁡xgs=\log_x g from variable α\alpha by multiplying by ν−1\nu^{-1}. See the second example for the general case, for which r≠2mr \neq 2^m. Screenshot 2025-12-10 at 16.38.37.png

​
Example: G=Z5×G = \mathbb{Z}_5^\times

For this specific demonstration, we choose G=Z5×G = \mathbb{Z}_5^\times, with g=3g=3 and x=2x=2. With this setting, log⁡gx=3\log_gx=3. We choose this specific example because the order of the group r=4r=4 is a power of 22, so we can get the exact discrete logarithm without continued-fractions postprocessing. In other cases, we use a larger quantum variable for the exponents so the continued fractions postprocessing converges. 
MODULU_NUM = 5
G_GENERATOR = 3
X_LOG_ARG = 2
ORDER = MODULU_NUM - 1  # as 5 is prime


@qfunc
def main(
    alpha: Output[QNum],
    beta: Output[QNum],
    func_res: Output[QNum],
) -> None:
    discrete_log(G_GENERATOR, X_LOG_ARG, MODULU_NUM, ORDER, alpha, beta, func_res)


qmod_Z5 = create_model(
    main,
    constraints=Constraints(max_width=13),
    preferences=Preferences(optimization_level=1),
    execution_preferences=ExecutionPreferences(num_shots=4000),
)

qprog_Z5 = synthesize(qmod_Z5)
show(qprog_Z5)
Output:

Quantum program link: https://platform.classiq.io/circuit/36zUExgS7ndoE6IUFudD7cLnNii
  


result_Z5 = execute(qprog_Z5).result_value()
result_Z5.dataframe
alphabetafunc_rescountprobabilitybitstring
02232680.067000111010
12222660.066500101010
20032660.066500110000
30012640.066000010000
43112590.064750010111
52212580.064500011010
61332560.064000111101
71312500.062500011101
81342500.062501001101
93132460.061500110111
101322440.061000101101
110042400.060001000000
123122380.059500100111
130022360.059000100000
142242300.057501001010
153142290.057251000111
Note that func_res is uncorrelated to the other variables, and we get uniform distribution, as expected. We take only the beta that are co-prime to r=4r=4, so they have a multiplicative-inverse. Hence beta=1,3 are the relevant results. So we get two relevant results (for all different λ\lambdas): ∣1⟩∣3⟩|1\rangle|3\rangle, ∣3⟩∣1⟩|3\rangle|1\rangle. All that remains to get the logarithm is to multiply alpha by the inverse of beta: 
for res in result_Z5.parsed_counts:
    if res.state["beta"] in [1, 3]:
        logarithm = res.state["beta"] * pow(res.state["alpha"], -1, 4)
        assert logarithm == 3
Verify we received the correct discrete logarithm: 
log_arg = (G_GENERATOR**logarithm) % MODULU_NUM
print(log_arg)
assert log_arg == X_LOG_ARG
Output:
2
And, indeed, both cases give the same result, which is exactly the discrete logarithm: log⁡32mod  5=3\log_32 \mod 5 = 3.

​
Generalization for the case where r≠2mr\neq 2^m

For an arbitrary r∈Nr\in \mathbb{N}, we employ two t=⌈log⁡r⌉+log⁡(1/ϵ)t=\lceil \log r\rceil + \log (1/\epsilon) qubit quantum variables and a third RR-qubit register, initialized to the state ∣ψ0⟩=∣0t⟩∣0t⟩∣0R⟩|\psi_0\rangle = |0^t\rangle|0^t\rangle|0^R\rangle , where R=⌈log⁡N⌉R = \lceil \log N\rceil. The derivation follows a similar structure, where the initial sums are now α,β∈ZT\alpha, \beta \in \mathbb Z_T, with T=2tT=2^t, while the period of ff remains rr, i.e., λ∈Zr\lambda\in \mathbb{Z}_r. Therefore, after the second step we obtain the state 1T∑α,β∈ZT∣α⟩∣β⟩∣0R⟩→Uf1T∑α,β∈ZT∣α⟩∣β⟩∣gαlog⁡gx+β⟩  .\frac{1}{T}\sum_{\alpha,\beta \in \mathbb{Z}_T}|\alpha\rangle|\beta\rangle|0^R\rangle\xrightarrow{U_f} \frac{1}{T}\sum_{\alpha,\beta\in \mathbb{Z}_T}|\alpha \rangle |\beta\rangle |g^{\alpha \log_g x + \beta}\rangle~~. Next, we express the periodic state in an alternative form. We introduce the operator UgU_g, which satisfies Ug∣h⟩=∣hg⟩U_g |h\rangle= |h g\rangle, where ee is the identity element of the group. Therefore ∣gλ⟩=Ugλ∣e⟩|{g^\lambda}\rangle = U_g^\lambda |e \rangle. The state ∣e⟩|e\rangle can be expressed as a uniform superposition of the eigenstates of UgU_g: ∣e⟩=1r∑k=0r−1∣Ψk⟩  ,|e \rangle = \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1} |\Psi_k\rangle~~, where ∣Ψk⟩=1r∑k=0r−1ei2πkk′/r∣gk′⟩| \Psi_k \rangle = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{i2\pi k k'/r}| g^{k'}\rangle, which satisfy Ug∣Ψk⟩=ei2πk/r∣Ψk⟩U_g | \Psi_k \rangle = e^{i 2\pi k/r}| \Psi_k\rangle. These relations allow writing the state after the oracle operation as 1T∑α,β∈ZT∣α⟩∣β⟩∑k=0r−1ei2π(αlog⁡gx+β)k/r∣Ψk⟩  .\frac{1}{T}\sum_{\alpha,\beta\in \mathbb{Z}_T}|\alpha \rangle |\beta\rangle \sum_{k=0}^{r-1}e^{i 2\pi (\alpha \log_g x + \beta) k/r }|\Psi_k\rangle~~. Applying the inverse quantum Fourier transform leads to the state →QFTT†×QFTT†1T∑k=0r−1∑μ,ν∈ZT(∑α∈ZTei2πα(log⁡gxk/r−μ/T)∣μ⟩)(∑β∈ZTei2πβ(k/r−ν/T)∣ν⟩)∣Ψk⟩  .\xrightarrow{\text{QFT}_T^\dagger \times \text{QFT}_T^\dagger}\frac{1}{T}\sum_{k=0}^{r-1}\sum_{\mu,\nu\in \mathbb{Z}_T}\left(\sum_{\alpha\in \mathbb{Z}_T} e^{i 2\pi \alpha (\log_g x k/r-\mu/T) }| \mu\rangle\right)\left( \sum_{\beta\in \mathbb{Z}_T} e^{i 2\pi \beta (k/r-\nu /T) }| \nu \rangle \right)|\Psi_k\rangle~~. Utilizing the geometric sum, one can show that the functions are peaked (with a width O(1/T)O(1/T)) around μ≈Tlog⁡gxk/r\mu \approx T\log_g x k/r and ν≈Tk/r \nu \approx T k / r, correspondingly. To evaluate the values of klog⁡gxk \log_g x and kk we measure alpha and beta registers, multiply by r/Nr/N and round. For sufficiently large TT, we obtain the correct value for s=log⁡gxs=\log_g x with high probability. Alternatively, one can apply the continued fraction algorithm [5] to evaluate ss, see [6] for further details. Note: Alternatively, you could implement the QFTZr\text{QFT}_{\mathbb{Z}_r} over general rr, and instead of the uniform superposition, prepare the states: 1r∑x∈r∣x⟩\frac{1}{\sqrt{r}}\sum_{x\in\mathbb{r}}|x\rangle in alpha, beta. Then, again, no continued fractions postprocessing is required.

​
Example: G=Z13×G = \mathbb{Z}_{13}^\times

MODULU_NUM = 13
G_GENERATOR = 7
X_LOG_ARG = 3
ORDER = 12


@qfunc
def discrete_log(
    g: CInt,
    x: CInt,
    N: CInt,
    order: CInt,
    alpha: Output[QNum],
    beta: Output[QNum],
    func_res: Output[QNum],
) -> None:
    reg_len = ceiling(log(order, 2)) + 1
    # we define the variables with fraction places to ease the postprocessing
    allocate(reg_len, False, reg_len, alpha)
    allocate(reg_len, False, reg_len, beta)

    hadamard_transform(alpha)
    hadamard_transform(beta)

    discrete_log_oracle(g, x, N, alpha, beta, func_res)

    invert(lambda: qft(alpha))
    invert(lambda: qft(beta))


@qfunc
def main(
    alpha: Output[QNum],
    beta: Output[QNum],
    func_res: Output[QNum],
) -> None:
    discrete_log(G_GENERATOR, X_LOG_ARG, MODULU_NUM, ORDER, alpha, beta, func_res)


constraints = Constraints(max_width=23)
preferences = Preferences(optimization_level=1)
execution_preferences = ExecutionPreferences(num_shots=10000)
qmod_Z13 = create_model(
    main,
    constraints=constraints,
    preferences=preferences,
    execution_preferences=execution_preferences,
    out_file="discrete_log",
)

qprog_Z13 = synthesize(qmod_Z13)
show(qprog_Z13)
Output:

Quantum program link: https://platform.classiq.io/circuit/36zUYE60MCHFZsP62VSX5yE4xi8
  


result_Z13 = execute(qprog_Z13).result_value()
df = result_Z13.dataframe
df.head(10)
alphabetafunc_rescountprobabilitybitstring
00.00.255910.009101010100000000
10.00.251880.008800010100000000
20.00.5010870.008710101000000000
30.00.501850.008500011000000000
40.00.504850.008501001000000000
50.00.004830.008301000000000000
60.00.754830.008301001100000000
70.00.505790.007901011000000000
80.00.0012790.007911000000000000
90.00.752780.007800101100000000

​

​
Postprocessing

We now have an additional step in postprocessing. We translate each result to the closest fraction with a denominator, which is the order: 
def closest_fraction(x, denominator):
    return round(x * denominator)


df["alpha_rounded"] = closest_fraction(df.alpha, ORDER)
df["beta_rounded"] = closest_fraction(df.beta, ORDER)
df.head(10)
alphabetafunc_rescountprobabilitybitstringalpha_roundedbeta_rounded
00.00.255910.0091010101000000000.03.0
10.00.251880.0088000101000000000.03.0
20.00.5010870.0087101010000000000.06.0
30.00.501850.0085000110000000000.06.0
40.00.504850.0085010010000000000.06.0
50.00.004830.0083010000000000000.00.0
60.00.754830.0083010011000000000.09.0
70.00.505790.0079010110000000000.06.0
80.00.0012790.0079110000000000000.00.0
90.00.752780.0078001011000000000.09.0
Now, we take a sample where beta is co-prime to the order, such that we can get the logarithm by multiplying alpha by the modular inverse. If the alpha, beta registers are large enough, we are guaranteed to sample it with a good probability: 
import numpy as np


def modular_inverse(x):
    return [pow(a, -1, ORDER) for a in x]


df = df[np.gcd(df.beta_rounded.astype(int), ORDER) == 1].copy()
df["beta_inverse"] = modular_inverse(df.beta_rounded.astype("int"))
df["logarithm"] = df.alpha_rounded * df.beta_inverse % ORDER
df.head(10)
alphabetafunc_rescountprobabilitybitstringalpha_roundedbeta_roundedbeta_inverselogarithm
480.656250.093759490.0049100100011101018.01.018.0
500.343750.906251480.0048000111101010114.011.0118.0
540.656250.093757420.0042011100011101018.01.018.0
550.343750.906258420.0042100011101010114.011.0118.0
580.656250.0937510410.0041101000011101018.01.018.0
590.343750.906254400.0040010011101010114.011.0118.0
610.343750.4062512400.0040110001101010114.05.058.0
660.656250.593752380.0038001010011101018.07.078.0
670.656250.593754380.0038010010011101018.07.078.0
690.343750.4062511380.0038101101101010114.05.058.0
print(f"The descrite logarithm is: {df.logarithm[:1].iloc[0]}")
Output:

The descrite logarithm is: 8.0
To verify the results, we check whether gsmod  N=glog⁡gxmod  N=xg^s \mod N = g^{\log_g x}\mod N = x. 
assert len(df.logarithm) > 0
assert np.allclose(G_GENERATOR ** df.logarithm[:10] % MODULU_NUM, X_LOG_ARG)

​

​
Measurement Distribution Heuristic Plots

The expected measurement results are showcased by plotting the theoretical probability distributions of the α\alpha and β\beta quantum variables, for varying numbers of qubits t=5t=5 and t=8t=8 and the specific case of k=5k=5. Since r=12≠2mr = 12 \neq 2^m for some integer mm, we obtain a probability distribution characterized by a narrow peak of width ∼1/T\sim 1/T. As the number of qubits is increased, T=2tT=2^t increases, the probability distribution becomes narrower. Therefore, improving the success probability. 
MODULU_NUM = 13
X_LOG_ARG = 3

t = np.ceil(np.log2(ORDER)) + 1
T = 2**t
k = 5
r = ORDER
x_data = np.arange(T)

x0_mu = k * T / r
amplitudes_mu = np.sin(np.pi * (x_data - x0_mu)) / (
    T * np.sin(np.pi * (x_data - x0_mu) / T)
)
probabilities_mu = np.abs(amplitudes_mu) ** 2
probabilities_mu /= np.sum(probabilities_mu)


s = np.log(X_LOG_ARG) / np.log(G_GENERATOR)  # same as log_{G_GENERATOR}(X_LOG_ARG)
x0_nu = s * k * T / r
amplitudes_nu = np.sin(np.pi * (x_data - x0_nu)) / (
    T * np.sin(np.pi * (x_data - x0_nu) / T)
)
probabilities_nu = np.abs(amplitudes_nu) ** 2
probabilities_nu /= np.sum(probabilities_nu)

# Create figure and axis with custom size
fig, ax = plt.subplots(figsize=(8, 6))

# Plot data
ax.plot(x_data, probabilities_mu, label="alpha")
ax.plot(x_data, probabilities_nu, "r-.", label="beta")


# Customize axis labels and title font sizes
ax.set_xlabel("mu and nu values", fontsize=14)
ax.set_ylabel("", fontsize=14)
ax.set_title("Probabilities the Alpha and Beta Registers, t=5", fontsize=16)

# Increase tick label (axis numbers) size
ax.tick_params(axis="both", labelsize=12)


# Show legend
ax.legend(fontsize=12, loc="best")

# Display the plot
plt.show()
output 
MODULU_NUM = 13
X_LOG_ARG = 3

t = np.ceil(np.log2(ORDER)) + 4
T = 2**t
k = 5
r = ORDER
x_data = np.arange(T)

x0_mu = k * T / r
amplitudes_mu = np.sin(np.pi * (x_data - x0_mu)) / (
    T * np.sin(np.pi * (x_data - x0_mu) / T)
)
probabilities_mu = np.abs(amplitudes_mu) ** 2
probabilities_mu /= np.sum(probabilities_mu)

s = np.log(X_LOG_ARG) / np.log(G_GENERATOR)  # same as log_{G_GENERATOR}(X_LOG_ARG)
x0_nu = s * k * T / r
amplitudes_nu = np.sin(np.pi * (x_data - x0_nu)) / (
    T * np.sin(np.pi * (x_data - x0_nu) / T)
)
probabilities_nu = np.abs(amplitudes_nu) ** 2
probabilities_nu /= np.sum(probabilities_nu)
# Create figure and axis with custom size
fig, ax = plt.subplots(figsize=(8, 6))

# Plot data
ax.plot(x_data, probabilities_mu, label="alpha")
ax.plot(x_data, probabilities_nu, "r-.", label="beta")


# Customize axis labels and title font sizes
ax.set_xlabel("mu and nu values", fontsize=14)
ax.set_ylabel("", fontsize=14)
ax.set_title("Probabilities the Alpha and Beta Registers, t=8", fontsize=16)

# Increase tick label (axis numbers) size
ax.tick_params(axis="both", labelsize=12)


# Show legend
ax.legend(fontsize=12, loc="best")

# Display the plot
plt.show()
output

​
Technical notes

​

​
Equivalence with the Abelian subgroup problem

The discrete logarithm is a specific case of the Hidden Subgroup Problem (HSP) [4]. The HSP can be stated as follows: Let GG be a group and HH is subgroup of GG. We are given an (oracle) function ff with the promise that:
  1. ff is constant on the right cosets of HH.
  2. Elements of distinct cosets produce different oracle values.
That is for g∈Gg\in G if and only if x,y∈gHx,y\in g H, f(x)=f(y)f(x) = f(y). Goal: Identify a generating set for the subgroup HH; that is, a collection of elements of HH whose products yield every element of the subgroup. Note, that there is always a generating set of size Ω(log⁡(∣H∣))\Omega(\log(|H|)). The discrete logarithm problem is an instance of the Abelian HSP.
  • GG is the additive group ZN×ZN\mathbb{Z}_N\times \mathbb{Z}_N.
  • The hidden subgroup is H={(0,0),(1,−log⁡gx),(2,−2log⁡gx,…,(r−1,−(r−1)log⁡gx))} H = \{(0,0), (1,-\log_g x),(2,-2\log_g x,\dots, (r-1,-(r-1)\log_g x))\}.
  • The cosets are the of the form {(α,λ+αlog⁡xg)}\{(\alpha, \lambda + \alpha\log_x g)\}, where λ∈Zr\lambda \in \mathbb{Z}_r. It is straightforward to show that f(α,β)=gλf(\alpha,\beta) = g^\lambda for all elements of the coset.
Solution of the HSP provides a generator of (ν,−νlog⁡gx)(\nu, -\nu\log_g x), for ν\nu coprime to rr, which allows evaluating the discrete logarithm by taking the modular inverse of ν\nu: s=−ν−1νlog⁡gxs=-\nu^{-1}\nu \log_g x.

​

​
Diffie-Hellman Secret Key Sharing Protocol

The Diffie–Hellman protocol enables two parties (“Alice” and “Bob”) to establish a shared secret key, and its security against an eavesdropper (“Eve”) relies on the computational hardness of the discrete logarithm problem. The protocol for sharing a secret key includes the following steps:
  1. A prime pp and a generator gg of a multiplicative group mod pp, are published publicly.
  2. Alice chooses a number a∈Zpa\in \mathbb{Z}_p and computes A=gaA = g^a. aa is known only to Alice.
  3. Bob chooses a number b∈Zpb\in \mathbb{Z}_p and computes B=gaB = g^a. bb is known only to Bob.
  4. Alice sends Bob a message over a public channel containing AA, Bob sends Alice a message containing BB.
  5. Alice computes the secret key K=Ba=gab K = B^a = g^{ab}, and similarly Bob computes the secret key K=Ab=gabK=A^b = g^{ab}.
If the eavesdropper, Eve, can implement the discrete logarithm algorithm, she can evaluate a=log⁡g(A)a=\log_g(A) and b=log⁡g(B)b = \log_g(B), by intercepting Alice’s and Bob’s messages, AA and BB. Knowledge of aa, bb and gg allows her a straightforward calculation of K=gabK=g^{ab}. a

​
References

[1]: Discrete Logarithm (Wikipedia) [2]: Shor, Peter W. “Algorithms for quantum computation: discrete logarithms and factoring.” Proceedings 35th annual symposium on foundations of computer science. IEEE, 1994. [3]: Diffie-Hellman Key Exchange (Wikipedia) [4]: Hidden Subgroup Problem (Wikipedia) [5]: Continued Fraction (Wikipedia) [6]: Proos, J., & Zalka, C. (2003). Shor’s discrete logarithm quantum algorithm for elliptic curves. arXiv preprint quant-ph/0301141