Qiskit Community Summer Jam – North Carolina, 2020
Matthew Gregoire, Richard Howe, Guan-Wen Chou, Anton Nikulsin, and Jessica Horn
In general terms, let $G$ be a group, let $a$ be an element of $G$ of order $r$, and let $b$ be some power of $a$, so that $b = a^m$, for $0 \leq m < r$. The discrete logarithm problem is as follows: given $a$ and $b$, find the integer $m$. While $b$ may be easy to calculate given $a$ and $m$, there is no known efficient classical algorithm for finding $m$ given $a$ and $b$. This fact is the basis of many modern cryptosystems. In this project, we present a novel method for solving the discrete logarithm problem.
For our purposes, we will always be dealing with integers in a certain modulus. So our goal is to solve the equation $b = a^m \mod N$ for $m$, assuming that we know $a$, $b$, and $N$.
The half-bit of $b$, denoted $HB_a(b)$, is defined as follows.
$$ HB_a(b) = \begin{cases} 0 & 0 \leq m < r/2 \\ 1 & r/2 \leq m < r \end{cases} $$Thus the half-bit of $b$ has a natural interpretation as the most significant bit of the binary expansion of $m$. In his 1988 thesis, Burton S. Kaliski Jr. proves that, given an oracle that is able to correctly predict the half-bit with probability at least $1/2 + \varepsilon$, the discrete logarithm problem can be solved in polynomial time, by making several oracle queries on values closely related to $b$. Kaliski also presents the algorithm to solve this problem on page 107 of his thesis.
In 2017, Kaliski published a paper describing a quantum circuit to efficiently implement this oracle, which he calls a quantum "magic box." While Shor's algorithm depends on solving the entire discrete logarithm problem using one quantum circuit, implementing this oracle would only require the quantum processor to calculate one bit of the discrete logarithm at a time, potentially saving valuable computational resources.
In this project, our first and foremost goal was to implement the quantum oracle in Qiskit. We also wanted to implement Kaliski's discrete logarithm algorithm, and test it on a simple example. We were able to complete both of these goals successfully.
Below is our code for the quantum "magic box" algorithm in its entirety. Here is a brief summary of the functions implemented:
invert(x, N): returns the integer $x^{-1}$ such that $x\cdot x^{-1} = 1$ in mod $N$.create_unitary(a, N): our quantum algorithm frequently uses unitary transformations $U_{a^x}$, which transform the state $\ket{x}\ket{y}$ into the state $\ket{x}\ket{a^x y}$. This function creates a unitary matrix to perform this transformation.create_controlled_unitary(a, N): similar to create_unitary, but returns a matrix which performs $U_{a^x}$ controlled by the first input qubit.oracle(a, b, N): Builds a quantum circuit to calculate the half-bit of $b$, then runs the circuit $\lceil 8\pi\rceil$ times, as specified in the paper, and outputs the half-bit estimation. The circuit can only run on a simulator, due to hardware limitations described later.import warnings
warnings.filterwarnings("ignore", category=UserWarning)
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute, Aer
from qiskit.circuit.library import QFT
from qiskit.extensions import UnitaryGate
from numpy import log, sqrt, floor, pi, arange
import numpy as np
# Calculates the multiplicative inverse of x mod N
# (the number y such that xy = 1 (mod N)) using
# the extended Euclidean algorithm.
def invert(x, N):
q = [0, 0]
r = [N, x]
t = [0, 1]
while r[-1] != 0:
q.append(r[-2]//r[-1])
r.append(r[-2] - (q[-1]*r[-1]))
t.append(t[-2] - (q[-1]*t[-1]))
return int(t[-2] % N)
# Returns a unitary matrix which has the effect of multiplying each
# input |x> by a in mod N, resulting in the state |ax>.
def create_controlled_unitary(a, N):
dim = 2**int(np.ceil(np.log(N)/np.log(2)) + 1)
U = np.zeros((dim, dim))
# Generate a permutation of the multiplicative group of Z_N.
for i in range(int(dim/2)):
U[i,i] = 1
for i in range(N):
U[int(dim/2) + i, ((a*i) % N)+int(dim/2)] = 1
# The remaining states are irrelevant.
for i in range(N, int(dim/2)):
U[int(dim/2) + i, int(dim/2) + i] = 1
return U
def create_unitary(a, N):
a = int(np.round(a) % N)
dim = 2**int(np.ceil(np.log(N)/np.log(2)))
U = np.zeros((dim, dim))
# Generate a permutation of the multiplicative group of Z_N.
for i in range(N):
U[i, ((a*i) % N)] = 1
# The remaining states are irrelevant.
for i in range(N, dim):
U[i, i] = 1
return U
# b is some power of a, and the oracle outputs the half-bit of m,
# where b = a^m (mod N), with >50% probability.
def oracle(a, b, N, verbose=False, draw=False):
# Calculate the order of a
r = 1
product = a
while product != 1:
product = (product * a) % N
r += 1
# Find number of bits needed to store a value from 0 to N-1 (n)
# and initialize 2 quantum registers of size n.
n = int(np.ceil(np.log(N)/np.log(2)))
qr1, qr2 = QuantumRegister(n), QuantumRegister(n)
cr1, cr2 = ClassicalRegister(n), ClassicalRegister(1)
qc = QuantumCircuit(qr1, qr2, cr1, cr2)
#Change second register to state |00...01>
qc.x(qr2[n-1])
#Add H gate to first register
for i in range(n):
qc.h(qr1[i])
# We need log_2(n) different matrices U_(a^(2^x))
for i in range(n):
U = create_controlled_unitary(a**(2**(n-i)) % N, N)
qubits = [qr1[i]] + [qr2[j] for j in range(n)]
qc.iso(U, qubits, [])
qc.append(QFT(n), [qr1[i] for i in range(n)])
for i in range(n):
qc.measure(qr1[i], cr1[i])
# Reuse the first quantum register, because we don't need it anymore.
for i in range(2**(n-1), 2**n):
qc.x(qr1[0]).c_if(cr1, i)
qc.h(qr1[0])
qc.barrier()
# Now cr1 is in state y. We define k to be the closest integer to y*r / 2**n.
# I don't think there's any way to get the result of the measurement mid-circuit
# in qiskit. So this is a stopgap measure for now.
for y in range(2**n):
k = int(np.round(y*r/(2**n))) % r
kInv = bin(invert(k, r))[2:]
# Pad kInv with initial zeros
while len(kInv) < n:
kInv = '0' + kInv
if '1' in kInv:
for i in range(len(kInv)):
bit = int(kInv[i])
if bit == 1:
# Apply U operation only if the value of cr1 is y.
U = create_unitary(b**(2**i) % N, N)
qc.iso(U, [qr2[i] for i in range(n)], []).c_if(cr1, y)
qc.barrier()
qc.rz(-np.pi/2 , qr1[0])
qc.h(qr1[0])
qc.measure(qr1[0], cr2[0])
if draw:
print(qc.draw(output='text'))
# Execute the circuit on a simulator.
backend_name = 'qasm_simulator'
backend = Aer.get_backend(backend_name)
if verbose:
print("Running circuit on", backend_name, "...")
result = execute(qc, backend, shots=int(np.ceil(8*np.pi))).result().get_counts(qc)
if verbose:
print(result)
# See whether the half-bit was measured more often as 0 or as 1, and
# return that value as the half-bit.
zeros_count = 0
ones_count = 0
for k in result.keys():
half_bit = k[0]
if half_bit == '0':
zeros_count += result[k]
else:
ones_count += result[k]
if verbose:
print('Zeros:', zeros_count, '\tOnes:', ones_count)
if zeros_count > ones_count:
return 0
else:
return 1
For example, here is the circuit associated with the oracle of a small half-bit calculation. The circuit below calculates $HB(2)$, in mod $5$ with base $3$. We can verify that $3$ has multiplicative order $4$, and $3^3 = 2 \;\mathrm{mod}\; 5$. Since $4/2 \leq 3 < 4$, we can see that $HB(2) = 1$.
oracle(3, 2, 5, draw=True)
print("The oracle predicts that HB(2) =", oracle(3, 2, 5))
Using the oracle above, the following code implements the algorithm Logarithm as specified on page 107 of Kaliski's thesis. This code also uses the notational conventions of that paper, for easy reference.
calculate_order(a, N) returns the order of $a$ in modulus $N$. In the previous notation, it returns the value of $r$.Before continuing, note that the below implementation of calculate_order actually solves the discrete logarithm problem! It calculates the order of $a$ by brute force repeated multiplication, which will cycle through all possible powers of $a$. This is fine for the examples we provide, because the numbers involved are so small. However, in almost all practical applications the order of $a$ is known beforehand, so there is no need to perform this calculation.
n1(c, n) is a straightforward helper function from Kaliski's thesis.EstimateCrossCorrelation(G, X, n, l, d, period) is another helper function, which transforms the discrete logarithm problem into a question about signal processing.Logarithm(G, X, n) solves the discrete logarithm problem. That is, given that $G^m = X \mod n$, this algorithm solves for $m$. It does so by making an initial guess for the logarithm, and then refining the guess through repeated calls to the oracle.from random import randint, random
e = 1/pi # Maximum epsilon from the Magic Box paper.
def calculate_order(a, N):
power = 1
curr = a
while (curr != 1):
curr = (curr * a) % N
power += 1
return power
'''
Determines the most significant bit of an integer c
@param c: Base
@param n: Modulus
@return: 1 or -1
'''
def n1(c, n):
val = -1
if c % n >= (n / 2):
val = 1
return val
'''
Decides which 'guess' best fits with the output of the oracle using cross correlation.
@param G: Base of the logarithm.
@param X: Power of the logarithm.
@param n: Modulus.
@param l: Lag of cross correlation - refers to how far the series are offset.
@param d: The probability that the estimation will be incorrect.
'''
def EstimateCrossCorrelation(G, X, n, l, d, period):
total = 0 # Called 'sum' in algorithm
# Compute the number of trials
m = int(round(2 / (sqrt(d) * e)))
# Compute estimate
for trial in range(m):
t = randint(1, n)
output = oracle(G, X, n) # Oracle(base, LHS, mod)
if output == 0:
output = -1
total = total + (output*(X * (G ** period)) * n1(t + l, n))
return (total / m)
'''
Main algorithm for computing the discrete logarithm.
@param G: Base of a the logarithm.
@param X: Power of the logarithm.
@param n: Modulus
'''
def Logarithm(G, X, n):
step = (n * e)
l = int(log(n)/log(2)) # Number of iterations
d = round((l / 4), 10) # Limit on probability error
period = calculate_order(G, n)
# Repeat until logarithm is found.
while True:
# Make initial guess.
for c in arange(0, n-1+step, step):
# Refine guess.
for i in range(l-1, -1, -1):
Xi = (X ** (2**i)) % n
if (EstimateCrossCorrelation(G, Xi, n, c/2, d, period)
> EstimateCrossCorrelation(G, Xi, n, c/2 + n/2, d, period)):
c = (c/2) % period
else:
c = (c/2 + n/2) % period
potential = int(floor(c))
if ((G ** potential) % n) == X:
return potential
Now we can calculate the discrete logarithm with just one call to the above function! Below, we calculate the value $m$ such that $7^m = 13 \mod 15$.
import time
start = time.time()
print(Logarithm(7, 13, 15))
end = time.time()
print("Execution took", int(round(end - start)), "seconds.")
Unfortunately, this algorithm is incredibly inefficient. There are a few reasons as to why this is the case.
The quantum oracle is structured in two phases. The first phase operates on two quantum registers, and results in (a) the measured value of register 1, and (b) the superposition on register 2. The measured value is then used to construct the rest of the circuit for the second phase of the quantum algorithm.
Unfortunately, Qiskit does not yet have a way to access the value in a classical register mid-circuit. Therefore, our algorithm has no way to know what value was measured in the first phase. Our stopgap solution is to build every possible circuit that could result from the first phase, and only execute them if the classical register is in a particular state. This starts on the line for y in range(2**n), which is clearly exponential in the input size.
In addition, Kaliski's 2017 paper relies heavily on the unitary transformations $U_{a^x}$, as defined above. However, we were unable to find an efficient way to build these unitary transformations. The current solution is to build them as needed, although the matrices used to represent these transformations have size exponential in the number of qubits used.
The discrete logarithm algorithm is directly from Kaliski's 1988 thesis, and is provably fast. Specifically, it requires $O\left(\log^{5/2}\left(n/\epsilon^3\right)\right)$ modular multiplications, and $O\left(\log^{3/2}\left(n/\epsilon^2\right)\right)$ calls to the oracle. Given the inefficiencies in the quantum oracle, this results in a significant speed disadvantage.
In addition, our quantum oracle relies on features that are currently only accessible on simulators, and not real quantum hardware. Specifically, the algorithm requires classically controlled quantum gates, which are not available on any current IBMQ processors.
However, if in the future IBM releases quantum hardware that supports modifying circuits on-the-fly, we could feasibly eliminate all classically controlled operations from our circuit. This would also significantly improve the efficiency of our oracle algorithm.
There is still much to be done on this project. Most pressingly, we would like to find an efficient implementation of the unitaries $U_{a^x}$ that the oracle algorithm depends on. In addition, a large amount of inefficiency comes from our implementation of the second phase of the quantum algorithm. While the best solution to this problem is additional functionality in Qiskit and quantum hardware, we could find other workarounds. For example, the QuantumCircuit.snapshot() method could greatly increase our code's efficiency on quantum simulators. And because our algorithm can only run on simulators for now, this could be a viable near-term solution. Kaliski's 2017 paper also describes a number of potential improvements to the quantum oracle, any one of which has the potential to improve its performance.
Finally, we would like to build off the algorithm presented here in the future. This is a raw mathematical model for solving the discrete logarithm problem, but this has many practical applications, most notably in the field of cryptography. In an attempt to make quantum algorithms more accessible, we would like to expand this notebook to explain in detail Kaliski's algorithm, and present a basic example of this technique being used to break RSA encryption. (Of course, the keys would likely be on the order of 3-4 bits!)
A huge thank you to IBM Quantum, the Qiskit community, and the organizers of the North Carolina Qiskit Community Summer Jam 2020 for making this project possible!
Kaliski Jr., B. (2017). A Quantum “Magic Box” for the Discrete Logarithm Problem. Cryptology EPrint Archive, Report 2017/745. https://eprint.iacr.org/2017/745
Kaliski Jr., B. (1988). Elliptic Curves and Cryptography: A Pseudorandom Bit Generator and Other Tools (Doctoral dissertation, MIT, Cambridge, USA). Retrieved from https://dspace.mit.edu/bitstream/handle/1721.1/14709/18494044-MIT.pdf
Mosca M., Ekert A. (1999) The Hidden Subgroup Problem and Eigenvalue Estimation on a Quantum Computer. In: Williams C.P. (eds) Quantum Computing and Quantum Communications. QCQC 1998. Lecture Notes in Computer Science, vol 1509. Springer, Berlin, Heidelberg