Added scripts

Scripts/Python_Scripts/run.sh

@@ -0,0 +1,7 @@
+# clone the estimator and setup in framework folder
+git clone https://github.com/lducas/leaky-LWE-Estimator.git tmp
+cp -r tmp/framework framework
+
+# run the script
+cd select_params
+sage select_params.py

Scripts/Python_Scripts/select_params/proba_util.py

@@ -0,0 +1,194 @@
+#/*---------------------------------------------------------------------
+#This file has been adapted from the implementation
+#(available at, Public Domain https://github.com/pq-crystals/kyber)
+#of "CRYSTALS - Kyber: a CCA-secure module-lattice-based KEM"
+#by : Joppe Bos, Leo Ducas, Eike Kiltz, Tancrede Lepoint,
+#Vadim Lyubashevsky, John M. Schanck, Peter Schwabe & Damien stehle
+#----------------------------------------------------------------------*/ 
+
+from math import factorial as fac
+from math import log, ceil, erf, sqrt
+
+
+def gaussian_center_weight(sigma, t):
+    """ Weight of the gaussian of std deviation s, on the interval [-t, t]
+    :param x: (float)
+    :param y: (float)
+    :returns: erf( t / (sigma*\sqrt 2) )
+    """
+    return erf(t / (sigma * sqrt(2.)))
+
+
+def binomial(x, y):
+    """ Binomial coefficient
+    :param x: (integer)
+    :param y: (integer)
+    :returns: y choose x
+    """    
+    try:
+        binom = fac(x) // fac(y) // fac(x - y)
+    except ValueError:
+        binom = 0
+    return binom
+
+
+def centered_binomial_pdf(k, x):
+    """ Probability density function of the centered binomial law of param k at x
+    :param k: (integer)
+    :param x: (integer)
+    :returns: p_k(x)
+    """    
+    return binomial(2*k, x+k) / 2.**(2*k)
+
+
+def build_centered_binomial_law(k):
+    """ Construct the binomial law as a dictionnary
+    :param k: (integer)
+    :param x: (integer)
+    :returns: A dictionnary {x:p_k(x) for x in {-k..k}}
+    """    
+    D = {}
+    for i in range(-k, k+1):
+        D[i] = centered_binomial_pdf(k, i)
+    return D
+
+
+def mod_switch(x, q, rq):
+    """ Modulus switching (rounding to a different discretization of the Torus)
+    :param x: value to round (integer)
+    :param q: input modulus (integer)
+    :param rq: output modulus (integer)
+    """    
+    return int(round(1.* rq * x / q) % rq)
+
+
+def mod_centered(x, q):
+    """ reduction mod q, centered (ie represented in -q/2 .. q/2)
+    :param x: value to round (integer)
+    :param q: input modulus (integer)
+    """        
+    a = x % q
+    if a < q/2:
+        return a
+    return a - q
+
+
+def build_mod_switching_error_law(q, rq):
+    """ Construct Error law: law of the difference introduced by switching from and back a uniform value mod q
+    :param q: original modulus (integer)
+    :param rq: intermediate modulus (integer)
+    """        
+    D = {}
+    V = {}
+    for x in range(q):
+        y = mod_switch(x, q, rq)
+        z = mod_switch(y, rq, q)
+        d = mod_centered(x - z, q)
+        D[d] = D.get(d, 0) + 1./q
+        V[y] = V.get(y, 0) + 1
+
+    return D
+
+
+def law_convolution(A, B):
+    """ Construct the convolution of two laws (sum of independent variables from two input laws)
+    :param A: first input law (dictionnary)
+    :param B: second input law (dictionnary)
+    """        
+
+    C = {}
+    for a in A:
+        for b in B:
+            c = a+b
+            C[c] = C.get(c, 0) + A[a] * B[b]
+    return C
+
+
+def law_product(A, B):
+    """ Construct the law of the product of independent variables from two input laws
+    :param A: first input law (dictionnary)
+    :param B: second input law (dictionnary)
+    """        
+    C = {}
+    for a in A:
+        for b in B:
+            c = a*b
+            C[c] = C.get(c, 0) + A[a] * B[b]
+    return C
+
+
+def clean_dist(A):
+    """ Clean a distribution to accelerate further computation (drop element of the support with proba less than 2^-300)
+    :param A: input law (dictionnary)
+    """        
+    B = {}
+    for (x, y) in A.items():
+        if y>2**(-300):
+            B[x] = y
+    return B
+
+
+def iter_law_convolution(A, i):
+    """ compute the -ith forld convolution of a distribution (using double-and-add)
+    :param A: first input law (dictionnary)
+    :param i: (integer)
+    """    
+    D = {0: 1.0}
+    i_bin = bin(i)[2:]  # binary representation of n
+    for ch in i_bin:
+        D = law_convolution(D, D)
+        D = clean_dist(D)
+        if ch == '1':
+            D = law_convolution(D, A)
+            D = clean_dist(D)
+    return D
+
+def convolution_remove_dependency(A, B, q, p):
+    normalizer={} 
+    maxa=int(q/p)
+    for a in A:
+    	normalizer[a%maxa]=normalizer.get(a%maxa,0)+A[a]
+
+
+    C={}
+    for a in A:
+    	for b in B:
+    		c=a-b
+    		if (c%maxa==0):
+    			C[c] = C.get(c, 0) + A[a] * B[b]/normalizer[a%maxa]
+    return C
+
+def tail_probability(D, t):
+    '''
+    Probability that an drawn from D is strictly greater than t in absolute value
+    :param D: Law (Dictionnary)
+    :param t: tail parameter (integer)
+    '''
+    s = 0
+    ma = max(D.keys())
+    if t >= ma:
+        return 0
+    for i in reversed(range(int(ceil(t)), ma)):  # Summing in reverse for better numerical precision (assuming tails are decreasing)
+        s += D.get(i, 0) + D.get(-i, 0)
+    return s
+
+
+
+## distribution analytics
+
+def distmean(dist):
+    # get mean of a distribution dist
+    res = 0.
+    for i in dist.keys():
+        res += dist[i] * i
+    return res
+
+
+def distvariance(dist):
+    # get variance of a distribution dist
+    mean = distmean(dist)
+    res = 0.
+    for i in dist.keys():
+        res += dist[i] * (i - mean)**2
+
+    return res

