Improved Bernoulli Sampling for Discrete Gaussian Distributions over the Integers

Abstract

Discrete Gaussian sampling is one of the fundamental mathematical tools for lattice-based cryptography. In this paper, we revisit the Bernoulli(-type) sampling for centered discrete Gaussian distributions over the integers, which was proposed by Ducas et al. in 2013. Combining the idea of Karney’s algorithm for sampling from the Bernoulli distribution ${\mathcal{B}}_{e^{-1/2}}$, we present an improved Bernoulli sampling algorithm. It does not require the use of floating-point arithmetic to generate a precomputed table, as the original Bernoulli sampling algorithm did. It only needs a fixed look-up table of very small size (e.g., 128 bits) that stores the binary expansion of $\ln 2$. We also propose a noncentered version of Bernoulli sampling algorithm for discrete Gaussian distributions with varying centers over the integers. It requires no floating-point arithmetic and can support centers of precision up to 52 bits. The experimental results show that our proposed algorithms have a significant improvement in the sampling efficiency as compared to other rejection algorithms.

1. Introduction

Lattice-based cryptography is one of the most widely and deeply studied post-quantum cryptography. Discrete Gaussian sampling is one of the important mathematical tools for lattice-based cryptography [1,2,3]. Many lattice-based cryptosystems, including their security proofs, involve discrete Gaussian sampling.

Generally, discrete Gaussian sampling is to sample from a (multivariate) discrete Gaussian distribution over an n-dimensional lattice $\mathsf{\Lambda }$. In particular, discrete Gaussian sampling over the integers, denoted by $\mathrm{Sample}\mathbb{Z}$, refers to the generation of random integers that conform to Gaussian distribution ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with parameter $\sigma >0$ and center $\mu \in \mathbb{R}$. Since $\mathrm{Sample}\mathbb{Z}$ is much more efficient and simpler than discrete Gaussian sampling over a general lattice, many lattice-based cryptosystems, such as [4,5], only involve discrete Gaussian sampling over the integers.

The methods of sampling from a continuous Gaussian distribution are not trivially applicable for the discrete case. Therefore, it is of great significance to study discrete Gaussian sampling and design sampling algorithms, especially for $\mathrm{Sample}\mathbb{Z}$. This issue has attracted a lot of attention in recent years.

1.1. Related Work

For the commonly used methods (techniques) for $\mathrm{Sample}\mathbb{Z}$, including inversion sampling, the Knuth-Yao method and the rejection sampling, we can refer to [6,7]. One can also consider the convolution technique (based on the convolution-like properties of discrete Gaussian developed by Peikert in [8]) [9,10]. In addition, Saarinen developed a new binary arithmetic coding sampler for discrete Gaussians (BAC Sampling) [11]. Moreover, Karmakar et al. designed a class of very efficient discrete Gaussian samplers based on computing a series of Boolean functions with a number of random bits as the input [12,13].

Sampling from a discrete Gaussian distribution over the integers usually involves floating-point arithmetic because of the exponential computation in Gaussian function. It is believed that high-precision floating-point arithmetic is required for attaining enough statistical distance to the target distribution. This implies that online high-precision floating-point computation may be an implementation bottleneck. Thus, we can see that most of the improved $\mathrm{Sample}\mathbb{Z}$ algorithms manage to use high-precision floating-point arithmetic only at offline time (such as algorithms based on CDT’s or Knuth-Yao) or even to avoid floating-point arithmetic (such as Karney’s sampling algorithm [14]).

Recently results show that the Bernoulli sampling, introduced by Ducas et al. in [4] for centered discrete Gaussian distributions over the integers, can be effectively improved by combining the polynomial approximation technique. Specifically, Bernoulli sampling is a kind of rejection sampling. Its rejection operation can be done very efficiently by calculating an approximated polynomial [15,16]. The polynomial is predetermined; it thus supports discrete Gaussians with a varying center $\mu \in [0,1)$ over the integers.

Generic discrete Gaussian sampling algorithms over the integers that support a varying center $\mu$ and even a varying parameter $\sigma$ are also very interesting. A discrete Gaussian sampling algorithm for an n-dimensional lattice usually requires such an algorithm as its subroutine [2,8,17,18]. However, most of the algorithms for sampling ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ and their improvements are designed only for the case where center $\mu$ is fixed in advance, such as [4,9,11,19,20]. When sampling from ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with varying $\mu \in [0,1)$ (e.g., double-precision numbers that are sampled from a specific distribution online), except [14], the existing algorithms available rely on floating-point arithmetic, such as the algorithms in [2,15] and the sampler used in Falcon signature [5] or require a large amount of offline computation and precomputation storage, such as the algorithms proposed by Micciancio et al. [10] and by Aguilar-Melchor et al. [21] respectively. Moreover, the sampler presented by Barthe et al. in [16] is not a generic algorithm, as it supports a varying $\mu \in \mathbb{R}$ but not a varying parameters $\sigma$. One needs to determine a new approximation polynomial for each new $\sigma$.

