# -*- coding: utf-8 -*-
from dataclasses import dataclass
from sage.all import binomial, ceil, exp, log, oo, parent, pi, RealField, RR, sqrt
[docs]def stddevf(sigma):
"""
Gaussian width parameter σ → standard deviation.
:param sigma: Gaussian width parameter σ
EXAMPLE::
>>> from estimator.nd import stddevf
>>> stddevf(64.0)
25.532...
>>> stddevf(64)
25.532...
>>> stddevf(RealField(256)(64)).prec()
256
"""
try:
prec = parent(sigma).prec()
except AttributeError:
prec = 0
if prec > 0:
FF = parent(sigma)
else:
FF = RR
return FF(sigma) / FF(sqrt(2 * pi))
[docs]def sigmaf(stddev):
"""
Standard deviation → Gaussian width parameter σ.
:param stddev: standard deviation
EXAMPLE::
>>> from estimator.nd import stddevf, sigmaf
>>> n = 64.0
>>> sigmaf(stddevf(n))
64.000...
>>> sigmaf(RealField(128)(1.0))
2.5066282746310005024157652848110452530
>>> sigmaf(1.0)
2.506628274631...
>>> sigmaf(1)
2.506628274631...
"""
RR = parent(stddev)
# check that we got ourselves a real number type
try:
if abs(RR(0.5) - 0.5) > 0.001:
RR = RealField(53) # hardcode something
except TypeError:
RR = RealField(53) # hardcode something
return RR(sqrt(2 * pi)) * stddev
[docs]@dataclass
class NoiseDistribution:
"""
All noise distributions are instances of this class.
"""
# cut-off for Gaussian distributions
gaussian_tail_bound = 2
# probability that a coefficient falls within the cut-off
gaussian_tail_prob = 1 - 2 * exp(-4 * pi)
stddev: float
mean: float = 0
n: int = None
bounds: tuple = (None, None)
density: float = 1.0 # Hamming weight / dimension.
tag: str = ""
def __lt__(self, other):
"""
We compare distributions by comparing their standard deviation.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.DiscreteGaussian(2.0) < ND.CenteredBinomial(18)
True
>>> ND.DiscreteGaussian(3.0) < ND.CenteredBinomial(18)
False
>>> ND.DiscreteGaussian(4.0) < ND.CenteredBinomial(18)
False
"""
try:
return self.stddev < other.stddev
except AttributeError:
return self.stddev < other
def __le__(self, other):
"""
We compare distributions by comparing their standard deviation.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.DiscreteGaussian(2.0) <= ND.CenteredBinomial(18)
True
>>> ND.DiscreteGaussian(3.0) <= ND.CenteredBinomial(18)
True
>>> ND.DiscreteGaussian(4.0) <= ND.CenteredBinomial(18)
False
"""
try:
return self.stddev <= other.stddev
except AttributeError:
return self.stddev <= other
def __str__(self):
"""
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.DiscreteGaussianAlpha(0.01, 7681)
D(σ=30.64)
"""
if self.n:
return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f}, n={int(self.n)})"
else:
return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f})"
def __repr__(self):
if self.mean == 0.0:
return f"D(σ={float(self.stddev):.2f})"
else:
return f"D(σ={float(self.stddev):.2f}, μ={float(self.mean):.2f})"
def __hash__(self):
"""
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> hash(ND(3.0, 1.0)) == hash((3.0, 1.0, None))
True
"""
return hash((self.stddev, self.mean, self.n))
def __len__(self):
"""
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> D = ND.SparseTernary(1024, p=128, m=128)
>>> len(D)
1024
>>> int(round(len(D) * float(D.density)))
256
"""
if self.n is not None:
return self.n
else:
raise ValueError("Distribution has no length.")
@property
def is_Gaussian_like(self):
return ("Gaussian" in self.tag) or ("CenteredBinomial" in self.tag)
@property
def is_bounded(self):
return (self.bounds[1] - self.bounds[0]) < oo
@property
def is_sparse(self):
"""
We consider a distribution "sparse" if its density is < 1/2.
"""
# NOTE: somewhat arbitrary
return self.density < 0.5
[docs] def support_size(self, n=None, fraction=1.0):
"""
Compute the size of the support covering the probability given as fraction.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> D = ND.Uniform(-3,3, 64)
>>> D.support_size(fraction=.99)
1207562882759477428726191443614714994252339953407098880
>>> D = ND.SparseTernary(64, 8)
>>> D.support_size()
32016101348447354880
"""
if not n:
if not self.n:
raise ValueError(f"Length required to determine support size, but n was {n}.")
n = self.n
if "SparseTernary" in self.tag:
h = self.get_hamming_weight(n)
# TODO: this is assuming that the non-zero entries are uniform over {-1,1}
# need p and m for more accurate calculation
size = 2**h * binomial(n, h) * RR(fraction)
elif self.is_bounded:
# TODO: this might be suboptimal/inaccurate for binomial distribution
a, b = self.bounds
size = RR(fraction) * (b - a + 1) ** n
else:
# Looks like nd is Gaussian
# -> we'll treat it as bounded (with failure probability)
t = self.gaussian_tail_bound
p = self.gaussian_tail_prob
if p**n < fraction:
raise NotImplementedError(
f"TODO(nd.support-size): raise t. {RR(p ** n)}, {n}, {fraction}"
)
b = 2 * t * sigmaf(self.stddev) + 1
return (2 * b + 1) ** n
return ceil(size)
[docs] def get_hamming_weight(self, n=None):
if not n:
if not self.n:
raise ValueError("Length required to determine hamming weight.")
n = self.n
return round(n * float(self.density))
[docs] @staticmethod
def DiscreteGaussian(stddev, mean=0, n=None):
"""
A discrete Gaussian distribution with standard deviation ``stddev`` per component.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.DiscreteGaussian(3.0, 1.0)
D(σ=3.00, μ=1.00)
"""
b_val = oo if n is None else ceil(log(n, 2) * stddev)
return NoiseDistribution(
stddev=RR(stddev),
mean=RR(mean),
n=n,
bounds=(-b_val, b_val),
density=1 - min(RR(1 / (sqrt(2 * pi) * stddev)), 1.0),
tag="DiscreteGaussian",
)
[docs] @staticmethod
def DiscreteGaussianAlpha(alpha, q, mean=0, n=None):
"""
A discrete Gaussian distribution with standard deviation α⋅q/√(2π) per component.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.DiscreteGaussianAlpha(0.001, 2048)
D(σ=0.82)
"""
stddev = stddevf(alpha * q)
return NoiseDistribution.DiscreteGaussian(stddev=RR(stddev), mean=RR(mean), n=n)
[docs] @staticmethod
def CenteredBinomial(eta, n=None):
"""
Sample a_1, …, a_η, b_1, …, b_η and return Σ(a_i - b_i).
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.CenteredBinomial(8)
D(σ=2.00)
"""
stddev = sqrt(eta / 2.0)
return NoiseDistribution(
stddev=RR(stddev),
density=1 - binomial(2 * eta, eta) * 2 ** (-2 * eta),
mean=RR(0),
n=n,
bounds=(-eta, eta),
tag="CenteredBinomial",
)
[docs] @staticmethod
def SparseTernary(n, p, m=None):
"""
Distribution of vectors of length ``n`` with ``p`` entries of 1 and ``m`` entries of -1, rest 0.
EXAMPLE::
>>> from estimator.nd import NoiseDistribution as ND
>>> ND.SparseTernary(100, p=10)
D(σ=0.45)
>>> ND.SparseTernary(100, p=10, m=10)
D(σ=0.45)
>>> ND.SparseTernary(100, p=10, m=8)
D(σ=0.42, μ=0.02)
"""
if m is None:
m = p
if n == 0:
# this might happen in the dual attack
return NoiseDistribution(
stddev=0, mean=0, density=0, bounds=(-1, 1), tag="SparseTernary", n=0
)
mean = RR(p / n - m / n)
stddev = sqrt(p / n * (1 - mean)**2 +
m / n * (-1 - mean)**2 +
(n - (p + m)) / n * (mean)**2)
density = RR((p + m) / n)
return NoiseDistribution(
stddev=stddev, mean=mean, density=density, bounds=(-1, 1), tag="SparseTernary", n=n
)