Source code for src.formulas

import math
from numpy import multiply

try:
    import numpy as np
except ImportError:
    print("Warning: Failed to import numpy.")
    np = None

try:
    from scipy.special import lambertw
except ImportError:
    print("Warning: Failed to import scipy.")
    lambertw = None

# Define global constants
PI = 3.14159265358979
e = 2.71828182845905
const = 2 * PI * e
ln2 = 0.693147180559945

# Auxiliary functions

# Low beta Gaussian Heuristic constant for use in NTRU Dense sublattice estimation.
gh_constant = {1: 0.00000, 2: -0.50511, 3: -0.46488, 4: -0.39100, 5: -0.29759, 6: -0.24880, 7: -0.21970, 8: -0.15748,
               9: -0.14673, 10: -0.07541, 11: -0.04870, 12: -0.01045, 13: 0.02298, 14: 0.04212, 15: 0.07014,
               16: 0.09205, 17: 0.12004, 18: 0.14988, 19: 0.17351, 20: 0.18659, 21: 0.20971, 22: 0.22728, 23: 0.24951,
               24: 0.26313, 25: 0.27662, 26: 0.29430, 27: 0.31399, 28: 0.32494, 29: 0.34796, 30: 0.36118, 31: 0.37531,
               32: 0.39056, 33: 0.39958, 34: 0.41473, 35: 0.42560, 36: 0.44222, 37: 0.45396, 38: 0.46275, 39: 0.47550,
               40: 0.48889, 41: 0.50009, 42: 0.51312, 43: 0.52463, 44: 0.52903, 45: 0.53930, 46: 0.55289, 47: 0.56343,
               48: 0.57204, 49: 0.58184, 50: 0.58852}

# Low beta \alpha_\beta quantity as defined in [AC:DucWoe21] for use in NTRU Dense subblattice estimation.
small_slope_t8 = {2: 0.04473, 3: 0.04472, 4: 0.04402, 5: 0.04407, 6: 0.04334, 7: 0.04326, 8: 0.04218, 9: 0.04237,
                  10: 0.04144, 11: 0.04054, 12: 0.03961, 13: 0.03862, 14: 0.03745, 15: 0.03673, 16: 0.03585,
                  17: 0.03477, 18: 0.03378, 19: 0.03298, 20: 0.03222, 21: 0.03155, 22: 0.03088, 23: 0.03029,
                  24: 0.02999, 25: 0.02954, 26: 0.02922, 27: 0.02891, 28: 0.02878, 29: 0.02850, 30: 0.02827,
                  31: 0.02801, 32: 0.02786, 33: 0.02761, 34: 0.02768, 35: 0.02744, 36: 0.02728, 37: 0.02713,
                  38: 0.02689, 39: 0.02678, 40: 0.02671, 41: 0.02647, 42: 0.02634, 43: 0.02614, 44: 0.02595,
                  45: 0.02583, 46: 0.02559, 47: 0.02534, 48: 0.02514, 49: 0.02506, 50: 0.02493, 51: 0.02475,
                  52: 0.02454, 53: 0.02441, 54: 0.02427, 55: 0.02407, 56: 0.02393, 57: 0.02371, 58: 0.02366,
                  59: 0.02341, 60: 0.02332}


