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.

Problem formulation

• 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 will also assume it in the demonstration.

The Discrete Log problem is a specific example for the Abelian Hidden Subgroup Problem [4], for the case of an additive group, with the function:

$f: \mathbb{Z}_N \times \mathbb{Z}_N \rightarrow G$
$f(\alpha, \beta) = x^\alpha g^\beta$

How to build the Algorithm with Classiq

The heart of the algorithm's logic is the implementation of the function:

$|x_1\rangle|x_2\rangle|1\rangle \rightarrow |x_1\rangle|x_2\rangle|x^{x_1} g^{x_2}\rangle$

This is done using 2 applications of the modular exponentiation function, which was described in detail in the Shor's Factoring Algorithm notebook. So here we will just import it from the classiq's library.

The function modular_exp accepts the following arguments:

• n: CInt - modulo number

• a: CInt - base of the exponentiation

• x: QArray[QBit] - unsigned integer to multiply be the exponentiation

• power: QArray[QBit]- power of the exponentiation

So that 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[QArray[QBit]],
) -> None:

allocate(ceiling(log(N_modulus, 2)), func_res)

inplace_prepare_int(1, func_res)
modular_exp(N_modulus, x_element, func_res, x1)
modular_exp(N_modulus, g_generator, func_res, x2)


The full algorithm:

1. Prepare uniform superposition over the first 2 quantum variables x1, x2. Each variable should be with size $$\lceil \log r\rceil + \log({1/{\epsilon}})$$. In the special case where $$r$$ is a power of 2, $$\log r$$ is enough.
2. Compute discrete_log_oracle on the func_res variable. func_res should be of size $$\lceil \log N\rceil$$.
3. Apply inverse Fourier transform x1, x2.
4. 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)

discrete_log_oracle(g, x, N, x1, x2, func_res)

invert(lambda: qft(x1))
invert(lambda: qft(x2))


After the inverse QFTs, we get in the variables (under the assumption of $$r=2^m$$ for some $$m$$):

$|\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 mutplicative inverse in $$\mathbb{Z}_r$$, we can extract $$s=\log_xg$$ by multiplying the first variable result by its inverse.

In the case where $$r$$ is not a power of 2, we get in the variables and 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 using the same technique to calculate $$\log_gx$$.

Note: Alternatively, one might 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 post-process 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$$. Under this setting, $$log_gx=3$$.

We choose this specific example as the order of of the group $$r=4$$ is a power of $$2$$, and so we can get exactly the discrete logarithm, without continued-fractions post processing. In other cases, one has to use larger quantum variable for the exponents so the continued fractions post-processing will converge.

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)

Opening: http://localhost:4200/circuit/a1507b50-7312-478f-ad92-0589a85ebf79?version=0.0.0

result_Z5 = execute(qprog_Z5).result_value()
result_Z5.parsed_counts


Notice 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 2 relevant results (for all different $$\delta$$s): $$|1\rangle|3\rangle$$, $$|3\rangle|1\rangle$$. All left to do 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 got the correct discrete logarithm:

print(3**3 % 5)

2


And indeed in both cases the same result, which is exactly the discrete logarithm: $$\log_32 \mod 5 = 3$$

Example: $$G = \mathbb{Z}_{13}^\times$$

Here we take the case were the order is not a power of 2. In the circuit creation, we will need to change the state preparation. Instead of creating the entire uniform distribution on the x1, x2 variables, we will load them with the uniform superposition of only the first #ORDER states.

In order to do that we can use the library function prepare_uniform_trimmed_state which prepare such a state efficiently.

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 in order to ease the post-processing
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)
execution_preferences = ExecutionPreferences(num_shots=10000)
qmod_Z13 = create_model(
main,
constraints=constraints,
out_file="discrete_log_large",
execution_preferences=execution_preferences,
)

qprog_Z13 = synthesize(qmod_Z13)
show(qprog_Z13)

Opening: http://localhost:4200/circuit/3c990661-c569-4c99-b534-225f30334935?version=0.0.0

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]


Post process

We now have additional step in post-process. We translate each result to the closest fraction with 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 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 the x1, x2 registers are large enough, we are guaranteed to sample it with a good probability:

import math

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