From f7f39e4db2f3e22a21e1dd635e0601caae2b4510 Mon Sep 17 00:00:00 2001 From: FreVercauteren Date: Fri, 4 Feb 2022 12:23:51 +0100 Subject: [PATCH] Added scripts --- Scripts/Python_Scripts/run.sh | 7 + .../__pycache__/proba_util.cpython-38.pyc | Bin 0 -> 5380 bytes .../select_params/proba_util.py | 194 ++++++++++++++++ .../select_params/select_params.py | 121 ++++++++++ Scripts/Sage_Scripts/sanity_check.sage | 218 ++++++++++++++++++ 5 files changed, 540 insertions(+) create mode 100644 Scripts/Python_Scripts/run.sh create mode 100644 Scripts/Python_Scripts/select_params/__pycache__/proba_util.cpython-38.pyc create mode 100644 Scripts/Python_Scripts/select_params/proba_util.py create mode 100644 Scripts/Python_Scripts/select_params/select_params.py create mode 100644 Scripts/Sage_Scripts/sanity_check.sage diff --git a/Scripts/Python_Scripts/run.sh b/Scripts/Python_Scripts/run.sh new file mode 100644 index 0000000..552c81d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Python_Scripts/select_params/__pycache__/proba_util.cpython-38.pyc b/Scripts/Python_Scripts/select_params/__pycache__/proba_util.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..7239ed3cb7060951b23f22f638144fc25e54fa03 GIT binary patch literal 5380 zcmcIo&vVdd66>Pc)_zNCpJW4BJL#^pfRrG%0M z$OWih0Z(dgJ$2g7bf)ddlapt9>!F9<+W!LA9x}f5+*7}|ASIHvJ56Q+4t5ud1+e?Q z?|ttre>^o+G;n;evR9j#H;jML#qeR`;x5kU9EC6hYZ@)KYcgt^Xj{8h%igtHx!qhV zzne$DCG1}sy9JRGdDM<52nY3qm=HzOMKLK#s3*mgm_}U^Gh!C?l(-`1P*01iVjlI3 zxF!}*&x%D+MtwyriDlGt;<{KteO25Luc4k7uM3X)nt0=pQMvgXvotC!ov6EB94g;! zR_xSjhWn}QNxzv|Qq@yC>ZsVGl?@JR=WyP~836_(HWHQ?F%#^J%^TN^GlsXDr+?y~ zN`L<#=3$-32a@l*-6-WZKa%LPU_+U1?=^hcWGpaVY8~2ylCvW2wc)f_@1zzo1WA!)# z_c^UY6;M>NMFKb0P&3p^5DejSZ7{XkU~09tsK1Fp1O8k?+avG?2l3|uq+7V~jUI8w zn8YaxdAa^Y+A0#}_3 zljv)Ja&P!u{PW|@3?*I<6+ec2@Xyy8>+6kPPk-lmpJO@Y24E_!ykaTRtdwo1%pV8u z`d-&>ioWYfaS0);OTa|rW=7s-^K6!>Wwbp4Z*YtNxq{wUt72k;Z<=$YjfuH1P{@Dq zk3-RIb|W4g`LTBZ<>0G|_KEnTP<5jb4+Ntipj}DN z>vz4}Q9GMJtKE%xt3QR8y0gNtup7U6Twg~2{Q-G096!zJZvtf^XUaOMsXF5@7xB$n zw&eFPYCuem_K5cqo^;;F8F3UOV+)Y51QZnPy3sIHnIwy@%jCc2>94tx0!*a1g3s{I zEk5uF03mU=tz<;nDu^y?o$77AzRpp711PGin6zGf9}q4Scw9t508TykWHholLMFZN z>zbt4Rs%{O;kk$cm65Yy+milGvoc#S)wj_*VrNQ9=+qvtjRx5{F<{?!dsaerzSXl~ z8+*~{hJ?0+Y$wY=+H%~xYv{s*9N4qnUc^GBTjU0lh=I}#h3T_>QE>c8N0Q~BZ+B{SCZV| zWiJ`5ZT5_*Hvz@7e9G>nY%68=)f>1@*kC?3v)TC0N5&CGa{+G)nzMEe=j2} z1xn9DQmx3@>sStmeSfb5)6k@ z+S}+GVJ6A`1%wZAA($|IWblxFWc#$G{FS+cP>@)FbVG*$1WZW13E?--Ffy@4u3@Sb z1OF;lE|BT8LF*o8o>LEG0{A;wYszcs-h zPnw4$naF9zBaB6}6JsKdi4NI7bMadpSnTY*XFQ3@wEXgOi1=#B;P*3tp zB+LlJ)JNGgo-F9mI~YyUzTI#Dd0^Lmv4cSI3)9fcP(W>ypT-b ztsi!_k^&`0JN^BZqP$e=XQenq&nz!#?-m0@V7KX?N|o9z_t@3aswnwYYgDXLF(w2z zhV-?jWGg(BwL!e_PJV+YB8sqv4tEZ_X$>Tz#5S1vKJL80@C!_LaaUyQHqO~o!jxpe zYUI?4!2TmPV=J~}BsXXHCDs`j4(x6sECS_o;GQ8xtx1jjqVYWh(a!CDDEzb5(BVxV z1RzFU!i92=f-J~wb00}MB*%-JCwyN?D2>XLf}$ID_aX|PGRs|tRe5`pe?*~d`*I63 z;t%mAk6pi+Z5$=Nbu_2R%vMM(7$dbC@+3+fB@d;Fq{tE$ZBXK5X4RY~y_!Hr%A(Y4 zxk|gjn9XMCtEAq8z3x*y5KW_?z*qv|ZMKYkwfX=(BN}d@<3fPGg$tPHQp4vVFM6AZ z`~&+QK+QUR{oeT!gY|)3><`d@@7kn$sfCpMrL%PCb_@FWR}9Z`Vv}Y!#}s+o$h>Gy z7dbIGYL*N4uv4I%iO2(Dtbjh`ev}`9NM-ypcNNZ~DKW1pVeP29cv|P~_;ZN9lg-vv zl$paVFdwx=-7f>SL{ru|QsBwIV^pT-CIMPv%cgpOzMaYfxmSFk^v<6C-cf|B(sra( z=>Mbie?9tN4ZU~OJI>1lRi}%Y@KK=E>ZcV-yOcIuHhcQwkk4u+WPhgI(TTKtpV$&{ c6eU)+oC#+NX9?$|bH$m))uc0n-<;$87wFLc3IG5A literal 0 HcmV?d00001 diff --git a/Scripts/Python_Scripts/select_params/proba_util.py b/Scripts/Python_Scripts/select_params/proba_util.py new file mode 100644 index 0000000..53f9cbc --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/Python_Scripts/select_params/select_params.py b/Scripts/Python_Scripts/select_params/select_params.py new file mode 100644 index 0000000..855a20c --- /dev/null +++ b/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() diff --git a/Scripts/Sage_Scripts/sanity_check.sage b/Scripts/Sage_Scripts/sanity_check.sage new file mode 100644 index 0000000..01b39b0 --- /dev/null +++ b/Scripts/Sage_Scripts/sanity_check.sage @@ -0,0 +1,218 @@ +#!/usr/bin/env sage + +n = int(sys.argv[1]) +q = int(sys.argv[2]) +numinputs = int(sys.argv[3]) + +set_random_seed(1) # reproducibility +proof.all(False) + + +try: + os_cores = len(os.sched_getaffinity(0)) +except: + os_cores = cpu_count() +os_cores = os.getenv('CORES',default=os_cores) +os_cores = int(os_cores) +if os_cores < 1: os_cores = 1 + + +import fpylll + +# not just for speed; also deals with "cannot open strategies file" fpylll bug +try: + fn = fpylll.BKZ.DEFAULT_STRATEGY + if not fn.startswith('/'): + fn = fpylll.BKZ.DEFAULT_STRATEGY_PATH+'/'+fn +except TypeError: + fn = fpylll.BKZ.DEFAULT_STRATEGY + if not fn.startswith(b'/'): + fn = fpylll.BKZ.DEFAULT_STRATEGY_PATH+b'/'+fn +precomputed_bkz_strategies = fpylll.load_strategies_json(fn) + +def bkz(M,beta): + return M.BKZ(block_size=beta,strategies=precomputed_bkz_strategies) + + +from sage.doctest.util import Timer +timer = Timer() + +Zx. = ZZ[] +R = Zx.quotient(x^n+1) + +def reduce(r): + assert r in R + return R(sum((r[i]%q)*x^i for i in range(n))) + # "r %= q" in Sage sometimes produces wrong reduction mod q + # for example, try this: R. = ZZ[]; (-x) % 3 + +def divq(r): + assert r in R + assert all(r[i]%q == 0 for i in range(n)) + return R(sum(ZZ(r[i]/q)*x^i for i in range(n))) + +def randomsubset(n,m): + L = [2*randrange(2^64)+1 for j in range(m)] + [2*randrange(2^64) for j in range(n-m)] + L.sort() + return tuple(j for j in range(n) if L[j]%2 == 1) + + +@parallel(ncpus=os_cores) +def doit(rounding,secreta1,secreta2,secrete1,secrete2,secretr1,secretr2,G11,G12,G21,G22,A1,A2): + result = [] + + for r in [G11,G12,G21,G22,A1,A2]: + assert all(r[i]>=0 for i in range(n)) + assert all(r[i]= 0 for j in S) + assert all(j < 2*n for j in S) + + def Extract(e1,e2): + assert e1 in R + assert e2 in R + e1e2 = list(e1)+list(e2) + return [e1e2[j] for j in S] + + def P(z,a1,a2,r1,r2): + assert z in ZZ + assert a1 in R + assert a2 in R + assert r1 in R + assert r2 in R + e1 = z*A1-a1*G11-a2*G21+r1*q + e2 = z*A2-a1*G12-a2*G22+r2*q + return Extract(e1,e2)+[z]+list(a1)+list(a2) + + secretvector = Extract(secrete1,secrete2)+[1]+list(secreta1)+list(secreta2) + assert secretvector == P(1,secreta1,secreta2,secretr1,secretr2) + + Cinbasis = [] + Cinbasis += [(1,R(0),R(0),R(0),R(0))] + for i in range(n): + Cinbasis += [(0,R(x^i),R(0),R(0),R(0))] + for i in range(n): + Cinbasis += [(0,R(0),R(x^i),R(0),R(0))] + for i in range(len(S)): + if S[i] < n: + Cinbasis += [(0,R(0),R(0),R(x^S[i]),R(0))] + else: + Cinbasis += [(0,R(0),R(0),R(0),R(x^(S[i]-n)))] + + Coutbasis = [P(z,a1,a2,r1,r2) for z,a1,a2,r1,r2 in Cinbasis] + + assert d == len(Cinbasis) + assert d == len(Coutbasis) + + M = matrix(Coutbasis) + assert abs(det(M)) == q^m + + success = False + time = 0 + + for beta in range(2,21,2): + if not success: + timer.start() + if beta == 2: + M = M.LLL() + else: + M = bkz(M,beta) + time += timer.stop().cputime + + # fpylll BKZ doesn't always put shortest basis vector first + v = M[0] + for i in range(1,M.nrows()): + if M[i]*M[i] < v*v: + v = M[i] + + if v[m] < 0: v = -v + + a1 = R(list(v[m+1:m+1+n])) + a2 = R(list(v[m+1+n:m+1+2*n])) + + if (a1,a2) == (secreta1,secreta2): + success = True + + result += [(rounding,m,d,beta,success,time)] + + return result + + +def randomsmallcoeff(): + return sum(randrange(2) for coin in range(4)) - 2 + +def randomsmall(): + return R([randomsmallcoeff() for i in range(n)]) + +def randomq(): + return R([randrange(q) for i in range(n)]) + +def round8offset(r): + result = R([((4-r[i])%8)-4 for i in range(n)]) + assert all(result[i] >= -4 for i in range(n)) + assert all(result[i] <= 3 for i in range(n)) + assert all((r[i]+result[i])%8 == 0 for i in range(n)) + return result + +def noise8(): + return R([randrange(-4,4) for i in range(n)]) + +def inputs(): + for i in range(numinputs): + secreta1 = randomsmall() + secreta2 = randomsmall() + + G11 = randomq() + G12 = randomq() + G21 = randomq() + G22 = randomq() + + for rounding in False,True: + Z1 = secreta1*G11+secreta2*G21 + Z2 = secreta1*G12+secreta2*G22 + + if rounding: + secrete1 = round8offset(Z1) + secrete2 = round8offset(Z2) + else: + secrete1 = noise8() + secrete2 = noise8() + + Z1 += secrete1 + Z2 += secrete2 + + A1 = reduce(Z1) + A2 = reduce(Z2) + secretr1 = divq(Z1-A1) + secretr2 = divq(Z2-A2) + + yield rounding,secreta1,secreta2,secrete1,secrete2,secretr1,secretr2,G11,G12,G21,G22,A1,A2 + +def strpoly(r): + return str(list(r)).replace(' ','') + +stats = {} + +for control,result in sorted(list(doit(inputs()))): + rounding,secreta1,secreta2,secrete1,secrete2,secretr1,secretr2,G11,G12,G21,G22,A1,A2 = control[0] + for rounding,m,d,beta,success,time in result: + sys.stdout.write( + 'n %d q %d rounding %s m %d d %d beta %d success %s time %.6f a1 %s a2 %s e1 %s e2 %s G11 %s G12 %s G21 %s G22 %s A1 %s A2 %s\n' + % (n,q,rounding,m,d,beta,success,time,strpoly(secreta1),strpoly(secreta2),strpoly(secrete1),strpoly(secrete2),strpoly(G11),strpoly(G12),strpoly(G21),strpoly(G22),strpoly(A1),strpoly(A2)) + ) + + if (m,d,beta,rounding) not in stats: + stats[m,d,beta,rounding] = [] + stats[m,d,beta,rounding] += [success] + +for m,d,beta,rounding in sorted(stats): + S = stats[m,d,beta,rounding] + assert len(S) == numinputs + print('n %d q %d rounding %s m %d d %d beta %d percentage %.9f numinputs %d' + % (n,q,rounding,m,d,beta,sum(S)/len(S),len(S)))