gaussian.go raw
1 // Discrete Gaussian sampler over Z and R_q for GPV signatures.
2 //
3 // The discrete Gaussian distribution D_{Z, σ, c} over the integers assigns
4 // probability proportional to exp(-π(x-c)²/σ²) to each x ∈ Z.
5 //
6 // For lattice-based cryptography, we need:
7 // - D_{Z, σ}: 1-dimensional Gaussian over integers (base case)
8 // - D_{Z^n, σ}: n-dimensional by sampling each coordinate independently
9 // - D_{Λ, σ, c}: Gaussian over lattice coset (used in GPV signing)
10 //
11 // This implements the Peikert convolution technique and the
12 // rejection-sampling method of Ducas, Durmus, Lepoint, Lyubashevsky (2013).
13 //
14 // The parameter σ (sigma) controls the width of the distribution:
15 // - Larger σ → more spread → more secure but larger signatures
16 // - σ must exceed the smoothing parameter η_ε(Λ) for statistical closeness
17 // - For Falcon-512: σ ≈ 165 (at ~128-bit security)
18 //
19 // Security: the sampler must produce outputs statistically close to a true
20 // discrete Gaussian. Any bias leaks information about the trapdoor.
21
22 package ring
23
24 import (
25 "crypto/rand"
26 "encoding/binary"
27 "io"
28 "math"
29
30 "golang.org/x/crypto/sha3"
31 )
32
33 // GaussianSampler samples from the discrete Gaussian D_{Z, σ, c}.
34 type GaussianSampler struct {
35 // Sigma is the Gaussian parameter (standard deviation ≈ σ/√(2π)).
36 Sigma float64
37
38 // TailBound is the cutoff: |x - c| ≤ TailBound * σ.
39 // Beyond this, probability is negligible (< 2^{-128} for TailBound=13).
40 TailBound float64
41
42 // rng is the randomness source.
43 rng io.Reader
44
45 // Precomputed table for the CDT (Cumulative Distribution Table) method.
46 // Only used for small σ (σ < 10). For larger σ, use rejection sampling.
47 cdt []uint64
48 }
49
50 // NewGaussianSampler creates a sampler for D_{Z, σ, 0}.
51 func NewGaussianSampler(sigma float64) *GaussianSampler {
52 return NewGaussianSamplerFrom(sigma, rand.Reader)
53 }
54
55 // NewGaussianSamplerFrom creates a sampler with the given randomness source.
56 func NewGaussianSamplerFrom(sigma float64, rng io.Reader) *GaussianSampler {
57 gs := &GaussianSampler{
58 Sigma: sigma,
59 TailBound: 13.0, // 2^{-128} security
60 rng: rng,
61 }
62
63 // Build CDT for σ up to ~200 (Falcon-512 range).
64 // The table has ceil(13·σ) entries — 424 for σ=32.6, 2145 for σ=165.
65 // Binary search makes lookup O(log n) ≈ 9-11 comparisons regardless.
66 // Only fall back to rejection for σ > 200 where the table exceeds ~2600 entries.
67 if sigma < 200 {
68 gs.buildCDT()
69 }
70
71 return gs
72 }
73
74 // buildCDT constructs the cumulative distribution table.
75 // For each z ≥ 0, CDF[z] = Σ_{x=0}^{z} ρ_σ(x), scaled to [0, 2^63).
76 func (gs *GaussianSampler) buildCDT() {
77 sigma := gs.Sigma
78 bound := int(math.Ceil(gs.TailBound * sigma))
79
80 // Compute unnormalized probabilities for z = 0, 1, ..., bound.
81 probs := make([]float64, bound+1)
82 total := 0.0
83 for z := 0; z <= bound; z++ {
84 p := math.Exp(-math.Pi * float64(z) * float64(z) / (sigma * sigma))
85 probs[z] = p
86 if z == 0 {
87 total += p
88 } else {
89 total += 2 * p // both +z and -z
90 }
91 }
92
93 // Build CDF scaled to 2^63.
94 gs.cdt = make([]uint64, bound+1)
95 cumulative := 0.0
96 scale := float64(uint64(1) << 63)
97 for z := 0; z <= bound; z++ {
98 // P(|X| ≤ z) = sum of probs for 0..z, where z>0 counts twice
99 if z == 0 {
100 cumulative += probs[0]
101 } else {
102 cumulative += 2 * probs[z]
103 }
104 gs.cdt[z] = uint64(cumulative / total * scale)
105 }
106 // Ensure last entry is max.
107 gs.cdt[bound] = 1<<63 - 1
108 }
109
110 // SampleZ samples from D_{Z, σ, c} (discrete Gaussian centered at c).
111 func (gs *GaussianSampler) SampleZ(center float64) int64 {
112 if gs.cdt != nil {
113 return gs.sampleCDT(center)
114 }
115 return gs.sampleRejection(center)
116 }
117
118 // sampleCDT uses the cumulative distribution table for small σ.
119 // Samples |z| from the table, then randomly assigns sign.
120 func (gs *GaussianSampler) sampleCDT(center float64) int64 {
121 // Sample from D_{Z, σ, 0} then shift by round(center).
122 cInt := int64(math.Round(center))
123 cFrac := center - float64(cInt)
124
125 // If fractional part is significant, fall back to rejection.
126 if math.Abs(cFrac) > 1e-10 {
127 return gs.sampleRejection(center)
128 }
129
130 // Sample |z| using CDT with binary search.
131 var buf [8]byte
132 io.ReadFull(gs.rng, buf[:])
133 u := binary.LittleEndian.Uint64(buf[:]) >> 1 // 63-bit uniform
134
135 // Binary search: find smallest z where cdt[z] > u.
136 lo, hi := 0, len(gs.cdt)-1
137 for lo < hi {
138 mid := (lo + hi) / 2
139 if gs.cdt[mid] <= u {
140 lo = mid + 1
141 } else {
142 hi = mid
143 }
144 }
145 z := int64(lo)
146
147 // Random sign for z > 0.
148 if z > 0 {
149 io.ReadFull(gs.rng, buf[:1])
150 if buf[0]&1 == 1 {
151 z = -z
152 }
153 }
154
155 return z + cInt
156 }
157
158 // sampleRejection uses rejection sampling for arbitrary σ and center.
159 // Samples a candidate from a uniform distribution over the tail bound,
160 // then accepts with probability proportional to ρ_σ(x - c).
161 func (gs *GaussianSampler) sampleRejection(center float64) int64 {
162 sigma := gs.Sigma
163 bound := int64(math.Ceil(gs.TailBound * sigma))
164 lo := int64(math.Floor(center)) - bound
165 hi := int64(math.Ceil(center)) + bound
166 width := hi - lo + 1
167
168 var buf [8]byte
169 piOverSigma2 := math.Pi / (sigma * sigma)
170
171 for {
172 io.ReadFull(gs.rng, buf[:])
173 // Sample candidate uniformly from [lo, hi].
174 u := binary.LittleEndian.Uint64(buf[:])
175 candidate := lo + int64(u%uint64(width))
176
177 // Accept with probability exp(-π(candidate - center)² / σ²).
178 diff := float64(candidate) - center
179 logProb := -piOverSigma2 * diff * diff
180
181 // Sample a uniform [0, 1) and check.
182 io.ReadFull(gs.rng, buf[:])
183 uFloat := float64(binary.LittleEndian.Uint64(buf[:])>>11) / float64(uint64(1)<<53)
184
185 if math.Log(uFloat) < logProb {
186 return candidate
187 }
188 }
189 }
190
191 // SamplePoly samples a polynomial with each coefficient drawn from D_{Z, σ, 0}.
192 //
193 // Uses a SHAKE256 DRBG seeded from one 32-byte crypto/rand read, with all
194 // randomness squeezed into a single buffer upfront. This reduces kernel
195 // crossings from 128+ to 1 and eliminates per-sample interface dispatch.
196 // SHAKE256 output is computationally indistinguishable from random
197 // (it's the PRF used in ML-KEM/Kyber).
198 func (gs *GaussianSampler) SamplePoly(p Params) *Poly {
199 if gs.cdt != nil {
200 return gs.samplePolyCDTFast(p)
201 }
202 return gs.samplePolySlow(p)
203 }
204
205 // samplePolyCDTFast is the optimized CDT path with bulk randomness.
206 func (gs *GaussianSampler) samplePolyCDTFast(p Params) *Poly {
207 // Seed a SHAKE256 XOF from one crypto/rand read.
208 var seed [32]byte
209 if _, err := io.ReadFull(gs.rng, seed[:]); err != nil {
210 panic("ring: randomness source failed: " + err.Error())
211 }
212 drbg := sha3.NewShake256()
213 drbg.Write(seed[:])
214
215 // Pre-squeeze all randomness: 9 bytes per coefficient (8 for CDT + 1 for sign).
216 // Over-allocate slightly for rejection fallback (rare with CDT).
217 bufSize := p.N * 9
218 buf := make([]byte, bufSize)
219 drbg.Read(buf)
220 pos := 0
221
222 poly := New(p)
223 cdtTable := gs.cdt
224 q := p.Q
225
226 cdtLen := len(cdtTable)
227 for i := range poly.Coeffs {
228 // Inline CDT sampling with binary search.
229 u := binary.LittleEndian.Uint64(buf[pos:pos+8]) >> 1 // 63-bit uniform
230 pos += 8
231
232 // Binary search: find smallest z where cdtTable[z] > u.
233 lo, hi := 0, cdtLen-1
234 for lo < hi {
235 mid := (lo + hi) / 2
236 if cdtTable[mid] <= u {
237 lo = mid + 1
238 } else {
239 hi = mid
240 }
241 }
242 z := int64(lo)
243
244 // Random sign for z > 0.
245 if z > 0 && buf[pos]&1 == 1 {
246 z = -z
247 }
248 pos++
249
250 // Convert to centered representation mod Q.
251 if z >= 0 {
252 poly.Coeffs[i] = uint32(z) % q
253 } else {
254 poly.Coeffs[i] = q - (uint32(-z) % q)
255 }
256 }
257 return poly
258 }
259
260 // samplePolySlow is the fallback for large σ (rejection sampling).
261 func (gs *GaussianSampler) samplePolySlow(p Params) *Poly {
262 poly := New(p)
263 for i := range poly.Coeffs {
264 z := gs.SampleZ(0)
265 if z >= 0 {
266 poly.Coeffs[i] = uint32(z) % p.Q
267 } else {
268 poly.Coeffs[i] = p.Q - (uint32(-z) % p.Q)
269 }
270 }
271 return poly
272 }
273
274 // SampleVec samples a vector of n independent Gaussian integers.
275 func (gs *GaussianSampler) SampleVec(n int) []int64 {
276 vec := make([]int64, n)
277 for i := range vec {
278 vec[i] = gs.SampleZ(0)
279 }
280 return vec
281 }
282
283 // SampleVecCentered samples a vector of n Gaussian integers centered at c.
284 func (gs *GaussianSampler) SampleVecCentered(n int, centers []float64) []int64 {
285 vec := make([]int64, n)
286 for i := range vec {
287 vec[i] = gs.SampleZ(centers[i])
288 }
289 return vec
290 }
291
292 // RingGaussianSigma returns the recommended σ for a given ring dimension
293 // to achieve security level λ bits.
294 //
295 // For GPV signatures over R_q with n coefficients:
296 // σ ≈ s1(R) · ω(√log n)
297 // where s1(R) is the largest singular value of the trapdoor basis.
298 //
299 // Practical values:
300 // n=512: σ ≈ 165 (Falcon-512 target)
301 // n=1024: σ ≈ 168 (Falcon-1024 target)
302 // n=64: σ ≈ 50 (for testing)
303 func RingGaussianSigma(n int) float64 {
304 // σ = √(n) * √(log(n)) * constant
305 // The constant is tuned so that σ/q gives appropriate noise ratio.
306 logN := math.Log(float64(n))
307 return math.Sqrt(float64(n)) * math.Sqrt(logN) * 2.0
308 }
309