diff --git a/main.sage b/main.sage
new file mode 100644
index 0000000000000000000000000000000000000000..3df4421bde94b57a7721e13d9849d6211694e58d
--- /dev/null
+++ b/main.sage
@@ -0,0 +1,130 @@
+from sage.all import *
+from sage.crypto.lwe import samples, LWE
+
+from tools import *
+from plots import *
+
+def nistlike(n, q):
+ """
+ Given n, q, generate a NIST-like LWE instance
+
+ args:
+ n (int): dimension
+ q (int): modulus
+
+ returns:
+ A (matrix): matrix
+ t (array): target vector
+ s (array): secret vector
+ e (array): error vector
+ """
+ m=1
+ DB = CenteredBinomialDistibution(2)
+ lwe_instance = LWE(n=n, q=q, D=DB, secret_dist='noise', m=m)
+ s = lwe_instance._LWE__s
+ s = array([int(x) for x in s])
+ s = array([int(x) if x<= q//2 else int(x-q) for x in s], dtype=int)
+ a = samples(m, n, lwe_instance, balanced=True)
+ a = a[0][0]
+ assert all(-q//2 <=x <= q // 2 for x in a)
+ A = zeros((n, n), dtype=int)
+ for i in range(n):
+ for j in range(n):
+ if j > i:
+ A[i][j] = int(-a[i-j])
+ else:
+ A[i][j] = int(a[i-j])
+ e = zeros(n, dtype=int)
+ for i in range(n):
+ e[i] = int(DB())
+ t = calculate_t(A, s, e, q)
+ return A, t, s, e
+
+def bkz_with_blocksize(instance):
+ A = instance[0]
+ t = instance[1]
+ s = instance[2]
+ e = instance[3]
+ q = instance[4]
+ b = instance[5]
+ B = SVPreduction(A, t, q)
+ params = BKZ.Param(b, strategies=BKZ.DEFAULT_STRATEGY, flags=BKZ.VERBOSE|BKZ.AUTO_ABORT)
+ bkz = BKZ2(GSO.Mat(B, float_type='mpfr'))
+ try:
+ bkz(params)
+ except:
+ return None
+ B_np = makenumpy(B)
+ s_best = B_np[0][:n]
+ e_best = B_np[0][n:-1]
+ sc = score(t, s_best, e_best, A, s, p)
+ if sc == 0:
+ return True
+ else: return False
+
+def time_bkz(task_function, *args, **kwargs):
+ start_time = time.time()
+ result = task_function(*args)
+ end_time = time.time()
+ execution_time = end_time - start_time
+ return result, execution_time
+
+
+# Run LLL
+tries = 128
+pp = [3, 17, 71, 277, 401, 521, 1031, 3329]
+nn = range(8, 81, 8)
+for p in pp:
+ x = zeros(len(nn))
+ k = 0
+ for n in nn:
+ for _ in range(tries):
+ A, t, s, e = nistlike(n, p)
+ B = SVPreduction(A.copy(), t, p)
+ L = LLL.reduction(B)
+ B_np = makenumpy(L)
+ s_best = B_np[0][:n]
+ e_best = B_np[0][n:2*n]
+ sc = score(t, s_best, e_best, A, s, p)
+ if sc == 0:
+ x[k] += 1
+ x[k] = x[k] / tries
+ output = open('lll_results.txt', 'a')
+ output.write(f'{p} {n} {x[k-1]}\n')
+ output.close()
+ k += 1
+
+plot_lll()
+
+
+
+# Run BKZ
+tries = 256
+pp = [71, 401, 3329]
+nn = range(8, 81, 8)
+blocks = [8, 16, 24, 32, 40, 48, 56]
+
+for p in pp:
+ for n in nn:
+ inst = []
+ for _ in range(tries):
+ A, t, s, e = nistlike(n, p)
+ i = (A, t, s, e, p)
+ inst.append(i)
+ for b in blocks:
+ b_inst = [(i[0].copy(), i[1].copy(), i[2].copy(), i[3].copy(), i[4].copy(), b) for i in inst]
+ with multiprocessing.Pool() as pool:
+ result = pool.starmap(time_bkz, [(bkz_with_blocksize, i) for i in b_inst])
+ res = [r for r in result if r[0] != None]
+ if len(res) == 0:
+ # fpylll infinite babei loop error for all instances
+ tt = nan
+ prob = 0
+ else:
+ tt = sum(r[1] for r in res)/len(res)
+ prob = sum(r[0] for r in res)/len(res)
+ output = open('bkz_results.txt', 'a')
+ output.write(f"{n} {p} {b} {prob} {tt} {len(res)}\n")
+ output.close()
+
+plot_bkz()
diff --git a/plots.py b/plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..d47a23f0274ade562a61799955722f80cc1fb6b4
--- /dev/null
+++ b/plots.py
@@ -0,0 +1,156 @@
+import matplotlib.pyplot as plt
+from numpy import log, exp, linspace, inf
+from scipy.optimize import curve_fit
+import scienceplots
+plt.style.use(["science", "grid"])
+
+def sigmoid(x, right, curve):
+ return 1 - 1/(1 + exp(right-curve*x))
+
+col = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22']
+marker = ["v", "p", "h", "8", "s", 'P', 'D', 'X', "*"]
+
+def plot_lll(filepath='results_lll.txt'):
+ d = open(filepath, 'r')
+ dt = d.readlines()
+ n = []
+ q = []
+ x = []
+ for l in dt:
+ ll = l.split()
+ q.append(int(ll[0]))
+ n.append(int(ll[1]))
+ x.append(float(ll[2]))
+ k=0
+ x_val = linspace(1, 80, 1001)
+ unique_q = list(set(q))
+ unique_q.sort()
+ s, q_s = [], []
+ plt.figure(figsize=(4,4))
+ for qs in unique_q:
+ indicies = [i for i in range(len(q)) if q[i] == qs]
+ nn = [n[i] for i in indicies]
+ xx = [x[i] for i in indicies]
+ assert len(nn) == len(xx)
+ popt, pcov = curve_fit(sigmoid, nn,xx, bounds=[(-20, 0), (120, inf)])
+ q_s.append(qs)
+ s.append(popt[0]/popt[1])
+ if qs in [71, 401, 3329]:
+ plt.plot(nn, xx, marker=marker[k%3], linestyle='', markersize=10,
+ color=col[k%3],
+ markerfacecolor=col[k%3],
+ markerfacecoloralt=col[k%3],
+ markeredgecolor=col[k%3], fillstyle='none', label=f"q={qs}")
+
+ yy = sigmoid(x_val, popt[0], popt[1])
+ plt.plot(x_val, yy, '--', color=col[k%3])
+ k += 1
+ plt.xlabel(r'Key length $n$', fontsize=12)
+ plt.ylabel(r'Probability $p$', fontsize=12)
+ legend = plt.legend(fancybox=False, edgecolor="black", loc=3)
+ legend.get_frame().set_linewidth(0.5)
+
+ plt.figure()
+ plt.plot(q_s[1:], s[1:], marker=marker[1], linestyle='', markersize=10, color=col[0], fillstyle='none',
+ markerfacecolor=col[0],
+ markerfacecoloralt=col[0],
+ markeredgecolor=col[0])
+ plt.xlabel(r'Modulus $q$', fontsize=12)
+ plt.ylabel(r'Key length $n$', fontsize=12)
+ plt.show()
+
+
+def plot_bkz(filename='bkz_results.txt'):
+ d = open(filename, 'r')
+ nn, pp, bb, score, time, sample = [], [], [], [], [], []
+ l = d.readline()
+ while l != '':
+ l = l.replace('\n', '')
+ strg = l.split(' ')
+ nn.append(int(strg[0]))
+ pp.append(int(strg[1]))
+ bb.append(int(strg[2]))
+ score.append(float(strg[3]))
+ if strg[4] == 'nan':
+ time.append(-1.0)
+ sample.append(-1)
+ else:
+ time.append(float(strg[4]))
+ sample.append(int(strg[5]))
+ l = d.readline()
+ assert len(nn) == len(pp) == len(bb) == len(score) == len(time) == len(sample)
+ n = list(set(nn))
+ n.sort()
+ p = set(pp)
+ b = list(set(bb))
+ b.sort()
+ results = dict()
+ for block in b:
+ results[block] = dict()
+ for module in p:
+ results[block][module] = dict()
+ for keylength in n:
+ results[block][module][keylength] = dict()
+ idx = [i for i in range(len(nn)) if nn[i] == keylength and pp[i] == module and bb[i] == block]
+ if len(idx) == 1:
+ results[block][module][keylength]['score'] = score[idx[0]]
+ results[block][module][keylength]['time'] = time[idx[0]]
+ results[block][module][keylength]['samplesize'] = sample[idx[0]]
+ elif len(idx) > 1:
+ z = len(idx)
+ samples = [sample[i] for i in idx]
+ times = [time[i] for i in idx]
+ scores = [score[i] for i in idx]
+ weighted_times = sum([times[i]*samples[i] for i in range(z)])/sum(samples)
+ weighted_scores = sum([scores[i]*samples[i] for i in range(z)])/sum(samples)
+ try:
+ results[block][module][keylength]['score'] = weighted_scores
+ results[block][module][keylength]['time'] = weighted_times
+ results[block][module][keylength]['samplesize'] = sum(samples)
+ except:
+ print(block, module, keylength)
+ block_plot = list(set(bb))
+ block_plot.sort()
+ primes = list(set(pp))
+ primes.sort()
+ k = 0
+ x_val = linspace(1, 80, 1001)
+ for prime in primes:
+ plt.figure(figsize=(4,4))
+ for block in block_plot:
+ nnn = []
+ scores = []
+ for keylength in n:
+ nnn.append(keylength)
+ scores.append(results[block][prime][keylength]['score'])
+ plt.plot(nnn, scores, marker=marker[k], linestyle='-', markersize=5,
+ color=col[k],
+ markerfacecolor=col[k],
+ markerfacecoloralt=col[k],
+ markeredgecolor=col[k], fillstyle='none', label=fr"$\beta={block}$")
+ k += 1
+ k %= 9
+ legend = plt.legend(fancybox=False, edgecolor="black", loc=3)
+ legend.get_frame().set_linewidth(0.5)
+ plt.xlabel(r'Key length $n$', fontsize=12)
+ plt.ylabel(r'Probability $p$', fontsize=12)
+ plt.title(f'Prime: {prime}')
+ plt.figure(figsize=(6,4))
+ k = 0
+ for keylength in n:
+ if keylength > 40:
+ b = []
+ times = []
+ for blocks in block_plot:
+ b.append(blocks)
+ times.append(results[blocks][401][keylength]['time'])
+ plt.plot(b, times, marker=marker[k], linestyle='-', markersize=5, color=col[k], fillstyle='none', label=fr"$n$ = {keylength}")
+
+ k += 1
+ legend = plt.legend(fancybox=False, edgecolor="black", bbox_to_anchor=(1.05, 1.0), loc='upper left')
+ legend.get_frame().set_linewidth(0.5)
+ plt.tight_layout()
+ plt.xlabel(r'Block size $\beta$', fontsize=12)
+ plt.ylabel(r'Time $t\: [s]$', fontsize=12)
+ plt.yscale('log')
+ plt.show()
diff --git a/tools.py b/tools.py
new file mode 100644
index 0000000000000000000000000000000000000000..c2aa883f651aa7d96d3f87ac5415e72ed7d2e074
--- /dev/null
+++ b/tools.py
@@ -0,0 +1,109 @@
+from fpylll import LLL, BKZ, IntegerMatrix, GSO
+from fpylll.algorithms.bkz2 import BKZReduction as BKZ2
+from numpy import zeros, array, matrix, dot, arange, nan
+from numpy.linalg import norm
+import random
+import multiprocessing
+import time
+
+class CenteredBinomialDistibution:
+
+ def __init__(self, n):
+ self.n = n
+ def __call__(self):
+ x, y = 0, 0
+ for k in range(self.n):
+ x += random.randint(0, 1)
+ y += random.randint(0, 1)
+ return x - y
+
+def SVPreduction(A, t, q):
+ """
+ Given A, t, q, return the SVP-reduced basis using an embedding technique
+ B = [I_n | -A^T | 0
+ 0 | qI_n | 0
+ 0 | t | 1]
+
+ args:
+ A (matrix): matrix
+ t (array): target vector
+ q (int): modulus
+
+ returns:
+ B (matrix): SVP-reduced basis
+ """
+ d, n = A.shape
+ B = IntegerMatrix(n+d+1,n+d+1)
+ for i in range(n):
+ for j in range(n, n+d):
+ B[i,j] = -A[j - n][i]
+
+ for i in range(n):
+ B[i,i] = 1
+
+ for i in range(n,n+d):
+ B[i,i] = q
+
+ B[-1,-1] = 1
+
+ for i in range(n, n+d):
+ B[-1,i] = t[i-n]
+
+ return B
+
+def calculate_t(A:matrix, s:array, e:array, q:int) -> array:
+ """
+ Given A, s, e, q, calculate t = (A*s + e) mod q
+
+ args:
+ A (matrix): matrix
+ s (array): secret vector
+ e (array): error vector
+ q (int): modulus
+ returns:
+ array t
+ """
+ t = array([(dot(A[j], s) + e[j]) % q for j in range(len(A))])
+ for i in range(len(t)):
+ t[i] = t[i] if t[i] <= q//2 else t[i]-q
+ return t
+
+def makenumpy(B):
+ """
+ Convert fpylll integer matrix to numpy array
+
+ args:
+ B (fpylll integer matrix): matrix
+
+ returns:
+ A (numpy array): matrix
+ """
+ A = zeros((B.nrows, B.ncols), dtype=int)
+ for i in range(B.nrows):
+ for j in range(B.ncols):
+ A[i][j] = int(B[i,j])
+ return A
+
+def score(t, ss, ee, A, s, p):
+ """
+ Given t, s, e, A, p, calculate the score
+ Calculate t' = (A*s + e) mod p
+ Returns ||t-t'||_1
+
+ args:
+ t (array): target vector
+ ss (array): secret vector
+ ee (array): error vector
+ A (matrix): matrix
+ p (int): modulus
+
+ returns:
+ score (int): score
+ """
+ tt = calculate_t(A, ss, ee, p)
+ for k in range(len(t)):
+ if t[k] >= p//2:
+ t[k] = t[k] - p
+ if tt[k] >= p//2:
+ tt[k] = tt[k] - p
+ return int(norm(t-tt, ord=1))