1.2. Our Contribution

In this paper, combining the idea of Karney’s algorithm for sampling the Bernoulli distribution ${\mathcal{B}}_{e^{-1/2}}$, we present a generalized version of Bernoulli sampling algorithm for both centered and noncentered discrete Gaussian distributions over the integers (see Algorithm 1).

Specifically, we present an improved implementation of the Bernoulli sampling algorithm for centered discrete Gaussian distributions over the integers (see Algorithm 2). We remove the precomputed table, that is required and obtained by using floating-point arithmetic in advance in the original Bernoulli sampling algorithm. Except a fixed look-up table of very small size (e.g., 128 bits) that stores the binary expansion of $\ln 2$, our proposed algorithm requires no precomputation storage and uses no floating-point arithmetic.

We also propose a noncentered Bernoulli sampling algorithm (see Algorithm 3) that supports discrete Gaussian distributions with varying centers. It uses no floating-point arithmetic, and it is suitable for sampling ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with $\mu \in [0,1)$ of precision up to 52 bits. We also show that our noncentered version of sampling algorithm has the same rejection rate as the original Bernoulli sampling algorithm. Furthermore, the experimental results show that our proposed algorithms have a significant improvement in the sampling efficiency as compared to other rejection algorithms for discrete Gaussian distributions over the integers [14].

Table 1. Comparison of sampling algorithms over the integers.

Algorithm	Memory	FPA.	$\exp (·)$	Varying $\mathit{\mu }$	Varying-$\mathit{\sigma }$
Rejection Sampling [2]	0	Yes	Yes	Yes	Yes
Inversion Sampling [8]	$O\left(\lambda \sigma \right)$	No	No	No	No
Discrete Ziggurat [22]	var.	Yes	Yes	No	No
Bernoulli Sampling [4]	$O\left(\lambda {log}_2\sigma \right)$	No	No	No	Yes
Knuth-Yao [7]	$O\left(\lambda \sigma \right)$	No	No	No	No
Karney [14]	0	No	No	Yes	Yes
Convolution-like [10]	var.	Yes	No	Yes	Yes
BAC Sampling [11]	$O\left(\lambda \sigma \right)$	No	No	No	No
Karmakar [12,13]	0	No	No	No	No
poly. approximation [15]	$O\left(\lambda \right)$	Yes	No	Yes	Yes
poly. approximation [16]	$O\left(\lambda \right)$	No	No	No	No
our Algorithm 2	$O\left(1\right)$	No	No	No	Yes
our Algorithm 3	$O\left(1\right)$	No	No	Yes	Yes

summarizes the some existing sampling algorithms for discrete Gaussian distributions over the integers and their properties in comparison to our work. Here, $\lambda$ is the security parameters. The column “FPA” and the column “$\exp (·)$” indicate if the algorithm requires to use floating-point arithmetic online and to evaluate $\exp (·)$ online respectively. The column “varying $\mu$” refers to the property of being able to produce samples (online) from discrete Gaussians with varying centers rather than a specific discrete Gaussian with a fixed center. Similarly, the column “varying $\sigma$” refers to the property of being able to produce samples (online) from discrete Gaussians with varying parameter $\sigma$.

2. Materials and Methods

2.1. Preliminaries

In this paper, $\mathbb{R}$ denotes the set of real numbers, $\mathbb{Z}$ and ${\mathbb{Z}}^{+}$ denotes the set of integers and non-negative integers respectively. We use $f\left(S\right)={\sum }_{x\in S}f\left(x\right)$ to define any real function $f(·)$ to a countable set S if it exists, and use ${\rho }_{\sigma ,\mu }\left(x\right)=\exp \left(-\frac{{(x-\mu )}^2}{2{\sigma }^2}\right)$ to define the Gaussian function evaluated at $x\in \mathbb{R}$, where parameter $\sigma >0$ and center $\mu \in \mathbb{R}$. Besides, ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }={\rho }_{\sigma ,\mu }\left(x\right)/{\rho }_{\sigma ,\mu }\left(\mathbb{Z}\right)$ defines the discrete Gaussian distribution evaluated at $x\in \mathbb{Z}$, where $x\in \mathbb{Z}$, $\mu \in \mathbb{R}$, $\sigma >0$.

Bernoulli(-type) sampling algorithm proposed by Ducas et al. in [4] was designed exclusive for a centered discrete Gaussian distribution over the integers, namely ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with $\sigma =k{\sigma }_2$ and $\mu =0$, where k is a positive integer and ${\sigma }_2=\sqrt{1/(2·\ln 2)}$. This algorithm uses rejection sampling from the binary discrete Gaussian distribution, denoted by ${\mathcal{D}}_{{\mathbb{Z}}^{+},{\sigma }_2}$. Sampling from ${\mathcal{D}}_{{\mathbb{Z}}^{+},{\sigma }_2}$ can be very efficient using only unbiased random bits (one may refer to Algorithm 10 in [4]). Firstly, it samples an integer $x\in {\mathbb{Z}}^{+}$ from ${\mathcal{D}}_{{\mathbb{Z}}^{+},{\sigma }_2}$, then samples $y\in \mathbb{Z}$ uniformly in $\{0,1,2,\cdots ,k-1\}$ and accepts $z=kx+y$ with probability $\exp (-(y^2+2kxy)/2k^2{\sigma }_2^2)$, finally, it outputs a signed integer, z or $-z$, with equal probabilities because of the symmetry of the target distribution (see Algorithm 11 and Algorithm 12 in [4]). It was shown that the average number of rejections is not more than $1.47$.

