Merge pull request #127 from ludopulles/change-ND-fix-ST

Thoroughly change NoiseDistribution

docs/algorithms/lwe-bkw.rst

@@ -6,7 +6,7 @@ Coded-BKW for LWE
 We construct an example LWE instance::
 
     from estimator import *
-    params = LWE.Parameters(n=400, q=7981, Xs=ND.SparseTernary(384, 16), Xe=ND.CenteredBinomial(4), m=800)
+    params = LWE.Parameters(n=400, q=7981, Xs=ND.SparseTernary(16), Xe=ND.CenteredBinomial(4), m=800)
     params
 
 and estimate the cost of Coded-BKW [C:GuoJohSta15]_, [C:KirFou15]_::

docs/algorithms/lwe-dual.rst

@@ -7,14 +7,14 @@ We construct an (easy) example LWE instance::
 
     from estimator import *
     from estimator.lwe_dual import dual_hybrid, matzov
-    params = LWE.Parameters(n=200, q=7981, Xs=ND.SparseTernary(384, 16), Xe=ND.CenteredBinomial(4))
+    params = LWE.Parameters(n=200, q=7981, Xs=ND.SparseTernary(16), Xe=ND.CenteredBinomial(4))
     params
 
-The simples (and quickest to estimate) algorithm is the "plain" dual attack as described in [PQCBook:MicReg09]_::
+The simplest (and quickest to estimate) algorithm is the "plain" dual attack as described in [PQCBook:MicReg09]_::
 
     LWE.dual(params)
 
-We can improve these results by considering a dual hybrid attack as in [EC:Albrecht17,INDOCRYPT:EspJouKha20]_::
+We can improve these results by considering a dual hybrid attack as in [EC:Albrecht17]_, [INDOCRYPT:EspJouKha20]_::
 
     dual_hybrid(params)
 

docs/algorithms/lwe-primal.rst

@@ -6,7 +6,7 @@ LWE Primal Attacks
 We construct an (easy) example LWE instance::
 
     from estimator import *
-    params = LWE.Parameters(n=200, q=7981, Xs=ND.SparseTernary(384, 16), Xe=ND.CenteredBinomial(4))
+    params = LWE.Parameters(n=200, q=7981, Xs=ND.SparseTernary(16), Xe=ND.CenteredBinomial(4))
     params
 
 The simplest (and quickest to estimate) model is solving via uSVP and assuming the Geometric Series
@@ -21,9 +21,7 @@ we optimize β and d separately::
 
     LWE.primal_usvp(params, red_shape_model=Simulator.GSA)
 
-To get a more precise answer we may use the CN11 simulator by Chen and Nguyen [AC:CheNgu11]_ (as
-`implemented in FPyLLL
-<https://github.com/fplll/fpylll/blob/master/src/fpylll/tools/bkz_simulator.py>_`)::
+To get a more precise answer we may use the CN11 simulator by Chen and Nguyen [AC:CheNgu11]_ (as `implemented in FPyLLL <https://github.com/fplll/fpylll/blob/master/src/fpylll/tools/bkz_simulator.py>`__)::
 
     LWE.primal_usvp(params, red_shape_model=Simulator.CN11)
 

docs/schemes/hes.rst

@@ -1,5 +1,5 @@
 Homomorphic Encryption Parameters
-===============================
+=================================
 
 ::
 

estimator/__init__.py

@@ -2,11 +2,11 @@
 
 __all__ = ['ND', 'Logging', 'RC', 'Simulator', 'LWE', 'NTRU', 'SIS', 'schemes']
 
-from .nd import NoiseDistribution as ND
 from .io import Logging
 from .reduction import RC
 from . import simulator as Simulator
 from . import lwe as LWE
 from . import ntru as NTRU
+from . import nd as ND
 from . import sis as SIS
 from . import schemes

estimator/gb.py

