falcon.go raw
1 // Falcon-like NTRU lattice signatures.
2 //
3 // Uses the NTRU-structured trapdoor to produce compact signatures (~700 B for
4 // Falcon-512). The verification key is a single polynomial h = g/f mod q, and
5 // the signature is a single short polynomial s2. The verifier derives s1 from
6 // the verification equation and checks both are short.
7 //
8 // Key generation:
9 // 1. Sample short f, g ← D_{Z^n, σ}
10 // 2. Compute h = g·f^{-1} mod q
11 // 3. Find short F, G such that f·G - g·F = q (NTRU equation)
12 // 4. Secret key: (f, g, F, G). Public key: h.
13 //
14 // Signing (ffSampling):
15 // 1. c = HashToPoint(m) ∈ R_q (sparse ±1 challenge)
16 // 2. Use the NTRU basis B = [(g, -f), (G, -F)] with Gram-Schmidt
17 // orthogonalization to sample a Gaussian preimage (s1, s2) such that
18 // s1 + h·s2 = c (mod q).
19 // 3. Output s2 as the signature (s1 is derived during verification).
20 //
21 // Verification:
22 // 1. c = HashToPoint(m)
23 // 2. s1 = c - h·s2 (mod q)
24 // 3. Check ||s1||_∞, ||s2||_∞ ≤ bound
25
26 package ring
27
28 import (
29 "crypto/rand"
30 "io"
31 "math"
32 )
33
34 // FalconPublicKey is the NTRU verification key.
35 type FalconPublicKey struct {
36 H *Poly // h = g·f^{-1} mod q (coefficient form)
37 P Params
38 }
39
40 // FalconSecretKey is the NTRU trapdoor basis.
41 type FalconSecretKey struct {
42 F, G, f, g *Poly // short NTRU basis (coefficient form)
43 Sigma float64
44 PK *FalconPublicKey
45 }
46
47 // FalconSignature is a compact NTRU signature (single short polynomial).
48 type FalconSignature struct {
49 S2 *Poly // signature component (short, compressed)
50 }
51
52 // DefaultFalconParams returns Falcon-512 parameters.
53 func DefaultFalconParams() Params {
54 return Falcon512()
55 }
56
57 // ─── Polynomial dot product and Gram-Schmidt ──────────────────────────────
58
59 // polyDot computes the Euclidean dot product Σ a_i·b_i (scalar).
60 // Both polynomials are in coefficient form. Used for Gram-Schmidt
61 // orthogonalization in the ffSampling algorithm.
62 func polyDot(a, b *Poly) int64 {
63 n := len(a.Coeffs)
64 var sum int64
65 half := int64(a.params.Q / 2)
66 for i := 0; i < n; i++ {
67 ai := int64(a.Coeffs[i])
68 bi := int64(b.Coeffs[i])
69 if ai > half {
70 ai -= int64(a.params.Q)
71 }
72 if bi > half {
73 bi -= int64(a.params.Q)
74 }
75 sum += ai * bi
76 }
77 return sum
78 }
79
80 // polyVectorDot computes the Euclidean dot product of two 2-vectors of
81 // polynomials: <(a1,a2), (b1,b2)> = dot(a1,b1) + dot(a2,b2).
82 func polyVectorDot(a1, a2, b1, b2 *Poly) int64 {
83 return polyDot(a1, b1) + polyDot(a2, b2)
84 }
85
86 // gramSchmidtNormSq returns ||b*||² for the Gram-Schmidt orthogonalization
87 // of the second basis vector. For basis B = [b1, b2] where b1 and b2 are
88 // 2-vectors of polynomials, computes:
89 //
90 // µ = <b2, b1> / <b1, b1>
91 // b2* = b2 - µ·b1
92 // return <b2*, b2*>
93 func gramSchmidtNormSq(b11, b12, b21, b22 *Poly) (normSq float64) {
94 b1NormSq := float64(polyVectorDot(b11, b12, b11, b12))
95 b2b1 := float64(polyVectorDot(b21, b22, b11, b12))
96 µ := b2b1 / b1NormSq
97
98 // b2*_i = b2_i - µ·b1_i for each coefficient of each polynomial
99 n := len(b11.Coeffs)
100 var sum float64
101 for i := 0; i < n; i++ {
102 b21c := float64(b21.Coeffs[i])
103 b11c := float64(b11.Coeffs[i])
104 b22c := float64(b22.Coeffs[i])
105 b12c := float64(b12.Coeffs[i])
106
107 diff1 := b21c - µ*b11c
108 diff2 := b22c - µ*b12c
109 sum += diff1*diff1 + diff2*diff2
110 }
111 return sum
112 }
113
114 // ─── NTRU Key Generation ─────────────────────────────────────────────────
115
116 // polyInverse computes f^{-1} mod (x^n+1, q) using NTT-based pointwise inverse.
117 // f is in coefficient form. Returns nil if f is not invertible.
118 // For Falcon-512 (q=12289, q≡1 mod 2n), the NTT splits the ring completely,
119 // so pointwise inversion in the NTT domain is correct.
120 func polyInverse(f *Poly) *Poly {
121 p := f.params
122 q := p.Q
123
124 fNTT := f.Clone()
125 NTT(fNTT)
126
127 resultNTT := New(p)
128 for i, c := range fNTT.Coeffs {
129 if c == 0 {
130 return nil // not invertible — has a zero in NTT domain
131 }
132 resultNTT.Coeffs[i] = powMod(c, q-2, q)
133 }
134 resultNTT.isNTT = true
135
136 INTT(resultNTT)
137 return resultNTT
138 }
139
140 // computeNTRUBasis finds short F, G such that f·G - g·F = q (mod x^n+1).
141 // Uses the XGCD in Z_q[x]/(x^n+1) followed by Babai reduction.
142 func computeNTRUBasis(f, g *Poly, sigma float64, rng io.Reader) (*Poly, *Poly) {
143 p := f.params
144
145 // G0 = q·f^{-1}, F0 = 0 satisfies f·G0 - g·F0 = q.
146 fInv := polyInverse(f)
147 if fInv == nil {
148 return New(p), New(p)
149 }
150 G0 := PolyScalarMul(fInv, uint32(p.Q))
151 F0 := New(p) // zero polynomial
152
153 // Babai reduction: iteratively subtract k·(f, g) from (F, G) to make them short.
154 // k = round(<(F,G), (f,g)> / ||(f,g)||²)
155 b1Norm := float64(polyVectorDot(f, g, f, g))
156 if b1Norm < 1 {
157 b1Norm = 1
158 }
159
160 Fcur, Gcur := F0, G0
161 for iter := 0; iter < 10; iter++ {
162 b2b1 := float64(polyVectorDot(Fcur, Gcur, f, g))
163 kFloat := b2b1 / b1Norm
164 k := int(math.Round(kFloat))
165 if k == 0 {
166 break
167 }
168 // Subtract k·(f, g) from (F, G).
169 for i := range Fcur.Coeffs {
170 Fcur.Coeffs[i] = subMod(Fcur.Coeffs[i], mulMod(uint32(k), f.Coeffs[i], p.Q), p.Q)
171 Gcur.Coeffs[i] = subMod(Gcur.Coeffs[i], mulMod(uint32(k), g.Coeffs[i], p.Q), p.Q)
172 }
173 }
174
175 return Fcur, Gcur
176 }
177
178 // PolyScalarMul multiplies each coefficient of a polynomial by a scalar mod q.
179 func PolyScalarMul(a *Poly, s uint32) *Poly {
180 c := a.Clone()
181 q := a.params.Q
182 for i := range c.Coeffs {
183 c.Coeffs[i] = mulMod(c.Coeffs[i], s, q)
184 }
185 c.isNTT = false
186 return c
187 }
188
189 // FalconKeyGen generates a Falcon NTRU key pair.
190 func FalconKeyGen(p Params) (*FalconPublicKey, *FalconSecretKey) {
191 return FalconKeyGenFrom(p, rand.Reader)
192 }
193
194 // FalconKeyGenFrom generates a Falcon NTRU key pair from the given RNG.
195 func FalconKeyGenFrom(p Params, rng io.Reader) (*FalconPublicKey, *FalconSecretKey) {
196 sigma := math.Sqrt(float64(p.N)) * 2.0
197
198 gs := NewGaussianSamplerFrom(sigma, rng)
199
200 // Sample f, g from discrete Gaussian. Resample f if not invertible
201 // (any NTT slot is zero — coefficient 0 check alone is insufficient).
202 var f, g, fInv *Poly
203 for {
204 f = gs.SamplePoly(p)
205 g = gs.SamplePoly(p)
206 fInv = polyInverse(f)
207 if fInv != nil {
208 break
209 }
210 }
211
212 // Compute h = g·f^{-1} mod q.
213 hCoeff := Mul(fInv, g) // coefficient form
214
215 // Find short F, G.
216 F, G := computeNTRUBasis(f, g, sigma, rng)
217
218 pk := &FalconPublicKey{H: hCoeff, P: p}
219 sk := &FalconSecretKey{F: F, G: G, f: f, g: g, Sigma: sigma, PK: pk}
220 return pk, sk
221 }
222
223 // ─── ffSampling (Gaussian preimage sampling) ─────────────────────────────
224
225 // ffSampling samples a short vector (s1, s2) such that s1 + h·s2 = target
226 // using the NTRU trapdoor basis with the ffSampling algorithm.
227 //
228 // The algorithm:
229 // 1. Set T = (0, target) — a non-short preimage
230 // 2. For each basis vector b_i (from bottom to top):
231 // c_i = <T, b_i*> / ||b_i*||²
232 // z_i ← D_{Z, σ/||b_i*||, c_i}
233 // T = T - z_i·b_i
234 // 3. Return T as the short preimage (s1, s2)
235 func ffSampling(p Params, f, g, F, G *Poly, target *Poly, sigma float64, rng io.Reader) (*Poly, *Poly) {
236 // NTRU basis vectors: b1 = (g, -f), b2 = (F, -G)
237 // Wait — verify the equation: s1 + h·s2 = target
238 // (g, -f): g + h·(-f) = g - h·f = g - g = 0 ✓
239 // (F, -G): F + h·(-G) = F - h·G ≠ 0
240 //
241 // Actually, the second basis vector is scaled: we need (G, -F)
242 // Check: G + h·(-F) = G - h·F. From the NTRU equation f·G - g·F = q,
243 // we have -g·F = q - f·G, so F = (f·G - q)/g. Then G - h·F = G - g/f·F = 0
244 // when multiplied by f: f·G - g·F = q. So f·(G - h·F) = q, not 0.
245 //
246 // For the ffSampling, both (g, -f) and ( -F, G) form the basis.
247 // The lattice generated by these two vectors has determinant q.
248
249 b1 := [2]*Poly{g, Neg(f)} // (g, -f)
250 b2 := [2]*Poly{Neg(G), F} // (-G, F) — columns of the NTRU basis
251
252 // Current target T = (0, target).
253 t1 := New(p) // zero polynomial
254 t2 := target.Clone()
255
256 // Gram-Schmidt norms.
257 b2NormSq := gramSchmidtNormSq(b1[0], b1[1], b2[0], b2[1])
258
259 // === Process b2 (second basis vector) ===
260
261 // c2 = <T, b2*> / ||b2*||²
262 // Compute the inner product <T, (b2 - µ·b1)> / ||b2*||²
263 b1NormSq := float64(polyVectorDot(b1[0], b1[1], b1[0], b1[1]))
264 if b1NormSq < 1 {
265 b1NormSq = 1
266 }
267 b2b1 := float64(polyVectorDot(b2[0], b2[1], b1[0], b1[1]))
268 µ := b2b1 / b1NormSq
269
270 // <T, b2*> = <T, b2> - µ·<T, b1>
271 tb2 := float64(polyVectorDot(t1, t2, b2[0], b2[1]))
272 tb1 := float64(polyVectorDot(t1, t2, b1[0], b1[1]))
273 c2 := (tb2 - µ*tb1) / b2NormSq
274
275 // Sample z2 from D_{Z, σ/||b2*||, c2}.
276 sigma2 := sigma / math.Sqrt(b2NormSq)
277 z2 := sampleZ(sigma2, c2, rng)
278
279 // T = T - z2·b2 (subtract z2 copies of the second basis vector).
280 // Subtracting z2 from the coefficient along b2 brings us closer to the lattice.
281 for i := 0; i < len(t1.Coeffs); i++ {
282 t1.Coeffs[i] = subMod(t1.Coeffs[i], mulMod(uint32(z2), b2[0].Coeffs[i], p.Q), p.Q)
283 t2.Coeffs[i] = subMod(t2.Coeffs[i], mulMod(uint32(z2), b2[1].Coeffs[i], p.Q), p.Q)
284 }
285
286 // === Process b1 (first basis vector) ===
287
288 // c1 = <T, b1*> / ||b1*||² = <T, b1> / ||b1||²
289 tb1_2 := float64(polyVectorDot(t1, t2, b1[0], b1[1]))
290 c1 := tb1_2 / b1NormSq
291
292 sigma1 := sigma / math.Sqrt(b1NormSq)
293 z1 := sampleZ(sigma1, c1, rng)
294
295 // T = T - z1·b1
296 for i := 0; i < len(t1.Coeffs); i++ {
297 t1.Coeffs[i] = subMod(t1.Coeffs[i], mulMod(uint32(z1), b1[0].Coeffs[i], p.Q), p.Q)
298 t2.Coeffs[i] = subMod(t2.Coeffs[i], mulMod(uint32(z1), b1[1].Coeffs[i], p.Q), p.Q)
299 }
300
301 return t1, t2
302 }
303
304 // sampleZ samples a single integer from D_{Z, σ, c} using rejection sampling.
305 // For sigma < 1, the distribution is concentrated at the nearest integer
306 // (probability of any other value is < exp(-π/1) ≈ 0.04 and decreasing
307 // rapidly). Returns the rounded center directly to avoid rejection loop hangs.
308 func sampleZ(sigma, center float64, rng io.Reader) int64 {
309 if sigma < 1.0 || math.IsNaN(sigma) {
310 return int64(math.Round(center))
311 }
312
313 tail := int(math.Ceil(13.0 * sigma))
314 if tail < 1 {
315 tail = 1
316 }
317 centerInt := int64(math.Round(center))
318
319 var buf [8]byte
320 piOverSigma2 := math.Pi / (sigma * sigma)
321
322 for {
323 io.ReadFull(rng, buf[:])
324 u := int64(readLEUint64(buf[:]))
325 candidate := centerInt + (u % int64(2*tail+1)) - int64(tail)
326
327 diff := float64(candidate) - center
328 logProb := -piOverSigma2 * diff * diff
329
330 io.ReadFull(rng, buf[:])
331 uFloat := float64(readLEUint64(buf[:])>>11) / float64(uint64(1)<<53)
332
333 if math.Log(uFloat) < logProb {
334 return candidate
335 }
336 }
337 }
338
339 // readLEUint64 reads a little-endian uint64 from a byte slice.
340 func readLEUint64(b []byte) uint64 {
341 _ = b[7]
342 return uint64(b[0]) | uint64(b[1])<<8 | uint64(b[2])<<16 | uint64(b[3])<<24 |
343 uint64(b[4])<<32 | uint64(b[5])<<40 | uint64(b[6])<<48 | uint64(b[7])<<56
344 }
345
346 // ─── Falcon signing and verification ─────────────────────────────────────
347
348 // hashToFalconTarget hashes a message to a short polynomial c for Falcon.
349 // Uses domain separation from the GPV and sigma functions.
350 func hashToFalconTarget(p Params, message []byte) *Poly {
351 return hashMessageToTarget(p, message) // same format: τ=40 (or τ=n/2) ±1 coeffs
352 }
353
354 // FalconSign signs a message using the NTRU secret key. Produces a single
355 // short polynomial s2 as the signature. The verifier derives s1 from the
356 // verification equation and checks both are short.
357 func FalconSign(sk *FalconSecretKey, message []byte) *FalconSignature {
358 return FalconSignFrom(sk, message, rand.Reader)
359 }
360
361 // FalconSignFrom signs with the given randomness source.
362 func FalconSignFrom(sk *FalconSecretKey, message []byte, rng io.Reader) *FalconSignature {
363 p := sk.PK.P
364 sigma := sk.Sigma
365
366 // c = HashToPoint(m)
367 c := hashToFalconTarget(p, message)
368
369 // ffSampling: find short (s1, s2) such that s1 + h·s2 = c.
370 s1, s2 := ffSampling(p, sk.f, sk.g, sk.F, sk.G, c, sigma, rng)
371 _ = s1
372
373 return &FalconSignature{S2: s2}
374 }
375
376 // FalconVerify verifies a Falcon signature.
377 //
378 // Checks:
379 // 1. s1 = c - h·s2 (mod q)
380 // 2. ||s1||_∞, ||s2||_∞ ≤ σ·1.5 + 13·σ (= bound for Gaussian preimage)
381 func FalconVerify(pk *FalconPublicKey, message []byte, sig *FalconSignature) bool {
382 p := pk.P
383 sigma := math.Sqrt(float64(p.N)) * 2.0
384
385 s2 := sig.S2
386
387 // c = HashToPoint(m).
388 c := hashToFalconTarget(p, message)
389
390 // s1 = c - h·s2 (mod q).
391 hs2 := Mul(pk.H, s2)
392 s1 := Sub(c, hs2)
393
394 // Norm check: ||s1||_∞, ||s2||_∞ ≤ bound.
395 // The bound accounts for the Gaussian preimage sampling width.
396 bound := uint32(sigma*1.5 + 13*sigma)
397 if Norm(s1) > bound {
398 return false
399 }
400 if Norm(s2) > bound {
401 return false
402 }
403
404 return true
405 }
406
407 // FalconSigBound returns the norm bound for Falcon signatures.
408 func FalconSigBound(p Params) uint32 {
409 sigma := math.Sqrt(float64(p.N)) * 2.0
410 return uint32(sigma*1.5 + 13*sigma)
411 }
412
413 // FalconSignFrom returns (s1, s2) pair for verification debugging.
414 func FalconSignWithS1(sk *FalconSecretKey, message []byte, rng io.Reader) (s1, s2 *Poly) {
415 p := sk.PK.P
416 sigma := sk.Sigma
417 c := hashToFalconTarget(p, message)
418 s1, s2 = ffSampling(p, sk.f, sk.g, sk.F, sk.G, c, sigma, rng)
419 return s1, s2
420 }
421
422 // Verify the s1 is derived correctly.
423 func verifyFalconEquation(pk *FalconPublicKey, c, s1, s2 *Poly) bool {
424 // c = s1 + h·s2 (mod q).
425 hs2 := Mul(pk.H, s2)
426 computed := Add(s1, hs2)
427 return Equal(c, computed)
428 }
429