A precomputed look-up table is required for computing the exponential function $\exp (-(y^2+2kxy)/2{\sigma }^2)$. The size of the table varies with the parameter $\sigma$. For examples, the size is $2.3$ kb for $\sigma \approx 107$ and 4 kb for $\sigma \approx 271$ according to the estimations given by Ducas et al. in [4], which is not very small but acceptable.

Algorithm 1 is adapted by Karney [14] from Von Neumann’s algorithm that samples from $e^{-x}$, where $x>0$. Algorithm 1 is used as one of the basic operations in Karney’s sampling algorithm for the standard normal distribution, which can be executed only by using integer arithmetic.

In Algorithm 1, by replacing $1/2$ with $p/q$, we can get a Bernoulli random value. The value is true with probability $e^{-p/q}$, namely sampling from ${\mathcal{B}}_{e^{-p/q}}$, where $p,q$ are positive integers and $p<q$. The above modification can be regarded as a generalized version of Algorithm 1, which is used implicitly in Karney’s algorithm. The discrete Gaussian distribution over integers in Karney’s algorithm also adopts the rejection sampling algorithm, which is a discretization of Karney’s algorithm for sampling from the standard normal distribution. It is fit for discrete Gaussians with varying rational parameters of precision up to 32 bits (including $\sigma$ and $\mu$). The algorithm does not require any floating-point arithmetic. Moreover, it needs no precomputation storage. One may refer to [14] for the full description of Karney’s sampling algorithm.

Algorithm 1: [14] Sample from Bernoulli distribution ${\mathcal{B}}_{1/\sqrt{e}}$

2.2. Removing Precomputed Tables in Bernoulli Sampling

One drawback of the Bernoulli sampling is that it needs a precomputed table to generate a Bernoulli random value from ${\mathcal{B}}_{\exp (-(y^2+2kxy)/2{\sigma }^2)}$. Although the table is not large, it has to be precomputed before sampling for every different $\sigma =k{\sigma }_2$. Therefore, it is interesting to sample ${\mathcal{D}}_{\mathbb{Z},k{\sigma }_2}$ without a precomputed table. In this subsection, we propose an alternative algorithm, which samples ${\mathcal{D}}_{\mathbb{Z},k{\sigma }_2}$ for a given $k\ge 1$ and does not require a precomputed table except a fixed look-up table with small size.

We use the notation of the sampling algorithm given by Ducas et al. (Algorithm 11 in [4]). Since ${\sigma }_2=\sqrt{1/(2·\ln 2)}$ and $\sigma =k{\sigma }_2$, it follows that $\exp (-(y^2+2kxy)/2{\sigma }^2)$ is equal to

$$\exp \left(-\ln 2\left(⌊\frac{y^2+2kxy}{k^2}⌋+\left\{\frac{y^2+2kxy}{k^2}\right\}\right)\right)$$

$$=\exp \left(-\ln 2⌊\frac{y^2+2kxy}{k^2}⌋\right)·\exp \left(-\ln 2\left\{\frac{y^2+2kxy}{k^2}\right\}\right)$$

where

$$\left\{\frac{y^2+2kxy}{k^2}\right\}=\frac{y^2+2kxy}{k^2}-⌊\frac{y^2+2kxy}{k^2}⌋=\frac{(y^2+2kxy)\mathrm{mod}{k}^2}{k^2}$$

is the fractional part of $(y^2+2kxy)/k^2$. Then the Bernoulli value can be obtained by generating $\lfloor y^2+2kxy/k^2\rfloor$ Bernoulli deviates according to ${\mathcal{B}}_{1/2}$ and one Bernoulli random value from ${\mathcal{B}}_{\exp (-\left(\ln 2\right)\{y^2+2kxy/k^2\})}$. If no false value is generated in this procedure, then we have true, otherwise we have false.

It is clear that a Bernoulli deviate according to ${\mathcal{B}}_{1/2}$ is equivalent to a uniform random bit. Thus, the only remaining problem is to generate a Bernoulli random value from ${\mathcal{B}}_{\exp (-\left(\ln 2\right)\{(y^2+2kxy)/k^2\})}$. We address this problem via Algorithm 2.

One can see that Algorithm 2 is adapted from Algorithm 1. It provides a technique for determining the relation between $\left(\ln 2\right)(z/k^2)$ and a uniform deviate $u_1\in [0,1)$, and it does not need floating-point operations during run-time, if the binary expansion of $\ln 2$, that is long enough (e.g., 128 bits), has been stored as a fixed look-up table.