[docs] def ball_log_vol(n): """ Calculate the logarithm of the volume of an n-dimensional ball. :param n: Dimension of the ball. :return: Logarithm of the volume of the ball. """ return ((n/2.) * math.log(PI) - math.lgamma(n/2. + 1))
[docs] def log_gh(d, logvol=0): """ Calculate the Gaussian Heuristic constant for a given dimension. :param d: Dimension. :param logvol: Logarithm of the volume. :return: Gaussian Heuristic constant. """ if d < 49: return (gh_constant[d] + logvol/d) return (1./d * (logvol - ball_log_vol(d)))
[docs] def delta(k): """ Calculate the root-Hermite factor for a given block size. :param k: Block size. :return: Root-Hermite factor. """ assert k >= 60 delta = math.exp(log_gh(k)/(k-1)) return (delta)
[docs] def delta_exact(beta): """ Borrowed from the Lattice Estimator :param beta: :return: """ small = ( (2, 1.02190), (5, 1.01862), (10, 1.01616), (15, 1.01485), (20, 1.01420), (25, 1.01342), (28, 1.01331), (40, 1.01295), ) if beta <= 2: return 1.0219 elif beta < 40: for i in range(1, len(small)): if small[i][0] > beta: return RR(small[i - 1][1]) elif beta == 40: return small[-1][1] else: return (beta / (2 * PI * e) * (PI * beta) ** (1 / beta)) ** (1 / (2 * (beta - 1)))
[docs] def check_overstreched(params): """ Check if the parameters are in the overstretched regime. :param params: Dictionary of parameters. :return: Beta value if overstretched, -1 otherwise. """ n = params['lwe_d'] lgq = params['logq'] stds = params['secret_dist'].stddev stde = params['error_dist'].stddev lnq = multiply(lgq, ln2) # dense_det_log = math.log(math.sqrt(stds**2*n)+math.sqrt(stde**2*n)) dense_det_log = math.log(math.sqrt(stds*n)+math.sqrt(stde*n)) for beta in range(50, 1000): for lgq_value in lgq: alpha_beta = delta_exact(beta)**2 alpha_beta_log = math.log(alpha_beta) # print(0.5+lgq_value/(2*alpha_beta_log)) # THIS FAILs. using python in-built round m = round(0.5+lgq_value/(2*alpha_beta_log)) rhs_log = 0.5*(m-1)*lgq_value-0.5*(m-1)**2*alpha_beta_log if dense_det_log < rhs_log: return beta return -1
[docs] def predicted_beta_bdd(n, q, sigma, zeta): """ Predict the beta value for the BDD model. :param n: Dimension. :param q: Modulus. :param sigma: Standard deviation of the error. :param zeta: Ratio of standard deviations. :return: Predicted beta value. """ lnq = math.log(q) beta_approx1 = 1 if n/(lnq) >= 1: beta_approx1 = 2*n/lnq*(math.log(n/(lnq))) # very rough approximation # we are assuming that n > lnq # print("beta approx", beta_approx1, "log(n/lnq)", # math.log(n/(lnq)), "n/lnq", n/lnq, "n", n, "lnq", lnq) A = 2*n*lnq B = 2*(lnq - math.log(sigma*math.sqrt(const)))+math.log(beta_approx1/const) C = n/(2*lnq) D = lnq - math.log(zeta) # approximates beta/(ln(beta/2*const)) Z = ((2*D*math.sqrt(C)+math.sqrt(A))/B)**2 ln_beta_const = -lambertw(-const/Z, -1) # approximates ln(beta/const) num = 2*n*lnq*ln_beta_const denom = ln_beta_const+2*(lnq-math.log(sigma*math.sqrt(const))) - \ 2*math.sqrt(n/(2*lnq*Z))*(lnq-math.log(zeta)) res = num/denom**2 return res
[docs] def predicted_beta_bdd_rev1(n, q, sigma, zeta): lnq = math.log(q) beta_approx1 = const+1 if n/(lnq) >= 1: beta_approx1 = max(const+1, 2*n/lnq*(math.log(n/(lnq)))) # very rough approximation A = 2*n*(lnq-math.log(zeta)) B = (math.log(beta_approx1)+math.log(math.log(beta_approx1/const)))/(0.292*2*math.log(2)) C = math.log(beta_approx1/const) + 2*math.log(q) - 2*math.log(sigma) - math.log(const) # approximates beta/(ln(beta/const)) Z = ((math.sqrt(A)+math.sqrt(A-C*(math.log(A)/(0.292*2*math.log(2)) - B )))/C)**2 ln_beta_const = -lambertw(-const/Z, -1) # approximates ln(beta/const) d_opt = math.sqrt(2*n*(lnq-math.log(zeta))*Z) num = (d_opt - (math.log2(d_opt)-math.log2(beta_approx1))/0.292)*ln_beta_const denom = ln_beta_const+2*(lnq-math.log(sigma*math.sqrt(const))) - \ 2*n/d_opt*(lnq-math.log(zeta)) res = num/denom return res
[docs] def predicted_beta_usvp(d, lnq, sig, chi): """ Predict the beta value for the USVP model. :param d: Dimension. :param lnq: Logarithm of the modulus. :param sig: Standard deviation of the error. :param chi: Ratio of standard deviations. :return: Predicted beta value. """ x = np.divide(d, lnq) # print("d", d, "lnq", lnq) f1 = np.divide(np.multiply(d, np.log(x)), np.multiply(const, lnq - np.log(sig))) # f2 = np.multiply(x, np.log(np.divide(d, lnq - np.log(sig)))) num = d den = lnq - np.log(sig) div = np.divide(num, den) log_div = np.log(div) f2 = np.multiply(x, log_div) beta = np.divide(np.multiply(2 * d, np.multiply(lnq - np.log(chi), np.log(f1))), np.power(lnq + 0.5 * np.log(f2) - np.log(const * sig), 2)) return beta
# Main functions # Eq. (15)
[docs] def model_lambda_usvp(d, logq, std_s, std_e, params): """ Model the lambda value for the USVP model. :param d: Dimension. :param logq: Logarithm of the modulus. :param std_s: Standard deviation of the secret. :param std_e: Standard deviation of the error. :param params: Model parameters. :return: Lambda value. """ sig = std_e chi = std_e/std_s lnq = np.multiply(logq, ln2) beta = predicted_beta_usvp(d, lnq, sig, chi) if isinstance(beta, list): for i, e in enumerate(beta): if e < const: # print("beta", beta[i], "const", const) beta[i] = const + 1 m2 = np.multiply(2 * d, beta * np.divide(lnq - np.log(chi), np.log(np.divide(beta, const)))) return np.multiply(params[0], beta) + np.multiply(params[1], np.log(m2)) + params[2]
# Eq. (17)
[docs] def model_lambda_usvp_s(d, logq, params): """ Model the lambda value for the simplified USVP model. :param d: Dimension. :param logq: Logarithm of the modulus. :param params: Model parameters. :return: Lambda value. """ lnq = np.multiply(logq, ln2) return np.multiply(np.multiply(params[0], np.divide(d, lnq)), np.log(np.divide(params[1] * d, lnq))) + np.multiply(params[2], np.log(d)) + params[3]
[docs] def model_lambda_bdd(d, logq, std_s, std_e, params): """ Model the lambda value for the BDD model. :param d: Dimension. :param logq: Logarithm of the modulus. :param std_s: Standard deviation of the secret. :param std_e: Standard deviation of the error. :param std_s_num: Numerical standard deviation of the secret. :param params: Model parameters. :return: Lambda value. """ chi = std_e/std_s lnq = np.multiply(logq, ln2) beta = [] if isinstance(logq, list): for lq in logq: beta.append(predicted_beta_bdd_rev1(d, 2 ** lq, std_e, std_e/std_s)) #beta.append(predicted_beta_bdd(d, 2 ** lq, std_e, std_e/std_s)) else: beta.append(predicted_beta_bdd_rev1(d, 2 ** logq, std_e, std_e/std_s)) #beta.append(predicted_beta_bdd(d, 2 ** logq, std_e, std_e/std_s)) # Intermediate calculations log_delta = np.log(np.divide(beta,const)) / (np.multiply(2, beta)) m2 = d * np.divide((lnq - np.log(chi)), log_delta) return np.real(np.multiply(params[0], beta) + params[1] * np.log(m2) + params[2])
# Eq. (21)
[docs] def model_lambda_bdd_s(d, logq, params): """ Model the lambda value for the simplified BDD model. :param d: Dimension. :param logq: Logarithm of the modulus. :param params: Model parameters. :return: Lambda value. """ lnq = np.multiply(logq, ln2) return np.multiply(np.divide(params[0] * d, lnq), np.log(params[1] * d / lnq)) + np.multiply(params[2], np.log(d)) + params[3]
# Eq. (22)
[docs] def model_n_usvp(l, logq, std_s, std_e, params): """ Model the n value for the USVP model. :param l: Security parameter. :param logq: Logarithm of the modulus. :param std_s: Standard deviation of the secret. :param std_e: Standard deviation of the error. :param params: Model parameters. :return: n value. """ beta_approx = l/0.292 sigma = std_e chi = std_e/std_s lnq = np.multiply(logq, ln2) num = (0.5*np.log(beta_approx)+lnq-np.log(const*sigma)+params[1])**2 denom = 2*(np.log(beta_approx/const + params[2])*(lnq-np.log(chi))) leading_order = (beta_approx*params[0])*num / denom return leading_order
# Eq. (23)
[docs] def model_n_usvp_s(l, logq, params): """ Model the n value for the simplified USVP model. :param l: Security parameter. :param logq: Logarithm of the modulus. :param params: Model parameters. :return: n value. """ lnq = np.multiply(logq, ln2) return np.multiply(np.divide(l + params[0] * np.log(lnq), params[1] * np.log(l) + params[2]) + params[3], lnq)
# Eq. (24)
[docs] def model_n_bdd_rev1(l, logq, std_s, std_e, params): """ Model the n value for the BDD model. :param l: Security parameter. :param logq: Logarithm of the modulus. :param std_s: Standard deviation of the secret. :param std_e: Standard deviation of the error. :param params: Model parameters. :return: n value. """ sigma = std_e zeta = std_e / std_s beta_approx = (l - np.log(l))/0.292 # approximate beta from lambda lnq = np.multiply(logq, ln2) n = l A = 2*lnq B = beta_approx C = np.log(beta_approx/const) D = lnq - np.log(zeta) denom = 9*(C + params[3])*D nom = (params[0] * B + params[1])*(A + C + params[2])**2 return nom/denom
[docs] def model_n_bdd_rev1_exact(l, logq, std_s, std_e, params): zeta = std_e / std_s beta_approx = (l - np.log(8*l)-16.4)/0.292+params[0] # approximate beta from lambda lnq = np.multiply(logq, ln2) logn_approx = np.log2(0.5*beta_approx*(lnq)/np.log(beta_approx/const)) A = 2*(lnq-np.log(zeta))*beta_approx/(np.log(beta_approx/const)) #A = (lnq)*beta_approx/(math.log(beta_approx/const)) B = np.log(beta_approx/const) C = np.log(beta_approx/const)+2*lnq - 2*np.log(std_e)-np.log(const) D = (lnq-np.log(zeta))*(np.log(beta_approx/const))/(2*beta_approx) #D = (lnq)*(math.log(beta_approx/const))/beta_approx E = 0.5*logn_approx+0.5*np.log2(A)-np.log2(beta_approx) nom = E*B + beta_approx*C denom = 2*beta_approx*np.sqrt(D)+B*np.sqrt(A) # nom = E*B+beta_approx*C # denom = 2*beta_approx*math.sqrt(D)+B*math.sqrt(A) return (params[2]*nom/denom+params[1])**2
# Eq. (25)
[docs] def model_n_bdd_s(l, logq, params): """ Model the n value for the simplified BDD model. :param l: Security parameter. :param logq: Logarithm of the modulus. :param std_s: Standard deviation of the secret. :param std_e: Standard deviation of the error. :param params: Model parameters. :return: n value. """ lnq = np.multiply(logq, ln2) term1 = params[0] * l + params[1] * np.log(l) term2 = params[2] * np.divide(lnq, np.log(l)) + \ params[3] * np.divide(np.log(l), lnq) return np.multiply(term1, term2)