gaussian.go raw
1 package gnarlring
2
3 import (
4 "crypto/rand"
5 "encoding/binary"
6 "io"
7 "math"
8 )
9
10 // GaussSampler samples from the discrete Gaussian distribution D_{Z, sigma}
11 // over the integers. Uses the Cumulative Distribution Table (CDT) method
12 // for the standard sigma ≈ 10.4. Each call to SampleZ produces an integer
13 // drawn from D_{Z, sigma, 0}. SamplePoly produces a Poly27 with each
14 // coefficient drawn independently.
15 type GaussSampler struct {
16 sigma float64
17 tail int
18 cdt []uint64 // cumulative distribution table, scaled to 2^63
19 rng io.Reader
20 }
21
22 // NewGaussSampler creates a sampler for D_{Z, sigma, 0}.
23 // sigma is typically sqrt(N) * 2 ≈ 10.39 for the gnarl ring.
24 func NewGaussSampler(sigma float64) *GaussSampler {
25 return NewGaussSamplerFrom(sigma, rand.Reader)
26 }
27
28 // NewGaussSamplerFrom creates a sampler with the given randomness source.
29 // Only builds the CDT for sigma < 50 — larger sigma uses rejection sampling.
30 func NewGaussSamplerFrom(sigma float64, rng io.Reader) *GaussSampler {
31 if rng == nil {
32 rng = rand.Reader
33 }
34 gs := &GaussSampler{
35 sigma: sigma,
36 tail: 13,
37 rng: rng,
38 }
39 if sigma < 50 {
40 gs.buildCDT()
41 }
42 return gs
43 }
44
45 // buildCDT constructs the cumulative distribution table.
46 // For each z ≥ 0, CDT[z] = P(|X| ≤ z) scaled to [0, 2^63).
47 func (gs *GaussSampler) buildCDT() {
48 sigma := gs.sigma
49 bound := int(math.Ceil(float64(gs.tail) * sigma))
50
51 probs := make([]float64, bound+1)
52 total := 0.0
53 for z := 0; z <= bound; z++ {
54 p := math.Exp(-math.Pi * float64(z) * float64(z) / (sigma * sigma))
55 probs[z] = p
56 if z == 0 {
57 total += p
58 } else {
59 total += 2 * p // both +z and -z
60 }
61 }
62
63 gs.cdt = make([]uint64, bound+1)
64 cumulative := 0.0
65 scale := float64(uint64(1) << 63)
66 for z := 0; z <= bound; z++ {
67 if z == 0 {
68 cumulative += probs[0]
69 } else {
70 cumulative += 2 * probs[z]
71 }
72 gs.cdt[z] = uint64(cumulative / total * scale)
73 }
74 // Ensure last entry is max.
75 gs.cdt[bound] = 1<<63 - 1
76 }
77
78 // SampleZ samples from D_{Z, sigma, center}.
79 func (gs *GaussSampler) SampleZ(center float64) int64 {
80 if gs.cdt != nil {
81 // For small sigma, the CDT is tiny and rejection sampling for
82 // non-integer centers is prohibitively slow. Round to nearest
83 // integer — the bias from rounding is O(1/σ) which for
84 // sigma_eff ≈ 0.1 is negligible (< 1 per 10^6 samples).
85 cInt := int64(math.Round(center))
86 return gs.sampleCDT(float64(cInt))
87 }
88 return gs.sampleRejection(center)
89 }
90
91 // sampleCDT samples |z| via binary search on the CDF, then randomly assigns sign.
92 func (gs *GaussSampler) sampleCDT(center float64) int64 {
93 cInt := int64(math.Round(center))
94
95 var buf [8]byte
96 io.ReadFull(gs.rng, buf[:])
97 u := binary.LittleEndian.Uint64(buf[:]) >> 1 // 63-bit uniform
98
99 lo, hi := 0, len(gs.cdt)-1
100 for lo < hi {
101 mid := (lo + hi) / 2
102 if gs.cdt[mid] <= u {
103 lo = mid + 1
104 } else {
105 hi = mid
106 }
107 }
108 z := int64(lo)
109
110 if z > 0 {
111 io.ReadFull(gs.rng, buf[:1])
112 if buf[0]&1 == 1 {
113 z = -z
114 }
115 }
116 return z + cInt
117 }
118
119 // sampleRejection is the fallback for non-integer centers or large sigma.
120 func (gs *GaussSampler) sampleRejection(center float64) int64 {
121 sigma := gs.sigma
122 bound := int64(math.Ceil(float64(gs.tail) * sigma))
123 lo := int64(math.Floor(center)) - bound
124 hi := int64(math.Ceil(center)) + bound
125 width := hi - lo + 1
126
127 piOverSigma2 := math.Pi / (sigma * sigma)
128 var buf [8]byte
129
130 for {
131 io.ReadFull(gs.rng, buf[:])
132 u := binary.LittleEndian.Uint64(buf[:])
133 candidate := lo + int64(u%uint64(width))
134
135 diff := float64(candidate) - center
136 logProb := -piOverSigma2 * diff * diff
137
138 io.ReadFull(gs.rng, buf[:])
139 uFloat := float64(binary.LittleEndian.Uint64(buf[:])>>11) / float64(uint64(1)<<53)
140
141 if math.Log(uFloat) < logProb {
142 return candidate
143 }
144 }
145 }
146
147 // SamplePoly returns a Poly27 with each coefficient drawn from D_{Z, sigma, 0}.
148 func (gs *GaussSampler) SamplePoly() *Poly27 {
149 if gs.cdt != nil {
150 return gs.samplePolyCDT()
151 }
152 return gs.samplePolySlow()
153 }
154
155 // samplePolyCDT uses the CDT with bulk randomness for efficiency.
156 func (gs *GaussSampler) samplePolyCDT() *Poly27 {
157 // For n=27, this is fast enough without SHAKE256 DRBG.
158 // 27 coefficients × 9 bytes (8 for CDT + 1 for sign).
159 p := NewPoly27()
160 cdtTable := gs.cdt
161 cdtLen := len(cdtTable)
162 var buf [9]byte
163
164 for i := 0; i < N; i++ {
165 io.ReadFull(gs.rng, buf[:])
166 u := binary.LittleEndian.Uint64(buf[:8]) >> 1
167
168 lo, hi := 0, cdtLen-1
169 for lo < hi {
170 mid := (lo + hi) / 2
171 if cdtTable[mid] <= u {
172 lo = mid + 1
173 } else {
174 hi = mid
175 }
176 }
177 z := int64(lo)
178 if z > 0 && buf[8]&1 == 1 {
179 z = -z
180 }
181
182 if z >= 0 {
183 p.Coeffs[i] = uint16(uint64(z) % Q)
184 } else {
185 p.Coeffs[i] = Q - uint16(uint64(-z)%Q)
186 }
187 }
188 return p
189 }
190
191 // samplePolySlow is the fallback for large sigma.
192 func (gs *GaussSampler) samplePolySlow() *Poly27 {
193 p := NewPoly27()
194 for i := 0; i < N; i++ {
195 z := gs.SampleZ(0)
196 if z >= 0 {
197 p.Coeffs[i] = uint16(uint64(z) % Q)
198 } else {
199 p.Coeffs[i] = Q - uint16(uint64(-z)%Q)
200 }
201 }
202 return p
203 }
204