In step 1, p is assigned to be the t most significant bits of the binary expansion of $\ln 2$. In steps 5 and 6, if $\lfloor p\rfloor >\lfloor q/z\rfloor =\lfloor \left(v{k}^2\right)/z\rfloor$, then we have $2^t\left(\ln 2\right)>\left(v{k}^2\right)/z$, which follows that $\left(\ln 2\right)(z/k^2)>v/2^t$. It is equivalent to obtaining a uniform deviate $u_1=(v/2^t)$ that is strictly less than the value of $\left(\ln 2\right)(z/k^2)$. On the contrary, in steps 7 and 8, if $\lfloor p\rfloor <\lfloor q/z\rfloor$, then it is equivalent to obtaining a uniform deviate $u_1$ that is strictly more than the value of $\left(\ln 2\right)(z/k^2)$. Otherwise, if $\lfloor p\rfloor =\lfloor q/z\rfloor$, i.e., $\lfloor 2^t\left(\ln 2\right)\rfloor =\lfloor \left(v{k}^2\right)/z\rfloor$. This means that we get a uniform deviate $u_1$, but in this moment we cannot determine whether it is less than the value of $\left(\ln 2\right)(z/k^2)$ or not.

Algorithm 2: Sample from Bernoulli distribution ${\mathcal{B}}_{\exp (-\left(\ln 2\right)(z/k^2))}$, where $k,z$ are positive integers such that $z<k^2$

In Algorithm 2, the while loop is designed to further compare the uniform deviate with the value of $\left(\ln 2\right)(z/k^2)$. Informally, if $\lfloor 2^t\left(\ln 2\right)\rfloor =\lfloor \left(v{k}^2\right)/z\rfloor$ for some $v_1\in [0,2^t)$, then we generate another uniformly random integer $v_2\in [0,2^t)$, which means that we get a uniform deviate $u_1\ge (v_1/2^t)+(v_2/2^{2t})$. When comparing $\left(\ln 2\right)(z/k^2)$ with $(v_1/2^t)+(v_2/2^{2t})$, if

$$2^t\left(\ln 2\right)>\frac{v_1{k}^2}{z}+\frac{v_2{k}^2}{z{2}^t}$$

then we have $\left(\ln 2\right)(z/k^2)>(v_1/2^t)+(v_2/2^{2t})$, which follows that $\left(\ln 2\right)(z/k^2)>u_1$. Since $\lfloor 2^t\left(\ln 2\right)\rfloor =\lfloor \left(v{k}^2\right)/z\rfloor$, we have $\left(\ln 2\right)(z/k^2)>u_1$ if

$$\lfloor 2^t\left(\ln 2\right)\rfloor +\left\{2^t\left(\ln 2\right)\right\}>\lfloor \frac{v_1{k}^2}{z}\rfloor +\left\{\frac{v_1{k}^2}{z}\right\}+\frac{v_2{k}^2}{z{2}^t},$$

i.e.,

$$\left\{2^t\left(\ln 2\right)\right\}>\left\{\frac{v_1{k}^2}{z}\right\}+\frac{v_2{k}^2}{z{2}^t},$$

which is equivalent to

$$2^t\left\{2^t\left(\ln 2\right)\right\}>2^t\left\{\frac{v_1{k}^2}{z}\right\}+\frac{v_2{k}^2}{z},$$

