gaussian_test.go raw
1 package ring
2
3 import (
4 "math"
5 "testing"
6 )
7
8 // TestGaussianSampleZDistribution verifies that SampleZ produces
9 // a distribution that looks Gaussian: concentrated near 0 with the
10 // right spread.
11 func TestGaussianSampleZDistribution(t *testing.T) {
12 sigma := 5.0
13 gs := NewGaussianSampler(sigma)
14
15 n := 100000
16 samples := make([]int64, n)
17 for i := range samples {
18 samples[i] = gs.SampleZ(0)
19 }
20
21 // Compute mean and variance.
22 sum := 0.0
23 for _, s := range samples {
24 sum += float64(s)
25 }
26 mean := sum / float64(n)
27
28 sumSq := 0.0
29 for _, s := range samples {
30 d := float64(s) - mean
31 sumSq += d * d
32 }
33 variance := sumSq / float64(n)
34
35 // Expected variance for D_{Z, σ, 0} ≈ σ²/(2π).
36 expectedVar := sigma * sigma / (2 * math.Pi)
37
38 t.Logf("σ=%.1f: mean=%.3f, variance=%.3f (expected ≈ %.3f)",
39 sigma, mean, variance, expectedVar)
40
41 // Mean should be near 0.
42 if math.Abs(mean) > 0.5 {
43 t.Errorf("mean too far from 0: %.3f", mean)
44 }
45
46 // Variance should be within 20% of expected.
47 ratio := variance / expectedVar
48 if ratio < 0.7 || ratio > 1.3 {
49 t.Errorf("variance ratio %.3f outside [0.7, 1.3]", ratio)
50 }
51 }
52
53 // TestGaussianSampleZCentered verifies centered sampling.
54 func TestGaussianSampleZCentered(t *testing.T) {
55 sigma := 8.0
56 center := 42.0
57 gs := NewGaussianSampler(sigma)
58
59 n := 50000
60 sum := 0.0
61 for range n {
62 s := gs.SampleZ(center)
63 sum += float64(s)
64 }
65 mean := sum / float64(n)
66
67 t.Logf("σ=%.1f, c=%.1f: mean=%.3f", sigma, center, mean)
68
69 if math.Abs(mean-center) > 1.0 {
70 t.Errorf("mean %.3f too far from center %.1f", mean, center)
71 }
72 }
73
74 // TestGaussianSampleZLargeSigma verifies the rejection sampler
75 // for larger σ values used in GPV.
76 func TestGaussianSampleZLargeSigma(t *testing.T) {
77 sigma := 50.0
78 gs := NewGaussianSampler(sigma)
79
80 n := 50000
81 sum := 0.0
82 maxAbs := int64(0)
83 for range n {
84 s := gs.SampleZ(0)
85 sum += float64(s)
86 if s < 0 {
87 s = -s
88 }
89 if s > maxAbs {
90 maxAbs = s
91 }
92 }
93 mean := sum / float64(n)
94
95 t.Logf("σ=%.1f: mean=%.3f, max|sample|=%d (tail bound=%.0f)",
96 sigma, mean, maxAbs, 13.0*sigma)
97
98 // Max should be within tail bound.
99 if float64(maxAbs) > 13.0*sigma {
100 t.Errorf("sample %d exceeds tail bound %.0f", maxAbs, 13.0*sigma)
101 }
102
103 if math.Abs(mean) > 2.0 {
104 t.Errorf("mean %.3f too far from 0", mean)
105 }
106 }
107
108 // TestGaussianSamplePoly verifies polynomial sampling.
109 func TestGaussianSamplePoly(t *testing.T) {
110 sigma := 5.0
111 gs := NewGaussianSampler(sigma)
112 p := Falcon512()
113
114 poly := gs.SamplePoly(p)
115
116 // All coefficients should be in [0, Q).
117 for i, c := range poly.Coeffs {
118 if c >= p.Q {
119 t.Fatalf("coeff[%d] = %d >= Q=%d", i, c, p.Q)
120 }
121 }
122
123 // Count small coefficients (centered representation).
124 // Most should be near 0 (i.e., small or near Q).
125 smallCount := 0
126 half := p.Q / 2
127 for _, c := range poly.Coeffs {
128 var v uint32
129 if c > half {
130 v = p.Q - c
131 } else {
132 v = c
133 }
134 if v < 20 {
135 smallCount++
136 }
137 }
138
139 fraction := float64(smallCount) / float64(p.N)
140 t.Logf("σ=%.1f: %.1f%% of coefficients within ±20", sigma, fraction*100)
141
142 // With σ=5, most coefficients should be small.
143 if fraction < 0.5 {
144 t.Errorf("too few small coefficients: %.1f%%", fraction*100)
145 }
146 }
147
148 // TestGaussianCDTvRejection verifies CDT and rejection produce
149 // similar distributions.
150 func TestGaussianCDTvRejection(t *testing.T) {
151 sigma := 5.0
152
153 gsCDT := NewGaussianSampler(sigma)
154 if gsCDT.cdt == nil {
155 t.Fatal("CDT not built for σ=5")
156 }
157
158 // Also create a rejection sampler by using larger sigma threshold.
159 gsRej := &GaussianSampler{
160 Sigma: sigma,
161 TailBound: 13.0,
162 rng: gsCDT.rng,
163 }
164 // Force rejection path by not building CDT.
165
166 n := 50000
167
168 // Build histograms.
169 histCDT := make(map[int64]int)
170 histRej := make(map[int64]int)
171
172 for range n {
173 s := gsCDT.sampleCDT(0)
174 histCDT[s]++
175 }
176 for range n {
177 s := gsRej.sampleRejection(0)
178 histRej[s]++
179 }
180
181 // Compare distributions using chi-squared-like test.
182 // Just verify both are peaked at 0 and symmetric.
183 if histCDT[0] < n/10 {
184 t.Errorf("CDT: too few zeros: %d/%d", histCDT[0], n)
185 }
186 if histRej[0] < n/10 {
187 t.Errorf("Rejection: too few zeros: %d/%d", histRej[0], n)
188 }
189
190 t.Logf("CDT[0]=%d, CDT[±1]=%d/%d, CDT[±2]=%d/%d",
191 histCDT[0], histCDT[1], histCDT[-1], histCDT[2], histCDT[-2])
192 t.Logf("Rej[0]=%d, Rej[±1]=%d/%d, Rej[±2]=%d/%d",
193 histRej[0], histRej[1], histRej[-1], histRej[2], histRej[-2])
194 }
195
196 // TestRingGaussianSigma verifies the recommended sigma values.
197 func TestRingGaussianSigma(t *testing.T) {
198 for _, n := range []int{64, 256, 512, 1024} {
199 sigma := RingGaussianSigma(n)
200 t.Logf("n=%d: σ=%.1f", n, sigma)
201 if sigma <= 0 {
202 t.Errorf("n=%d: σ=%.1f should be positive", n, sigma)
203 }
204 }
205 }
206
207 func BenchmarkGaussianSampleZ(b *testing.B) {
208 gs := NewGaussianSampler(5.0)
209 b.ResetTimer()
210 for range b.N {
211 gs.SampleZ(0)
212 }
213 }
214
215 func BenchmarkGaussianSampleZLarge(b *testing.B) {
216 gs := NewGaussianSampler(165.0)
217 b.ResetTimer()
218 for range b.N {
219 gs.SampleZ(0)
220 }
221 }
222
223 func BenchmarkGaussianSamplePoly(b *testing.B) {
224 gs := NewGaussianSampler(165.0)
225 p := Falcon512()
226 b.ResetTimer()
227 for range b.N {
228 gs.SamplePoly(p)
229 }
230 }
231