Discrete Logarithm
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.
Formulating the Problem
-
Input: A cyclic group \(G = \langle g \rangle\) with \(g\) as a generator, and an element \(x\in G\).
-
Promise: There is a number \(s\) such that \(g^s = x\).
-
Output: \(s\), the discrete logarithm: \(s = \log_gx\).
In Shor's implementation, the order of \(g\) is assumed to be known beforehand (for example, using the order finding algorithm). We also assume it in the demonstration.
The Discrete Log problem is a specific example of the Abelian Hidden Subgroup Problem [4] for the case of an additive group, with this function:
Building the Algorithm with Classiq
The heart of the algorithm's logic is the implementation of the function
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_exp function accepts these arguments:
-
n: CInt- modulo number -
a: CInt- base of the exponentiation -
x: QArray[QBit]- unsigned integer to multiply by the exponentiation -
power: QArray[QBit]- power of the exponentiation
So the function implements \(|power\rangle|x\rangle \rightarrow |power\rangle|x \cdot a ^ {power}\mod n\rangle\).
from classiq import *
@qfunc
def discrete_log_oracle(
g_generator: CInt,
x_element: CInt,
N_modulus: CInt,
x1: QArray[QBit],
x2: QArray[QBit],
func_res: Output[QNum],
) -> None:
allocate(ceiling(log(N_modulus, 2)), func_res)
func_res ^= 1
modular_exp(N_modulus, x_element, func_res, x1)
modular_exp(N_modulus, g_generator, func_res, x2)
Full Algorithm
-
Prepare uniform superposition over the first two quantum variables
x1,x2. Each variable has size \(\lceil \log r\rceil + \log({1/{\epsilon}})\). In the special case where \(r\) is a power of 2, \(\log r\) is enough. -
Compute
discrete_log_oracleon thefunc_resvariable.func_resis of size \(\lceil \log N\rceil\). -
Apply the inverse Fourier transform
x1,x2. -
Measure.
from classiq.qmod.symbolic import ceiling, log
@qfunc
def discrete_log(
g: CInt,
x: CInt,
N: CInt,
order: CInt,
x1: Output[QArray[QBit]],
x2: Output[QArray[QBit]],
func_res: Output[QArray[QBit]],
) -> None:
reg_len = ceiling(log(order, 2))
allocate(reg_len, x1)
allocate(reg_len, x2)
hadamard_transform(x1)
hadamard_transform(x2)
discrete_log_oracle(g, x, N, x1, x2, func_res)
invert(lambda: qft(x1))
invert(lambda: qft(x2))
After the inverse QFTs (under the assumption of \(r=2^m\) for some \(m\)), the variables become
\(\(|\psi\rangle = \sum_{\nu\in\mathbb{Z}_r, \delta\in G}\omega^{\nu\delta}|\nu\cdot log_gx\rangle_{x_1}|\nu\rangle_{x_2}|\delta>_{func\_res}\)\) .
For every \(\nu\) that has a multiplicative inverse in \(\mathbb{Z}_r\), we can extract \(s=\log_xg\) by multiplying the first variable result by its inverse.
If \(r\) is not a power of 2, the variables get an approximation of |\(\log_g(x)\cdot \nu/ r\rangle_{x_1} |\nu / r\rangle_{x_2}\). So we can use the continued fractions algorithm [5] to compute \(\nu/r\), then use the same technique to calculate \(\log_gx\).
Note: Alternatively, you could implement the \(QFT_{\mathbb{Z}_r}\) over general \(r\), and instead of the uniform superposition, prepare the states: \(\frac{1}{\sqrt{r}}\sum_{x\in\mathbb{r}}|x\rangle\) in x1, x2. Then, again, no continued fractions postprocessing is required.
Example: \(G = \mathbb{Z}_5^\times\)
For this specific demonstration, we choose \(G = \mathbb{Z}_5^\times\), with \(g=3\) and \(x=2\). With this setting, \(log_gx=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.
MODULU_NUM = 5
G_GENERATOR = 3
X_LOGARITHM = 2
ORDER = MODULU_NUM - 1 # as 5 is prime
@qfunc
def main(
x1: Output[QNum],
x2: Output[QNum],
func_res: Output[QNum],
) -> None:
discrete_log(G_GENERATOR, X_LOGARITHM, MODULU_NUM, ORDER, x1, x2, func_res)
from classiq.execution import ExecutionPreferences
constraints = Constraints(max_width=13)
qmod_Z5 = create_model(
main,
constraints=constraints,
execution_preferences=ExecutionPreferences(num_shots=4000),
out_file="discrete_log",
)
qprog_Z5 = synthesize(qmod_Z5)
show(qprog_Z5)
result_Z5 = execute(qprog_Z5).result_value()
result_Z5.parsed_counts
Note that func_res is uncorrelated to the other variables, and we get uniform distribution, as expected.
We take only the x2 that are co-prime to \(r=4\), so they have a multiplicative-inverse. Hence x2=1,3 are the relevant results.
So we get two relevant results (for all different \(\delta\)s): \(|1\rangle|3\rangle\), \(|3\rangle|1\rangle\). All that remains to get the logarithm is to multiply x1 by the inverse of x2:
print(3 * pow(1, -1, 4))
3
print(1 * pow(3, -1, 4))
3
Verify we received the correct discrete logarithm:
print(3**3 % 5)
2
And, indeed, both cases give the same result, which is exactly the discrete logarithm: \(\log_32 \mod 5 = 3\).
Example: \(G = \mathbb{Z}_{13}^\times\)
In this case the order is not a power of 2. During circuit creation, we change the state preparation: instead of creating the entire uniform distribution on the x1, x2 variables, we load them with the uniform superposition of only the first #ORDER states.
We do that using the prepare_uniform_trimmed_state library function, which efficiently prepares such a state.
MODULU_NUM = 13
G_GENERATOR = 7
X_LOGARITHM = 3
ORDER = 12
@qfunc
def discrete_log(
g: CInt,
x: CInt,
N: CInt,
order: CInt,
x1: Output[QNum],
x2: 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_num(reg_len, False, reg_len, x1)
allocate_num(reg_len, False, reg_len, x2)
prepare_uniform_trimmed_state(ORDER, x1)
prepare_uniform_trimmed_state(ORDER, x2)
discrete_log_oracle(g, x, N, x1, x2, func_res)
invert(lambda: qft(x1))
invert(lambda: qft(x2))
@qfunc
def main(
x1: Output[QNum],
x2: Output[QNum],
func_res: Output[QNum],
) -> None:
discrete_log(G_GENERATOR, X_LOGARITHM, MODULU_NUM, ORDER, x1, x2, 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,
out_file="discrete_log_large",
execution_preferences=execution_preferences,
)
qprog_Z13 = synthesize(qmod_Z13)
show(qprog_Z13)
result_Z13 = execute(qprog_Z13).result_value()
result_Z13.parsed_counts[:10]
[{'x1': 0.65625, 'x2': 0.34375, 'func_res': 8.0}: 19,
{'x1': 0.0, 'x2': 0.5, 'func_res': 3.0}: 19,
{'x1': 0.34375, 'x2': 0.65625, 'func_res': 4.0}: 18,
{'x1': 0.3125, 'x2': 0.40625, 'func_res': 7.0}: 17,
{'x1': 0.6875, 'x2': 0.34375, 'func_res': 2.0}: 16,
{'x1': 0.0, 'x2': 0.0, 'func_res': 6.0}: 16,
{'x1': 0.65625, 'x2': 0.84375, 'func_res': 3.0}: 16,
{'x1': 0.6875, 'x2': 0.59375, 'func_res': 12.0}: 16,
{'x1': 0.65625, 'x2': 0.0625, 'func_res': 7.0}: 15,
{'x1': 0.0, 'x2': 0.5, 'func_res': 12.0}: 15]
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)
parsed = []
for r in result_Z13.parsed_counts[:30]:
x1 = closest_fraction(r.state["x1"], ORDER)
x2 = closest_fraction(r.state["x2"], ORDER)
rounded = {"x1": x1, "x2": x2, "shots": r.shots}
parsed.append(rounded)
print(rounded)
{'x1': 8, 'x2': 4, 'shots': 19}
{'x1': 0, 'x2': 6, 'shots': 19}
{'x1': 4, 'x2': 8, 'shots': 18}
{'x1': 4, 'x2': 5, 'shots': 17}
{'x1': 8, 'x2': 4, 'shots': 16}
{'x1': 0, 'x2': 0, 'shots': 16}
{'x1': 8, 'x2': 10, 'shots': 16}
{'x1': 8, 'x2': 7, 'shots': 16}
{'x1': 8, 'x2': 1, 'shots': 15}
{'x1': 0, 'x2': 6, 'shots': 15}
{'x1': 4, 'x2': 8, 'shots': 15}
{'x1': 0, 'x2': 0, 'shots': 15}
{'x1': 0, 'x2': 12, 'shots': 15}
{'x1': 4, 'x2': 8, 'shots': 15}
{'x1': 0, 'x2': 3, 'shots': 15}
{'x1': 4, 'x2': 8, 'shots': 15}
{'x1': 0, 'x2': 6, 'shots': 14}
{'x1': 8, 'x2': 1, 'shots': 14}
{'x1': 8, 'x2': 7, 'shots': 14}
{'x1': 0, 'x2': 6, 'shots': 14}
{'x1': 0, 'x2': 9, 'shots': 14}
{'x1': 4, 'x2': 5, 'shots': 14}
{'x1': 4, 'x2': 8, 'shots': 13}
{'x1': 8, 'x2': 10, 'shots': 13}
{'x1': 0, 'x2': 6, 'shots': 13}
{'x1': 0, 'x2': 6, 'shots': 13}
{'x1': 0, 'x2': 9, 'shots': 13}
{'x1': 8, 'x2': 7, 'shots': 13}
{'x1': 4, 'x2': 2, 'shots': 13}
{'x1': 0, 'x2': 3, 'shots': 13}
Now, we take a sample where x2 is co-prime to the order, such that we can get the logarithm by multiplying x1 by the modular inverse. If the x1, x2 registers are large enough, we are guaranteed to sample it with a good probability:
logarithm = None
for sample in parsed:
try:
inverse = pow(sample["x2"], -1, ORDER)
logarithm = sample["x1"] * inverse % ORDER
print(sample)
print("Inverse:", inverse, "Logarithm:", logarithm)
break
except:
continue
{'x1': 4, 'x2': 5, 'shots': 17}
Inverse: 5 Logarithm: 8
assert logarithm is not None
assert G_GENERATOR**logarithm % MODULU_NUM == X_LOGARITHM
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)