Use this file to discover all available pages before exploring further.
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⟩ with a generator g∈G, an element x∈G and the order of the group is known r=∣G∣. A quantum oracle Uf applying the transformation
Uf∣α⟩∣β⟩∣0⟩=∣α⟩∣β⟩∣f(α,β)⟩,where f(α,β)=xαgβ with α,β∈Zr.
Promise: There is a positive integer s∈N such that gs=x.
Output: The least positive integer s=loggx, i.e, the discrete logarithm.
Complexity: A single application of the algorithm succeeds with probability O(loglogr) and requires O((logt)2) time, where t=⌈logr+log(1/ϵ)⌉ and 1−ϵ is success probability. Hence, boosting the success probability to a constant O(1) via repetition yields a total expected runtime of
O((logt)2loglogr).
Keywords: Abelain Hidden subgroup problem, quantum Fourier transform, period finding.
We first consider the case where r=2m and later generalize.
We begin by introducing an input state∣ψ0⟩=∣0m⟩∣0m⟩∣0R⟩,where R is the number of bits required to represent f‘s image. In addition, we note that f(α,β)=xαgβ=gαloggx+β is constant on lines satisfying{(α,β)∈Zr2,αloggx+β=λmodr}.The algorithm includes four main steps, and has a similar structure to the other quantum algorithms for the hidden subgroup problem.
We first prepare a uniform superposition over the input states of the first two registers:
∣0m⟩∣0m⟩∣0m⟩H⊗mr1α,β∈Zr∑∣α⟩∣β⟩∣0m⟩.
Perform a query to the oracle
Ufr1α,β∈Zr∑∣α⟩∣β⟩∣f(α,β)⟩.Utilizing the periodicity of f, we can express the state as=r1α,λ∈Zr∑∣α⟩∣λ−αloggx⟩∣gλ⟩.
Perform a double inverse Fourier transform over the group Zr
QFTZr†×QFTZr†r21λ,μ,ν∈Zr∑(α∑e−i2πα(μ−νloggx)/r)e−i2πλν/r∣μ⟩∣ν⟩∣gλ⟩=r1λ,ν∈Zr∑ei2πλν/r∣νloggx⟩∣ν⟩∣gλ⟩,where the Fourier transform over Zr is QFTZr∣α⟩=r1∑μ∈Zrei2παμ/r∣μ⟩ and in the second line we utilized the identity ∑αei2παη=rδ0η (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λ⟩ with a uniform distribution.After the collapse, the resulting state is independent of λ, therefore, we can discard this measurement outcome.From the measurement of the second quantum variable, we obtain with a uniform distribution over ν the outcome νloggx.With a probablility of order O(loglogr), the obtained result ν is co-prime to r, and there exists a modular r multiplicative inverse ν−1.Under this condition, we can multiply the outcome of the second register to obtain the discrete log s=loggx. Hence, repeating the experiment O(loglogr) times leads to a success probability of order O(1).
import matplotlib.pyplot as pltimport numpy as npfrom classiq import *from classiq.qmod.symbolic import ceiling, log
We consider two exemplary cases. First, we study the group Z5×, where the order is a power of r=4=22.For this case, quantum variables with only logr=m qubits are required. Following, the group Z13× exemplifies the alternative case where the order r=2m for an integer m.For these examples, we denote the modulus by N. Since, these are cyclic groups, we have N=r+1, therefore R=⌈logN⌉. In the following examples, the third register, containing the function f(α,β) after the second step, is represented by ⌈logN⌉.The heart of the algorithm’s logic is the implementation of the function∣α⟩∣β⟩∣0⟩→∣α⟩∣β⟩∣xαgβ⟩.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⋅apwmodN⟩.
After the inverse QFTs (under the assumption of r=2m for some m, therefore t=2r), the variables become∣ψ⟩=r1λ,ν∈Zr∑ei2πλν/r∣νloggx⟩α∣ν⟩β∣gλ⟩reswhere we added a subscript to the quantum variables for clarity, and ∣⟩res designates the func_res variable.If ν∈Zr has an inverse, ν−1∈Zr we can extract s=logxg from variable α by multiplying by ν−1.See the second example for the general case, for which r=2m.
For this specific demonstration, we choose G=Z5×, with g=3 and x=2.With this setting, loggx=3.We choose this specific example because the order of the group r=4 is a power of 2, 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.
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=4, so they have a multiplicative-inverse.Hence beta=1,3 are the relevant results.
So we get two relevant results (for all different λs): ∣1⟩∣3⟩, ∣3⟩∣1⟩.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:
For an arbitrary r∈N, we employ two t=⌈logr⌉+log(1/ϵ) qubit quantum variables and a third R-qubit register, initialized to the state ∣ψ0⟩=∣0t⟩∣0t⟩∣0R⟩, where R=⌈logN⌉.The derivation follows a similar structure, where the initial sums are now α,β∈ZT, with T=2t, while the period of f remains r, i.e., λ∈Zr. Therefore, after the second step we obtain the stateT1α,β∈ZT∑∣α⟩∣β⟩∣0R⟩UfT1α,β∈ZT∑∣α⟩∣β⟩∣gαloggx+β⟩.Next, we express the periodic state in an alternative form. We introduce the operator Ug, which satisfies Ug∣h⟩=∣hg⟩, where e is the identity element of the group.Therefore ∣gλ⟩=Ugλ∣e⟩.The state ∣e⟩ can be expressed as a uniform superposition of the eigenstates of Ug:∣e⟩=r1k=0∑r−1∣Ψk⟩,where ∣Ψk⟩=r1∑k=0r−1ei2πkk′/r∣gk′⟩, which satisfy Ug∣Ψk⟩=ei2πk/r∣Ψk⟩.These relations allow writing the state after the oracle operation asT1α,β∈ZT∑∣α⟩∣β⟩k=0∑r−1ei2π(αloggx+β)k/r∣Ψk⟩.Applying the inverse quantum Fourier transform leads to the stateQFTT†×QFTT†T1k=0∑r−1μ,ν∈ZT∑(α∈ZT∑ei2πα(loggxk/r−μ/T)∣μ⟩)β∈ZT∑ei2πβ(k/r−ν/T)∣ν⟩∣Ψk⟩.Utilizing the geometric sum, one can show that the functions are peaked (with a width O(1/T)) around μ≈Tloggxk/r and ν≈Tk/r, correspondingly.
To evaluate the values of kloggx and k we measure alpha and beta registers, multiply by r/N and round.For sufficiently large T, we obtain the correct value for s=loggx with high probability. Alternatively, one can apply the continued fraction algorithm [5] to evaluate s, see [6] for further details.Note: Alternatively, you could implement the QFTZr over general r, and instead of the uniform superposition, prepare the states: r1∑x∈r∣x⟩ in alpha, beta. Then, again, no continued fractions postprocessing is required.
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 npdef 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 % ORDERdf.head(10)
alpha
beta
func_res
count
probability
bitstring
alpha_rounded
beta_rounded
beta_inverse
logarithm
48
0.65625
0.09375
9
49
0.0049
10010001110101
8.0
1.0
1
8.0
50
0.34375
0.90625
1
48
0.0048
00011110101011
4.0
11.0
11
8.0
54
0.65625
0.09375
7
42
0.0042
01110001110101
8.0
1.0
1
8.0
55
0.34375
0.90625
8
42
0.0042
10001110101011
4.0
11.0
11
8.0
58
0.65625
0.09375
10
41
0.0041
10100001110101
8.0
1.0
1
8.0
59
0.34375
0.90625
4
40
0.0040
01001110101011
4.0
11.0
11
8.0
61
0.34375
0.40625
12
40
0.0040
11000110101011
4.0
5.0
5
8.0
66
0.65625
0.59375
2
38
0.0038
00101001110101
8.0
7.0
7
8.0
67
0.65625
0.59375
4
38
0.0038
01001001110101
8.0
7.0
7
8.0
69
0.34375
0.40625
11
38
0.0038
10110110101011
4.0
5.0
5
8.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 gsmodN=gloggxmodN=x.
The expected measurement results are showcased by plotting the theoretical probability distributions of the α and β quantum variables, for varying numbers of qubits t=5 and t=8 and the specific case of k=5.Since r=12=2m for some integer m, we obtain a probability distribution characterized by a narrow peak of width ∼1/T.
As the number of qubits is increased, T=2t increases, the probability distribution becomes narrower. Therefore, improving the success probability.
MODULU_NUM = 13X_LOG_ARG = 3t = np.ceil(np.log2(ORDER)) + 1T = 2**tk = 5r = ORDERx_data = np.arange(T)x0_mu = k * T / ramplitudes_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) ** 2probabilities_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 / ramplitudes_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) ** 2probabilities_nu /= np.sum(probabilities_nu)# Create figure and axis with custom sizefig, ax = plt.subplots(figsize=(8, 6))# Plot dataax.plot(x_data, probabilities_mu, label="alpha")ax.plot(x_data, probabilities_nu, "r-.", label="beta")# Customize axis labels and title font sizesax.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) sizeax.tick_params(axis="both", labelsize=12)# Show legendax.legend(fontsize=12, loc="best")# Display the plotplt.show()
MODULU_NUM = 13X_LOG_ARG = 3t = np.ceil(np.log2(ORDER)) + 4T = 2**tk = 5r = ORDERx_data = np.arange(T)x0_mu = k * T / ramplitudes_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) ** 2probabilities_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 / ramplitudes_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) ** 2probabilities_nu /= np.sum(probabilities_nu)# Create figure and axis with custom sizefig, ax = plt.subplots(figsize=(8, 6))# Plot dataax.plot(x_data, probabilities_mu, label="alpha")ax.plot(x_data, probabilities_nu, "r-.", label="beta")# Customize axis labels and title font sizesax.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) sizeax.tick_params(axis="both", labelsize=12)# Show legendax.legend(fontsize=12, loc="best")# Display the plotplt.show()
The discrete logarithm is a specific case of the Hidden Subgroup Problem (HSP) [4].The HSP can be stated as follows:Let G be a group and H is subgroup of G. We are given an (oracle) function f with the promise that:
f is constant on the right cosets of H.
Elements of distinct cosets produce different oracle values.
That is for g∈G if and only if x,y∈gH, f(x)=f(y).Goal: Identify a generating set for the subgroup H; that is, a collection of elements of H whose products yield every element of the subgroup.Note, that there is always a generating set of size Ω(log(∣H∣)).The discrete logarithm problem is an instance of the Abelian HSP.
G is the additive group ZN×ZN.
The hidden subgroup is H={(0,0),(1,−loggx),(2,−2loggx,…,(r−1,−(r−1)loggx))}.
The cosets are the of the form {(α,λ+αlogxg)}, where λ∈Zr. It is straightforward to show that f(α,β)=gλ for all elements of the coset.
Solution of the HSP provides a generator of (ν,−νloggx), for ν coprime to r, which allows evaluating the discrete logarithm by taking the modular inverse of ν: s=−ν−1νloggx.
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:
A prime p and a generator g of a multiplicative group mod p, are published publicly.
Alice chooses a number a∈Zp and computes A=ga. a is known only to Alice.
Bob chooses a number b∈Zp and computes B=ga. b is known only to Bob.
Alice sends Bob a message over a public channel containing A, Bob sends Alice a message containing B.
Alice computes the secret key K=Ba=gab, and similarly Bob computes the secret key K=Ab=gab.
If the eavesdropper, Eve, can implement the discrete logarithm algorithm, she can evaluate a=logg(A) and b=logg(B), by intercepting Alice’s and Bob’s messages, A and B.Knowledge of a, b and g allows her a straightforward calculation of K=gab.a