@@ -213,8 +213,7 @@ def __call__(
             rop: ≈2^227.2, dreg: 54, mem: ≈2^227.2, t: 4, m: 1024, tag: arora-gb
             >>> LWE.arora_gb(params.updated(Xs=ND.UniformMod(3), Xe=ND.CenteredBinomial(4), m=1024))
             rop: ≈2^189.9, dreg: 39, mem: ≈2^189.9, t: 4, m: 1024, tag: arora-gb
-            >>> Xs, Xe =ND.SparseTernary(1024, 64, 0), ND.DiscreteGaussian(2**10)
-            >>> LWE.arora_gb(LWE.Parameters(n=1024, q=2**40, Xs=Xs, Xe=Xe))
+            >>> LWE.arora_gb(LWE.Parameters(n=1024, q=2**40, Xs=ND.SparseBinary(64), Xe=ND.DiscreteGaussian(2**10)))
             rop: ≈2^inf, dreg: ≈2^inf, tag: arora-gb
 
         ..  [EPRINT:ACFP14] Martin R. Albrecht, Carlos Cid, Jean-Charles Faugère & Ludovic Perret. (2014).

estimator/lwe_bkw.py

@@ -239,7 +239,7 @@ def sf(x, best):
         # the search cannot fail. It just outputs some X with X["oracle"]>m.
         if best["m"] > params.m:
             raise InsufficientSamplesError(
-                f"Got m≈2^{float(log(params.m, 2.0)):.1f} samples, but require ≈2^{float(log(best['m'],2.0)):.1f}.",
+                f"Got m≈2^{float(log(params.m, 2.0)):.1f} samples, but require ≈2^{float(log(best['m'], 2.0)):.1f}.",
                 best["m"],
             )
         return best

estimator/lwe_dual.py

@@ -7,7 +7,6 @@
 """
 
 from functools import partial
-from dataclasses import replace
 
 from sage.all import oo, ceil, sqrt, log, cached_function, RR, exp, pi, e, coth, tanh
 
@@ -19,7 +18,7 @@
 from .io import Logging
 from .conf import red_cost_model as red_cost_model_default, mitm_opt as mitm_opt_default
 from .errors import OutOfBoundsError, InsufficientSamplesError
-from .nd import NoiseDistribution
+from .nd import DiscreteGaussian, SparseTernary
 from .lwe_guess import exhaustive_search, mitm, distinguish
 
 
@@ -61,22 +60,26 @@ def dual_reduce(
 
         # Compute new secret distribution
         if params.Xs.is_sparse:
-            h = params.Xs.get_hamming_weight(params.n)
+            h = params.Xs.hamming_weight
             if not 0 <= h1 <= h:
                 raise OutOfBoundsError(f"Splitting weight {h1} must be between 0 and h={h}.")
-            # assuming the non-zero entries are uniform
-            p = h1 / 2
-            red_Xs = NoiseDistribution.SparseTernary(params.n - zeta, h / 2 - p)
-            slv_Xs = NoiseDistribution.SparseTernary(zeta, p)
+
+            if type(params.Xs) is SparseTernary:
+                # split the +1 and -1 entries in a balanced way.
+                slv_Xs, red_Xs = params.Xs.split_balanced(zeta, h1)
+            else:
+                # TODO: Implement this for sparse secret that are not SparseTernary,
+                # i.e. DiscreteGaussian with extremely small stddev.
+                raise NotImplementedError(f"Unknown how to exploit sparsity of {params.Xs}")
 
             if h1 == h:
                 # no reason to do lattice reduction if we assume
                 # that the hw on the reduction part is 0
-                return replace(params, Xs=slv_Xs, m=oo), 1
+                return params.updated(Xs=slv_Xs, m=oo), 1
         else:
             # distribution is i.i.d. for each coordinate
-            red_Xs = replace(params.Xs, n=params.n - zeta)
-            slv_Xs = replace(params.Xs, n=zeta)
+            red_Xs = params.Xs.resize(params.n - zeta)
+            slv_Xs = params.Xs.resize(zeta)
 
         c = red_Xs.stddev * params.q / params.Xe.stddev
 
@@ -91,7 +94,7 @@ def dual_reduce(
         # Compute new noise as in [INDOCRYPT:EspJouKha20]
         # ~ sigma_ = rho * red_Xs.stddev * delta ** (m_ + red_Xs.n) / c ** (m_ / (m_ + red_Xs.n))
         sigma_ = rho * red_Xs.stddev * delta**d / c ** (m_ / d)
-        slv_Xe = NoiseDistribution.DiscreteGaussian(params.q * sigma_)
+        slv_Xe = DiscreteGaussian(params.q * sigma_)
 
         slv_params = LWEParameters(
             n=zeta,
@@ -177,7 +180,7 @@ def cost(
 
         rep = 1
         if params.Xs.is_sparse:
-            h = params.Xs.get_hamming_weight(params.n)
+            h = params.Xs.hamming_weight
             probability = RR(prob_drop(params.n, h, zeta, h1))
             rep = prob_amplify(success_probability, probability)
         # don't need more samples to re-run attack, since we may
@@ -210,7 +213,7 @@ def fft_solver(params, success_probability, t=0):
         probability = sqrt(success_probability)
 
         try:
-            size = params.Xs.support_size(n=params.n, fraction=probability)
+            size = params.Xs.support_size(probability)
             size_fft = 2**t
         except NotImplementedError:
             # not achieving required probability with search space
@@ -367,7 +370,7 @@ def __call__(
 
             >>> from estimator import *
             >>> from estimator.lwe_dual import dual_hybrid
-            >>> params = LWE.Parameters(n=1024, q = 2**32, Xs=ND.Uniform(0,1), Xe=ND.DiscreteGaussian(3.0))
+            >>> params = LWE.Parameters(n=1024, q = 2**32, Xs=ND.Binary, Xe=ND.DiscreteGaussian(3.0))
             >>> LWE.dual(params)
             rop: ≈2^107.0, mem: ≈2^66.4, m: 970, β: 264, d: 1994, ↻: 1, tag: dual
             >>> dual_hybrid(params)
@@ -377,13 +380,13 @@ def __call__(
             >>> dual_hybrid(params, mitm_optimization="numerical")
             rop: ≈2^129.0, m: 1145, k: 1, mem: ≈2^131.0, ↻: 1, β: 346, d: 2044, ζ: 125, tag: dual_mitm_hybrid
 
-            >>> params = params.updated(Xs=ND.SparseTernary(params.n, 32))
+            >>> params = params.updated(Xs=ND.SparseTernary(32))
             >>> LWE.dual(params)
             rop: ≈2^103.4, mem: ≈2^63.9, m: 904, β: 251, d: 1928, ↻: 1, tag: dual
             >>> dual_hybrid(params)
-            rop: ≈2^92.1, mem: ≈2^78.2, m: 716, β: 170, d: 1464, ↻: 1989, ζ: 276, h1: 8, tag: dual_hybrid
+            rop: ≈2^91.6, mem: ≈2^77.2, m: 711, β: 168, d: 1456, ↻: ≈2^11.2, ζ: 279, h1: 8, tag: dual_hybrid
             >>> dual_hybrid(params, mitm_optimization=True)
-            rop: ≈2^98.2, mem: ≈2^78.6, m: 728, k: 292, ↻: ≈2^18.7, β: 180, d: 1267, ζ: 485, h1: 17, tag: ...
+            rop: ≈2^98.7, mem: ≈2^78.6, m: 737, k: 288, ↻: ≈2^19.6, β: 184, d: 1284, ζ: 477, h1: 17, tag: dual_mitm_...
 
             >>> params = params.updated(Xs=ND.CenteredBinomial(8))
             >>> LWE.dual(params)
@@ -402,7 +405,7 @@ def __call__(
             rop: ≈2^160.7, mem: ≈2^156.8, m: 1473, k: 25, ↻: 1, β: 456, d: 2472, ζ: 25, tag: dual_mitm_hybrid
 
             >>> dual_hybrid(schemes.NTRUHPS2048509Enc)
-            rop: ≈2^131.7, mem: ≈2^128.5, m: 436, β: 358, d: 906, ↻: 1, ζ: 38, tag: dual_hybrid
+            rop: ≈2^136.2, mem: ≈2^127.8, m: 434, β: 356, d: 902, ↻: 35, ζ: 40, h1: 19, tag: dual_hybrid
 
             >>> LWE.dual(schemes.CHHS_4096_67)
             rop: ≈2^206.9, mem: ≈2^137.5, m: ≈2^11.8, β: 616, d: 7779, ↻: 1, tag: dual
@@ -440,7 +443,7 @@ def _optimize_blocksize(
                 log_level=None,
                 fft=False,
             ):
-                h = params.Xs.get_hamming_weight(params.n)
+                h = params.Xs.hamming_weight
                 h1_min = max(0, h - (params.n - zeta))
                 h1_max = min(zeta, h)
                 if h1_min == h1_max:

estimator/lwe_guess.py

@@ -15,7 +15,7 @@
 from .lwe_parameters import LWEParameters
 from .prob import amplify as prob_amplify, drop as prob_drop, amplify_sigma
 from .util import local_minimum, log2
-from .nd import sigmaf
+from .nd import sigmaf, SparseTernary
 
 
 class guess_composition:
@@ -102,7 +102,7 @@ def sparse_solve(cls, f, params, log_level=5, **kwds):
         :param params: LWE parameters.
         """
         base = params.Xs.bounds[1] - params.Xs.bounds[0]  # we exclude zero
-        h = ceil(len(params.Xs) * params.Xs.density)  # nr of non-zero entries
+        h = params.Xs.hamming_weight
 
         with local_minimum(0, params.n - 40, log_level=log_level) as it:
             for zeta in it:
@@ -127,13 +127,13 @@ def __call__(self, params, log_level=5, **kwds):
 
             >>> from estimator import *
             >>> from estimator.lwe_guess import guess_composition
-            >>> guess_composition(LWE.primal_usvp)(schemes.Kyber512.updated(Xs=ND.SparseTernary(512, 16)))
-            rop: ≈2^99.4, red: ≈2^99.4, δ: 1.008705, β: 113, d: 421, tag: usvp, ↻: ≈2^37.5, ζ: 265, |S|: 1, ...
+            >>> guess_composition(LWE.primal_usvp)(schemes.Kyber512.updated(Xs=ND.SparseTernary(16)))
+            rop: ≈2^102.2, red: ≈2^102.2, δ: 1.008011, β: 132, d: 461, tag: usvp, ↻: ≈2^34.9, ζ: 252, |S|: 1, ...
 
         Compare::
 
-            >>> LWE.primal_hybrid(schemes.Kyber512.updated(Xs=ND.SparseTernary(512, 16)))
-            rop: ≈2^85.8, red: ≈2^84.8, svp: ≈2^84.8, β: 105, η: 2, ζ: 366, |S|: ≈2^85.1, d: 315, prob: ≈2^-23.4, ...
+            >>> LWE.primal_hybrid(schemes.Kyber512.updated(Xs=ND.SparseTernary(16)))
+            rop: ≈2^85.8, red: ≈2^84.8, svp: ≈2^84.8, β: 105, η: 2, ζ: 366, |S|: ≈2^85.1, d: 315, prob: ≈2^-23.4, ↻:...
 
         """
         params = LWEParameters.normalize(params)
@@ -161,12 +161,12 @@ def __call__(self, params: LWEParameters, success_probability=0.99, quantum: boo
 
             >>> from estimator import *
             >>> from estimator.lwe_guess import exhaustive_search
-            >>> params = LWE.Parameters(n=64, q=2**40, Xs=ND.UniformMod(2), Xe=ND.DiscreteGaussian(3.2))
+            >>> params = LWE.Parameters(n=64, q=2**40, Xs=ND.Binary, Xe=ND.DiscreteGaussian(3.2))
             >>> exhaustive_search(params)
             rop: ≈2^73.6, mem: ≈2^72.6, m: 397.198
-            >>> params = LWE.Parameters(n=1024, q=2**40, Xs=ND.SparseTernary(n=1024, p=32), Xe=ND.DiscreteGaussian(3.2))
+            >>> params = LWE.Parameters(n=1024, q=2**40, Xs=ND.SparseTernary(32), Xe=ND.DiscreteGaussian(3.2))
             >>> exhaustive_search(params)
-            rop: ≈2^417.3, mem: ≈2^416.3, m: ≈2^11.2
+            rop: ≈2^413.9, mem: ≈2^412.9, m: ≈2^11.1
 
         """
         params = LWEParameters.normalize(params)
@@ -175,7 +175,7 @@ def __call__(self, params: LWEParameters, success_probability=0.99, quantum: boo
         probability = sqrt(success_probability)
 
         try:
-            size = params.Xs.support_size(n=params.n, fraction=probability)
+            size = params.Xs.support_size(probability)
         except NotImplementedError:
             # not achieving required probability with search space
             # given our settings that means the search space is huge
@@ -221,7 +221,7 @@ def X_range(self, nd):
         else:
             # setting fraction=0 to ensure that support size does not
             # throw error. we'll take the probability into account later
-            rng = nd.support_size(n=1, fraction=0.0)
+            rng = nd.resize(1).support_size(0.0)
             return rng, nd.gaussian_tail_prob
 
     def local_range(self, center):
@@ -240,11 +240,16 @@ def mitm_analytical(self, params: LWEParameters, success_probability=0.99):
         # about 3x faster and reasonably accurate
 
         if params.Xs.is_sparse:
-            h = params.Xs.get_hamming_weight(n=params.n)
-            split_h = round(h * k / n)
-            success_probability_ = (
-                binomial(k, split_h) * binomial(n - k, h - split_h) / binomial(n, h)
-            )
+            h = params.Xs.hamming_weight
+            if type(params.Xs) is SparseTernary:
+                # split optimally and compute the probability of this event
+                success_probability_ = params.Xs.split_probability(k)
+            else:
+                split_h = (h * k / n).round('down')
+                # Assume each coefficient is sampled i.i.d.:
+                success_probability_ = (
+                    binomial(k, split_h) * binomial(n - k, h - split_h) / binomial(n, h)
+                )
 
             logT = RR(h * (log2(n) - log2(h) + log2(sd_rng - 1) + log2(e))) / (2 - delta)
             logT -= RR(log2(h) / 2)
@@ -279,16 +284,20 @@ def cost(
         n = params.n
 
         if params.Xs.is_sparse:
-            h = params.Xs.get_hamming_weight(n=n)
-
+            h = params.Xs.hamming_weight
             # we assume the hamming weight to be distributed evenly across the two parts
             # if not we can rerandomize on the coordinates and try again -> repeat
-            split_h = round(h * k / n)
-            size_tab = RR((sd_rng - 1) ** split_h * binomial(k, split_h))
-            size_sea = RR((sd_rng - 1) ** (h - split_h) * binomial(n - k, h - split_h))
-            success_probability_ = (
-                binomial(k, split_h) * binomial(n - k, h - split_h) / binomial(n, h)
-            )
+            if type(params.Xs) is SparseTernary:
+                sec_tab, sec_sea = params.Xs.split_balanced(k)
+                size_tab = sec_tab.support_size()
+                size_sea = sec_sea.support_size()
+            else:
+                # Assume each coefficient is sampled i.i.d.:
+                split_h = (h * k / n).round('down')
+                size_tab = RR((sd_rng - 1) ** split_h * binomial(k, split_h))
+                size_sea = RR((sd_rng - 1) ** (h - split_h) * binomial(n - k, h - split_h))
+
+            success_probability_ = size_tab * size_sea / params.Xs.support_size()
         else:
             size_tab = sd_rng**k
             size_sea = sd_rng ** (n - k)
@@ -338,16 +347,16 @@ def __call__(self, params: LWEParameters, success_probability=0.99, optimization
 
             >>> from estimator import *
             >>> from estimator.lwe_guess import mitm
-            >>> params = LWE.Parameters(n=64, q=2**40, Xs=ND.UniformMod(2), Xe=ND.DiscreteGaussian(3.2))
+            >>> params = LWE.Parameters(n=64, q=2**40, Xs=ND.Binary, Xe=ND.DiscreteGaussian(3.2))
             >>> mitm(params)
             rop: ≈2^37.0, mem: ≈2^37.2, m: 37, k: 32, ↻: 1
             >>> mitm(params, optimization="numerical")
             rop: ≈2^39.2, m: 36, k: 32, mem: ≈2^39.1, ↻: 1
-            >>> params = LWE.Parameters(n=1024, q=2**40, Xs=ND.SparseTernary(n=1024, p=32), Xe=ND.DiscreteGaussian(3.2))
+            >>> params = LWE.Parameters(n=1024, q=2**40, Xs=ND.SparseTernary(32), Xe=ND.DiscreteGaussian(3.2))
             >>> mitm(params)
-            rop: ≈2^215.4, mem: ≈2^210.2, m: ≈2^13.1, k: 512, ↻: 43
+            rop: ≈2^217.8, mem: ≈2^210.2, m: ≈2^15.5, k: 512, ↻: 226
             >>> mitm(params, optimization="numerical")
-            rop: ≈2^216.0, m: ≈2^13.1, k: 512, mem: ≈2^211.4, ↻: 43
+            rop: ≈2^215.6, m: ≈2^15.5, k: 512, mem: ≈2^208.6, ↻: 226
 
         """
         Cost.register_impermanent(rop=True, mem=False, m=True, k=False)
@@ -400,7 +409,7 @@ def __call__(self, params: LWEParameters, success_probability=0.99):
 
             >>> from estimator import *
             >>> from estimator.lwe_guess import distinguish
-            >>> params = LWE.Parameters(n=0, q=2 ** 32, Xs=ND.UniformMod(2), Xe=ND.DiscreteGaussian(2 ** 32))
+            >>> params = LWE.Parameters(n=0, q=2 ** 32, Xs=ND.Binary, Xe=ND.DiscreteGaussian(2 ** 32))
             >>> distinguish(params)
             rop: ≈2^60.0, mem: ≈2^60.0, m: ≈2^60.0