Scripts/Python_Scripts/select_params/select_params.py

@@ -0,0 +1,121 @@
+#/*---------------------------------------------------------------------
+#This file has been adapted from the implementation
+#(available at, Public Domain https://github.com/pq-crystals/kyber)
+#of "CRYSTALS - Kyber: a CCA-secure module-lattice-based KEM"
+#by : Joppe Bos, Leo Ducas, Eike Kiltz, Tancrede Lepoint,
+#Vadim Lyubashevsky, John M. Schanck, Peter Schwabe & Damien stehle
+#----------------------------------------------------------------------*/ 
+
+import numpy as np
+import matplotlib.pyplot as plt
+from math import sqrt, exp, log, floor
+from proba_util import *
+from sage.all import *
+
+load("../framework/instance_gen.sage")
+
+# core SVP cost models
+cost_model_c = lambda beta: 0.292 * beta
+cost_model_q = lambda beta: 0.265 * beta
+
+# estimate cost using leaky-LWE-Estimator https://eprint.iacr.org/2020/292.pdf
+def estimate(n, q, m, D_e, D_s):
+    A, b, dbdd = initialize_from_LWE_instance(DBDD_predict, n, q, n, D_e, D_s, verbosity=0)
+    _ = dbdd.integrate_q_vectors(q, report_every=100)
+    beta, delta = dbdd.estimate_attack()
+    return beta
+
+
+# search through parameter sets given D_s, q, k, n
+def search_params(D_s, q, k, n, minsec, maxerror):
+    logq=int(log(q,2))
+
+    # loop over all p values that are a power of 2 (from 10, others are too low, by experiment)
+    for logp in range(10,int(floor(log(q,2)))):
+        p=int(2**logp)
+        D_e = build_mod_switching_error_law(q, p)
+
+        # security estimate
+        beta = estimate(k*n, q, k*n, D_e, D_s)
+        sec_c = cost_model_c(beta)
+
+        # stop and break if under minimum security
+        if sec_c < minsec:
+            break
+
+        # failure calculation part 1
+        se = law_product(D_s, D_e)
+        se2 = iter_law_convolution(se, k*n)
+        se2 = convolution_remove_dependency(se2, se2, q, p)
+
+        ### loop over all reconciliation values (note that p - T < q - p so that the security proof works)
+        for logT in range(1,2*logp-logq+1):
+            T=2**logT
+
+            # failure calculation part 2
+            e2 = build_mod_switching_error_law(q, T)
+            D = law_convolution(se2, e2)
+
+            prob = tail_probability(D, q/4.)
+            if prob!=0:
+    	        prob = log(256*prob,2)
+
+            # if too low, search for a bigger T
+            if prob > maxerror:
+                continue
+
+            ######################################################
+            ### We have found good parameters, give a printout ###
+
+            ### bandwidth calculation
+            # size of b in bytes
+            BWb=logp*k*n/8
+            # size of c in bytes
+            BWc = n*logT/8
+
+            BWsfull = logq*k*n/8
+
+            print('-- parameters --')
+            print('q: ',log(q,2),'p: ',logp,'T: ',logT, 'l: ', k, 'n: ', n, 'D_s: ', D_s)
+
+            print('-- bandwidth --')
+            print('bandwidth b', BWb, 'bandwidth c', BWc, 'bandwidth total', 2*BWb + BWc + 256/8)
+            print('KEM')
+            #pk: b, seedA; sk: b, s, (z,seedA,H(pk))
+            print('pk: ',BWb + 32, 'send: ', BWb+BWc)
+            print('PKE')
+            print('pk: ',BWb + 32, 'send: ', BWb+BWc)
+
+            print('-- correctness --')
+            print('failure probability: ', float(prob))
+
+            print('-- security --')
+            print('quantum: ', cost_model_q(beta), 'classical: ', cost_model_c(beta))
+
+            print('\n\n')
+
+            # succes, we found a b that is secure, and don't need to increase b 
+            break
+
+
+def main():
+    # search_params(binomial coins, q, k, n, log2(classical_security), log2(failure probability)
+    # script looks for optimal p and t
+    # parameter for binomial distribution in Kyber is (binomial coins)/2
+    print('-lightsaber-')
+    D_s = build_centered_binomial_law(5)
+    search_params(D_s, q = 2**13, k = 2, n = 256, minsec = 113, maxerror = -100)
+
+    print('-saber-')
+    D_s = build_centered_binomial_law(4)
+    search_params(D_s, q = 2**13, k = 3, n = 256, minsec = 177, maxerror = -128)
+
+    print('-firesaber-')
+    D_s = build_centered_binomial_law(3)
+    search_params(D_s, q = 2**13, k = 4, n = 256, minsec = 241, maxerror = -160)
+
+
+
+
+if __name__ == '__main__':
+    main()