where $\{v_1{k}^2/z\}$ is the fractional part of the rational number $v_1{k}^2/z$, which can be determined by computing $v_1{k}^2\left(\mathrm{mod}z\right)$. Therefore, as presented in Algorithm 2, after steps 9, 10, and 11 what we get is actually $p=(p-\lfloor p\rfloor )·2^t=(2^t\left(\ln 2\right)-\lfloor (2^t\left(\ln 2\right)\rfloor )2^t=2^t\left\{2^t\left(\ln 2\right)\right\}$ and $q=(\left(v_1{k}^2\mathrm{mod}z\right)·2^t+v_2{k}^2)$. So, we can continue to compare the uniform deviate $u_1$ with the value of $\left(\ln 2\right)(z/k^2)$. For the full proof of the correctness, one may refer to Appendix A.

2.3. Bernoulli Sampling for Discrete Gaussians with Varying Centers

In this subsection, we propose a new noncentered Bernoulli sampling algorithm for ${\mathcal{D}}_{\mathbb{Z},\sigma ,i/b}$ with $\sigma =k{\sigma }_2$, where b is a power of two, $i\in [0,b-1]\cap \mathbb{Z}$, k is a positive integer, ${\sigma }_2=\sqrt{1/(2·\ln 2)}$. Our proposed algorithm does not require any floating-point arithmetic, and it needs no precomputed table except a fixed and small size look-up table for storing the binary expansion of $\ln 2$, as mentioned in Section 2.2.

Informally, the sketch of the proposed algorithm is as follows. We sample an integer $x\in {\mathbb{Z}}^{+}$ from ${\mathcal{D}}_{{\mathbb{Z}}^{+},{\sigma }_2}$ firstly, and then sample $y\in \mathbb{Z}$ uniformly in $\{0,1,2,\cdots ,k\}$. Since the symmetry of a noncentered discrete Gaussian distribution does not hold, we need to determine the sign of prospective output before the acceptance-rejection operation, by setting $s\leftarrow \pm 1$ with equal probabilities. Finally, we accept $s(kx+y)$ as the output with probability

$$\exp \left(-\frac{{(y-s\mu )}^2+2kx(y-s\mu )}{2k^2{\sigma }_2^2}\right),$$

and restart the whole algorithm if otherwise, where $\mu =i/b$. The following equation

$$\exp \left(-\frac{x^2}{2{\sigma }_2^2}\right)·\exp \left(-\frac{{(y-s\mu )}^2+2kx(y-s\mu )}{2k^2{\sigma }_2^2}\right)=\exp \left(-\frac{{(s(kx+y)-\mu )}^2}{2k^2{\sigma }_2^2}\right)$$

guarantees that the output value has the desired relative probability density ${\rho }_{k{\sigma }_2,\mu }\left(s(kx+y)\right)$. The sketch of the algorithm seems to work, but there are several subtleties we have to address.

Bernoulli sampling for centered discrete Gaussian distributions needs to restart its whole sampling procedure with probability $1/2$ when the algorithm is about to return $z=0$ (see Algorithm 12 in [4]). In our algorithm, in order to avoid double counting 0, we need to set $y\leftarrow y+1$ if $s=1$, so that there is only one case where $z=0$. This also guarantees that $y-s\mu$ is always non-negative, and thus $\exp \left(-({(y-s\mu )}^2+2kx(y-s\mu ))/\left(2k^2{\sigma }_2^2\right)\right)$ is never more than one.

Another problem is how to compute the value of the exponential function efficiently. Since ${\sigma }_2=\sqrt{1/(2·\ln 2)}$ and $\mu =i/b$ with $i\in [0,b-1]\cap \mathbb{Z}$, it follows that

$$\exp \left(-\frac{{(y-s\mu )}^2+2kx(y-s\mu )}{2k^2{\sigma }_2^2}\right)=2^{-\lfloor \delta \rfloor }+2^{-\left\{\delta \right\}}$$

where $\delta =({(by-si)}^2+2bkx(by-si))/\left(b^2{k}^2\right)$ and

$$\left\{\delta \right\}=\delta -\lfloor \delta \rfloor =\frac{\left({(by-si)}^2+2bkx(by-si)\right)\mathrm{mod}\left(b^2{k}^2\right)}{b^2{k}^2}$$

is the fractional part of $\delta$. Then, when computing the value of exponential function $\exp \left(-({(y-s\mu )}^2+2kx(y-s\mu ))/\left(2k^2{\sigma }_2^2\right)\right)$, we generate $\lfloor \delta \rfloor$ Bernoulli deviates according to ${\mathcal{B}}_{1/2}$ and one Bernoulli random value which is true with probability $\exp (-\ln 2·\left\{\delta \right\})=2^{-\left\{\delta \right\}}$. If no false value is generated, then we accept $s(kx+y)$ as the output, otherwise we restart the algorithm. In summary, our proposed algorithm can be formally described as Algorithm 3, which summarizes the above discussion.

Recall that our goal is to give an algorithm that requires no floating-point arithmetic. It is clear that a Bernoulli deviate according to ${\mathcal{B}}_{1/2}$ is equivalent to a uniform random bit, and sampling from ${\mathcal{D}}_{{\mathbb{Z}}^{+},{\sigma }_2}$ can be efficiently implemented by using only unbiased random bits (Algorithm 10 in [4]). Thus, the only remaining problem is to generate a Bernoulli random value with probability

$$\exp (-\ln 2·\left\{\delta \right\})=\exp \left(-\ln 2·\left\{\frac{{(by-si)}^2+2bkx(by-si)}{b^2{k}^2}\right\}\right)$$

without using floating-point arithmetic. Since ${(by-si)}^2+2bkx(by-si)$ is a non-negative integer, what we need essentially is a solution for generating a Bernoulli random value from ${\mathcal{B}}_{\exp (-\left(\ln 2\right)(z/b^2{k}^2))}$, where $z={(by-si)}^2+2bkx(by-si)$ and k are positive integers such that $z<b^2{k}^2$. To end this, we can apply Algorithm 2 by replacing k with $b^2{k}^2$, which involves no floating-point arithmetic and requires no precomputation storage.

Algorithm 3: Sample from discrete Gaussian distribution ${\mathcal{D}}_{\mathbb{Z},k{\sigma }_2,i/b}$ for $k\in {\mathbb{Z}}^{+}$, $0\le i<b$ and $b\ge 2$

Moreover, let $k<2^8$ and $b=2^{20}$, then we have ${(by-si)}^2\le b^2{k}^2<2^{56}$ and $2bkx(by-si)<2^{60}$ with overwhelming probability. This means that Algorithm 3 supports a $\mu$ of precision up to 20 bits if k is restricted to be less than $2^8$, and it can be implemented only by using 64-bit integer arithmetic.

A sampler for centers of low precision can be used in q-ary lattices, where q is relatively large but not very large, and sampling from discrete Gaussian distributions with q different centers may be required in this case. It involves a large storage requirement if one uses precomputed tables. So, one could consider our sampler as described in Algorithm 3. It requires no precomputed table.

2.4. Increasing the Precision of $\mu$

Algorithm 3 only supports centers of finite precision. The precision of $\mu$ depends on the value of b. A larger b means a $\mu$ of higher precision. For example, let $b=2^{52}$. We remark that the bound of 52 bits is significant. We note that float-pointing with 53 bits of precision corresponds to the double precision type in the IEEE 754 standard, and it is usable in software. When $\mu \in [0,1)$ is given by a double-precision floating-point number, we could represent $\mu$ as a rational number of denominator $2^{52}$ and use Algorithm 3 directly.

However, b cannot be $2^{52}$ in practice, since in this case we need to compute

$$\frac{{(by-si)}^2+2bkx(by-si)}{b^2{k}^2},$$

where $b^2{k}^2>2^{64}$, which leads to the integer overflow (for 64-bit C/C++ compilers). A big integer library works but compromises the performance at runtime. In this subsection, therefore, we try to improve Algorithm 3 so that it supports a $\mu$ of precision up to 52 bits.

The basic idea is to avoid compute $b^2{k}^2$. We note that generating a Bernoulli random value which is true with probability

$$\exp \left(-\ln 2\left(\frac{{(by-si)}^2+2bkx(by-si)}{b^2{k}^2}\right)\right)$$

can be decomposed to generating two Bernoulli random values which are true with probability $\exp \left(-\ln 2{\left(\frac{(by-si)}{bk}\right)}^2\right)$ and $\exp \left(-\ln 2\left(\frac{2x(by-si)}{bk}\right)\right)$ respectively. It is not hard to see that a Bernoulli random value which is true with probability $\exp \left(-\ln 2\left(\frac{2x(by-si)}{bk}\right)\right)$ can be generated by using Algorithm 2 with $z=2x(by-si)$ and $k^2=bk$. In particular, Algorithm 2 requires $2x(by-si)<bk$. When $2x(by-si)>bk$, we need to compute $\lfloor 2x(by-si)/bk\rfloor$ and invoke Algorithm 1 firstly, and then use Algorithm 2 with $z=2x(by-si)\mathrm{mod}\left(bk\right)$ and $k^2=bk$.

The following algorithm is designed to generate a Bernoulli random value from Bernoulli distribution ${\mathcal{B}}_{\theta }$ with $\theta =\exp \left(-\ln 2{\left(\frac{z}{K}\right)}^2\right)$. As the inequality $by-si<bk$ holds except $i=0$ and $y=0$, $\exp \left(-\ln 2{\left(\frac{(by-si)}{bk}\right)}^2\right)$ can be calculated by using Algorithm 4 with $z=(by-si)$ and $K=bk$.

Algorithm 4: Sample from Bernoulli distribution ${\mathcal{B}}_{\theta }$ with $\theta =\exp \left(-\ln 2{\left(\frac{z}{K}\right)}^2\right)$ and $0<z<K$

The main idea behind Algorithm 4 is to sample two sets of uniform deviates $u_1,u_2,\dots$ and $v_1,v_2,\dots$, and to decide the maximum value $n\ge 0$ such that $\ln 2·(z/K)>u_1>u_2>\dots >u_n$ and $v_i<z/k$. Then, the probability that n is even, which means it returns true, is given by

$$\begin{array}{c}{\displaystyle 1-\frac{\left(\right(z/K)·\ln 2)}{1!}\left(\frac{z}{K}\right)+\frac{{((z/K)·\ln 2)}^2}{2!}{\left(\frac{z}{K}\right)}^2-\frac{{((z/K)·\ln 2)}^3}{3!}{\left(\frac{z}{K}\right)}^3+\dots }\hfill \\ =1-\frac{\ln 2}{1!}{\left(\frac{z}{K}\right)}^2+\frac{{\left(\ln 2\right)}^2}{2!}{\left(\frac{z}{K}\right)}^4-\frac{({\left(\ln 2\right)}^3}{3!}{\left(\frac{z}{K}\right)}^6+\dots =\exp \left(-\left(\ln 2\right){\left(\frac{z}{K}\right)}^2\right).\hfill \end{array}$$

Furthermore, when implementing Algorithm 4, we need to apply the technique presented in Algorithm 2 for comparing u with $(z/K)·\ln 2$ and avoiding using floating-point arithmetic.

Let us go back to the problem of increasing the precision of $\mu$. If k is restricted to be less than $2^8$, we can take $b=2^{52}$. In this case, we have $(by-si)\le bk<2^{60}$ and $2x(by-si)<2^{64}$ with overwhelming probability. This means that Algorithm 3, incorporated with Algorithm 4 as well as Algorithm 2, supports a $\mu$ of precision up to 52 bits if k is restricted to be less than $2^8$, and it can be implemented only by using 64-bit integer arithmetic.

3. Results

On a laptop computer (Intel i7-8550U, 16 GB RAM, the g++ compiler and enabling -O3 optimization option), we implemented GPV algorithm in [2], Karney’s algorithm in [14], the original Bernoulli sampling algorithm in [4], and our Algorithm 3 respectively. Except Karney’s algorithm, the uniformly pseudorandom numbers are generated by using AES256-CTR with AES-NI (Intel Advanced Encryption Standard New Instructions). Figure 1 shows the performance comparison of these sampling algorithms.

We select Karney’s algorithm and the GPV algorithm to compare with our generalized Bernoulli sampling because they also use rejection sampling and requires no precomputation storage. Although algorithms based on the inversion sampling or Knuth-Yao method are usually more efficient, they require a large amount precomputation storage, especially for sampling from distributions with variable $\mu$ or large $\sigma$.

For centered discrete Gaussian distributions with $\sigma$ roughly from 10 to 190, one can get about $12.7\times {10}^6$ samples per second by using Algorithm 3. This means that the performance of our improved Bernoulli sampling algorithm is significantly better than the original Bernoulli sampling as well as Karney’s sampling algorithm. When $b=2^{20}$, for the discrete Gaussian distribution ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with $\mu =i/b$ and i uniformly from $[0,b-1]$, we also tested that Algorithm 3 has almost the same performance as the case of centered discrete Gaussian distributions. When we take $b=2^{52}$, we need combine Algorithm 4 to implement Algorithm 3. The performance in this case decreases to about $6.4\times {10}^6$ per second and is slightly worse than the original Bernoulli sampling algorithm, but the latter does not support noncentered discrete Gaussian distributions.

Here, the GPV algorithm has the worst performance, since its proposal distribution is poor and its rejection rate is too high, which is close to 10 [2], while the Bernoulli sampling algorithm has a rejection rate of about $1.47$, and the rejection rate of Karney’s algorithm is even as lower as $1.39$ [14]. We will prove in Section 4 of this paper that the rejection rate of Algorithm 3 is still about $1.47$.

Our Algorithm 3 has a better performance than Karney’s algorithm, this is due to the fact that it uses the binary discrete Gaussian distribution ${\mathcal{D}}_{\mathbb{Z},{\sigma }_2}$ as its base sampler instead of ${\mathcal{D}}_{\mathbb{Z},1}$, and there exists a very efficient sampling algorithm for such a discrete Gaussian distribution. Algorithm 3 is faster than the original Bernoulli sampling algorithm, mainly because its rejection operation is simplified to sample from only one Bernoulli distribution with an exponential function as its parameter, while the original Bernoulli sampling algorithm requires a lookup table to sample a relatively large number of Bernoulli distributions with exponential functions as parameters in order to perform the rejection operation, which affects its practical performance.

Note that our Algorithm 4 requires Algorithm 2 to be fully implemented, and the basic idea of Algorithm 4 follows from that of Algorithm 2. We perform algorithm 4 once roughly as the same amount of computation as Algorithm 2 twice. This helps us intuitively explain why the performance of Algorithm 3 drops by about half after using Algorithm 4.

4. Discussion

In this section, we show that the noncentered version of Bernoulli sampling presented by Algorithm 3 has the same rejection rate as the Bernoulli sampling for centered discrete Gaussians.

Intuitively, evaluating the rejection rate of Algorithm 3 could be done just as evaluating the rejection rate of Bernoulli sampling for centered discrete Gaussian distributions, which has been presented by Ducas et al. in [4]. Notwithstanding, we may need the notion of the smoothing parameter of lattices, which was introduced by Micciancio and Regev in [1], because of the noncentered discrete Gaussian distributions. In fact, we only need the notion of smoothing parameter of lattice $\mathbb{Z}$ and its two important properties, since this work only involves the (one-dimensional) integer lattice $\mathbb{Z}$, which is the simplest lattice.

Definition 1

(a special case of [1], Definition 3.1). Let $ϵ>0$ be a positive real. The smoothing parameter of lattice $\mathbb{Z}$, denoted by ${\eta }_{ϵ}\left(\mathbb{Z}\right)$, is defined to be the smallest real s such that ${\rho }_{1/s}(\mathbb{Z}\backslash \left\{0\right\})\le ϵ$.

Lemma 1

(a special case of [1], Lemma 3.3). For any real $ϵ>0$, the smoothing parameter of lattice $\mathbb{Z}$ satisfies ${\eta }_{ϵ}\left(\mathbb{Z}\right)\le \sqrt{\ln \left(2\right(1+1/ϵ\left)\right)/2}/\pi$.

Lemma 2

(Ref. [1], Lemma 2.8, implicit). For any (small) real $ϵ>0$ and $\mu \in \mathbb{R}$, if $\sigma \ge {\eta }_{ϵ}\left(\mathbb{Z}\right)$, then ${\rho }_{\sigma ,\mu }\left(\mathbb{Z}\right)=\left[\frac{1-ϵ}{1+ϵ},1\right]·{\rho }_{\sigma }\left(\mathbb{Z}\right)$, and ${\rho }_{\sigma }\left(\mathbb{Z}\right)\approx {\int }_{-\infty }^{+\infty }{\rho }_{\sigma }\left(x\right)dx=\sqrt{2\pi }\sigma$.

Theorem 1.

For Algorithm 3, if $\sigma =k{\sigma }_2\ge {\eta }_{ϵ}\left(\mathbb{Z}\right)$, then the expected number of rejections that it requires is not more than ${\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)/{\sigma }_2\sqrt{\pi /2}\approx 1.47$.

Proof. 

From the perspective of rejection sampling, for Algorithm 3, the proposal distribution density function is $g\left(z\right)={\rho }_{{\sigma }_2}\left(x\right)/\left(2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)\right)$ and the target distribution density function f can be written as

$$f\left(z\right)=f\left(s(kx+y+(1+s)/2)\right)=\frac{{\rho }_{k{\sigma }_2,\mu }\left(s(kx+y+(1+s)/2)\right)}{{\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)},$$

where $s\in \{0,1\}$, $\mu =i/b$ and $i\in [0,b/2]\cap \mathbb{Z}$. The expected number of rejections that the algorithm requires can be evaluated by

$$\max \frac{f\left(z\right)}{g\left(z\right)}=\max \left(\frac{{\rho }_{k{\sigma }_2,\mu }\left(s(kx+y+(1+s)/2)\right)}{{\rho }_{{\sigma }_2}\left(x\right)}·\frac{2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{{\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)}\right)\le \frac{2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{{\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)}.$$

If $\sigma =k{\sigma }_2\ge {\eta }_{ϵ}\left(\mathbb{Z}\right)$, then ${\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)\approx {\rho }_{k{\sigma }_2}\left(\mathbb{Z}\right)$ (or ${\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)={\rho }_{k{\sigma }_2}\left(\mathbb{Z}\right)$ when $\mu =0$). Then we have

$$\frac{2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{{\rho }_{k{\sigma }_2,\mu }\left(\mathbb{Z}\right)}\approx \frac{2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{{\rho }_{k{\sigma }_2}\left(\mathbb{Z}\right)}\approx \frac{2k{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{2k{\sigma }_2\sqrt{\pi /2}}=\frac{{\rho }_{{\sigma }_2}\left({\mathbb{Z}}^{+}\right)}{{\sigma }_2\sqrt{\pi /2}}\approx 1.47,$$

where the two approximate equality symbols follow from Lemma 2. This completes the proof. □

We remark that Theorem 1 holds for $k\ge 1$ in practice. In the proof of Theorem 1, taking $ϵ=2^{-19}$ we have $(1-ϵ)/(1+ϵ)\approx 1$, which guarantees that ${\rho }_{\sigma ,\mu }\left(\mathbb{Z}\right)$ gets close enough to ${\rho }_{\sigma }\left(\mathbb{Z}\right)$ by Lemma 2. According to Lemma 1, we get ${\eta }_{ϵ}\left(\mathbb{Z}\right)\le \sqrt{\ln \left(2\right(1+1/ϵ\left)\right)/2}/\pi \approx 0.838036$ with $ϵ=2^{-19}$.

We can also see that the rejection rate of $1.47$ also represents the average number of times that Algorithm 3 calls Algorithm 2 or Algorithm 4. This is because the rejection rate of Algorithm 3 is actually equal to the inverse of the probability that the algorithm executes step 7 and step 10 without restarting. Algorithm 3 needs to perform step 7 whenever it successfully outputs the sample or restarts, and step 7 requires a call to Algorithm 2 or Algorithm 4. Thus, the average number of times that step 7 is executed is the rejection rate, and the average number of calls of Algorithm 2 or Algorithm 4 is thus the rejection rate, which is about $1.47$ and independent from the parameter $\sigma$ and the center $\mu$, according to Theorem 1.

5. Conclusions

In this paper, we show that the Bernoulli(-type) sampling for centered discrete Gaussian distributions over the integers can be implemented without floating-point arithmetic. Except a fixed look-up table of very small size (e.g., 128 bits) that stores the binary expansion of $\ln 2$, our proposed algorithm does not need any precomputation storage and requires no floating-point arithmetic. We also present a noncentered Bernoulli sampling algorithm that satisfies the need of discrete Gaussian distributions with varying centers. The algorithm needs no floating-point arithmetic, and it is suitable for sampling ${\mathcal{D}}_{\mathbb{Z},\sigma ,\mu }$ with $\mu \in [0,1)$ of precision up to 52 bits. The noncentered version of algorithm has the same rejection rate as the original Bernoulli sampling algorithm. The experiment shows that our improved algorithms have better performance as compared to other rejection algorithms for discrete Gaussian distributions over the integers.