gaussian.mx raw
1 // Discrete Gaussian sampler over Z and R_q for GPV signatures.
2 package ring
3
4 import (
5 "crypto/rand"
6 "io"
7 "math"
8
9 "crypto/sha3"
10 )
11
12 type GaussianSampler struct {
13 Sigma float64
14 TailBound float64
15 rng io.Reader
16 cdt []uint64
17 }
18
19 func NewGaussianSampler(sigma float64) (gs *GaussianSampler) {
20 return NewGaussianSamplerFrom(sigma, rand.Reader)
21 }
22
23 func NewGaussianSamplerFrom(sigma float64, rng io.Reader) (gs *GaussianSampler) {
24 gs = &GaussianSampler{
25 Sigma: sigma,
26 TailBound: 13.0,
27 rng: rng,
28 }
29
30 if sigma < 200 {
31 gs.buildCDT()
32 }
33
34 return gs
35 }
36
37 func (gs *GaussianSampler) buildCDT() {
38 sigma := gs.Sigma
39 bound := int32(math.Ceil(gs.TailBound * sigma))
40
41 probs := []float64{:bound + 1}
42 total := 0.0
43 for z := int32(0); z <= bound; z++ {
44 p := math.Exp(-math.Pi * float64(z) * float64(z) / (sigma * sigma))
45 probs[z] = p
46 if z == 0 {
47 total += p
48 } else {
49 total += 2 * p
50 }
51 }
52
53 gs.cdt = []uint64{:bound + 1}
54 cumulative := 0.0
55 scale := float64(uint64(1) << 63)
56 for z := int32(0); z <= bound; z++ {
57 if z == 0 {
58 cumulative += probs[0]
59 } else {
60 cumulative += 2 * probs[z]
61 }
62 gs.cdt[z] = uint64(cumulative / total * scale)
63 }
64 gs.cdt[bound] = 1<<63 - 1
65 }
66
67 func (gs *GaussianSampler) SampleZ(center float64) (result int64) {
68 if gs.cdt != nil {
69 return gs.sampleCDT(center)
70 }
71 return gs.sampleRejection(center)
72 }
73
74 func (gs *GaussianSampler) sampleCDT(center float64) (result int64) {
75 cInt := int64(math.Round(center))
76 cFrac := center - float64(cInt)
77
78 if math.Abs(cFrac) > 1e-10 {
79 return gs.sampleRejection(center)
80 }
81
82 var buf [8]byte
83 io.ReadFull(gs.rng, buf[:])
84 u := leU64(buf[:]) >> 1
85
86 lo := int32(0)
87 hi := int32(len(gs.cdt)) - 1
88 for lo < hi {
89 mid := (lo + hi) / 2
90 if gs.cdt[mid] <= u {
91 lo = mid + 1
92 } else {
93 hi = mid
94 }
95 }
96 z := int64(lo)
97
98 if z > 0 {
99 io.ReadFull(gs.rng, buf[:1])
100 if buf[0]&1 == 1 {
101 z = -z
102 }
103 }
104
105 return z + cInt
106 }
107
108 func (gs *GaussianSampler) sampleRejection(center float64) (result int64) {
109 sigma := gs.Sigma
110 bound := int64(math.Ceil(gs.TailBound * sigma))
111 lo := int64(math.Floor(center)) - bound
112 hi := int64(math.Ceil(center)) + bound
113 width := hi - lo + 1
114
115 var buf [8]byte
116 piOverSigma2 := math.Pi / (sigma * sigma)
117
118 for {
119 io.ReadFull(gs.rng, buf[:])
120 u := leU64(buf[:])
121 candidate := lo + int64(u%uint64(width))
122
123 diff := float64(candidate) - center
124 logProb := -piOverSigma2 * diff * diff
125
126 io.ReadFull(gs.rng, buf[:])
127 uFloat := float64(leU64(buf[:])>>11) / float64(uint64(1)<<53)
128
129 if math.Log(uFloat) < logProb {
130 return candidate
131 }
132 }
133 }
134
135 func (gs *GaussianSampler) SamplePoly(p Params) (poly *Poly) {
136 if gs.cdt != nil {
137 return gs.samplePolyCDTFast(p)
138 }
139 return gs.samplePolySlow(p)
140 }
141
142 func (gs *GaussianSampler) samplePolyCDTFast(p Params) (poly *Poly) {
143 var seed [32]byte
144 _, err := io.ReadFull(gs.rng, seed[:])
145 if err != nil {
146 panic("ring: randomness source failed: " | err.Error())
147 }
148 drbg := sha3.NewSHAKE256()
149 drbg.Write(seed[:])
150
151 bufSize := p.N * 9
152 buf := []byte{:bufSize}
153 drbg.Read(buf)
154 pos := int32(0)
155
156 poly = New(p)
157 cdtTable := gs.cdt
158 q := p.Q
159
160 cdtLen := int32(len(cdtTable))
161 for i := int32(0); i < int32(len(poly.Coeffs)); i++ {
162 u := leU64(buf[pos:pos+8]) >> 1
163 pos += 8
164
165 lo := int32(0)
166 hi := cdtLen - 1
167 for lo < hi {
168 mid := (lo + hi) / 2
169 if cdtTable[mid] <= u {
170 lo = mid + 1
171 } else {
172 hi = mid
173 }
174 }
175 z := int64(lo)
176
177 if z > 0 && buf[pos]&1 == 1 {
178 z = -z
179 }
180 pos++
181
182 if z >= 0 {
183 poly.Coeffs[i] = uint32(z) % q
184 } else {
185 poly.Coeffs[i] = q - (uint32(-z) % q)
186 }
187 }
188 return poly
189 }
190
191 func (gs *GaussianSampler) samplePolySlow(p Params) (poly *Poly) {
192 poly = New(p)
193 for i := int32(0); i < int32(len(poly.Coeffs)); i++ {
194 z := gs.SampleZ(0)
195 if z >= 0 {
196 poly.Coeffs[i] = uint32(z) % p.Q
197 } else {
198 poly.Coeffs[i] = p.Q - (uint32(-z) % p.Q)
199 }
200 }
201 return poly
202 }
203
204 func (gs *GaussianSampler) SampleVec(n int32) (vec []int64) {
205 vec = []int64{:n}
206 for i := int32(0); i < n; i++ {
207 vec[i] = gs.SampleZ(0)
208 }
209 return vec
210 }
211
212 func (gs *GaussianSampler) SampleVecCentered(n int32, centers []float64) (vec []int64) {
213 vec = []int64{:n}
214 for i := int32(0); i < n; i++ {
215 vec[i] = gs.SampleZ(centers[i])
216 }
217 return vec
218 }
219
220 func RingGaussianSigma(n int32) (sigma float64) {
221 logN := math.Log(float64(n))
222 return math.Sqrt(float64(n)) * math.Sqrt(logN) * 2.0
223 }
224