WIP Support for CRYSTALS-Dilithium (NIST PQC Round2)

https://pq-crystals.org/dilithium/data/dilithium-specification-round2.pdf

Notes

- At the moment, instead of using templates like other parts of circl,
  the dilithium code uses duplicated source files using symlinks.  This
  makes development a bit easier, but the layout a bit awkward.
- This implementation is based on the spec with the reference
  implementaiton used as, what is in a name, reference.  It turns out to
  be pretty close to the reference implementation.  I did not simply
  mimic their choices, but checked them (for instance, whether
  coefficients stay bounded.)  I did not find an issue except for some
  errors in comments.

To do

- [ ] Add public API
- [ ] Compare performance against C implementation
- [ ] Check some implementation details
- [ ] Document Make/UseHint
- [ ] Check whether the optimized checks on the signature (as also used
      in the reference implementation) correspond to those prescribed in
      the specification.

Nice to have

- [ ] Add AVX2 optimized NTT

Investigate

- [ ] Check how stupid the current Go compiler is --- it might be worthwhile
      to do one of the following.
    - [ ] Add explicit constants for literals like 1<<(D-1)
    - [ ] Manually unroll loops in Poly, VecK and VecL
    - [ ] Eliminate branches
- [ ] Check if it's an issue that a single public key has multiple valid
      packed representations, (see PackT1/UnpackT1.)

README.md

@@ -41,6 +41,7 @@ Version numbers are [Semvers](https://semver.org/). We release a minor version f
 | Key Exchange | FourQ | One of the fastest elliptic curves at 128-bit security level. | Experimental for key agreement and digital signatures. |
 | Key Exchange / Digital signatures | P-384 | Our optimizations reduce the burden when moving from P-256 to P-384. |  ECDSA and ECDH using Suite B at top secret level. |
 | Digital Signatures | Ed25519 | RFC-8032 provides new signature schemes based on Edwards curves. | Digital certificates and authentication. |
+| Digital Signatures | Dilithium | Lattice (Module LWE) based signature scheme | Post-Quantum PKI |
 
 ### Work in Progress
 

sign/dilithium/README.md

@@ -0,0 +1,6 @@
+dilithium
+=========
+
+Implementation of CRYSTALS-Dilithium
+[as proposed](https://pq-crystals.org/dilithium/data/dilithium-specification-round2.pdf)
+for the second round of the NIST PQC project.

sign/dilithium/internal/aes.go

@@ -0,0 +1,46 @@
+package internal
+
+import (
+	"crypto/aes"
+	"crypto/cipher"
+	"encoding/binary"
+)
+
+// AES CTR stream used as a replacement for SHAKE in Dilithium[1234]-AES.
+type AesStream struct {
+	c       cipher.Block
+	counter uint64
+	nonce   uint16
+}
+
+// Create a new AesStream as a replacement of SHAKE128.  (Well... it's not
+// replaced consistently.)
+func NewAesStream128(key *[32]byte, nonce uint16) AesStream {
+	c, _ := aes.NewCipher(key[:])
+	return AesStream{c: c, nonce: nonce}
+}
+
+// Create a new AesStream as a replacement of SHAKE256.  (Well... it's not
+// replaced consistently.)
+//
+// Yes, in an AES mode, Dilithium throws away the last 16 bytes of a seed ...
+// See the remark at the end of the caption of Figure 4 in the Round 2 spec.
+func NewAesStream256(key *[48]byte, nonce uint16) AesStream {
+	c, _ := aes.NewCipher(key[:32])
+	return AesStream{c: c, nonce: nonce}
+}
+
+// Squeeze some more blocks from the AES CTR stream into buf.
+//
+// Assumes length of buf is a multiple of 16.
+func (s *AesStream) SqueezeInto(buf []byte) {
+	var tmp [16]byte
+	binary.LittleEndian.PutUint16(tmp[:], s.nonce)
+
+	for len(buf) != 0 {
+		binary.BigEndian.PutUint64(tmp[8:], s.counter)
+		s.counter++
+		s.c.Encrypt(buf, tmp[:])
+		buf = buf[16:]
+	}
+}

sign/dilithium/internal/field.go

@@ -0,0 +1,105 @@
+package internal
+
+// Returns y with y < 2q and y = x mod q.
+func reduceLe2Q(x uint32) uint32 {
+	// Note 2^23 = 2^13 - 1 mod q. So, writing  x = x1 2^23 + x2 with x2 < 2^23
+	// and x1 < 2^9, we have x = y (mod q) where
+	// y = x2 + x1 2^13 - x1 ≤ 2^23 + 2^13 < 2q.
+	x1 := x >> 23
+	x2 := x & 0x7FFFFF // 2^23-1
+	return x2 + (x1 << 13) - x1
+}
+
+// Returns x mod q
+func modQ(x uint32) uint32 {
+	return le2qModQ(reduceLe2Q(x))
+}
+
+// For x R ≤ q 2^32, find y ≤ 2q with y = x mod q.
+func montReduceLe2Q(x uint64) uint32 {
+	// 4236238847 = -(q^-1) mod R = 2^32
+	m := (x * 4236238847) & 0xffffffff
+	return uint32((x + m*uint64(Q)) >> 32)
+}
+
+// Returns x mod q for 0 < x < 2q
+func le2qModQ(x uint32) uint32 {
+	x -= Q
+	mask := uint32(int32(x) >> 31) // mask is 2^32-1 if x was neg.; 0 otherwise
+	return x + (mask & Q)
+}
+
+// Splits 0 ≤ a < Q into a0 and a1 with a = a1*2^D + a0
+// and -2^{D-1} < a0 < 2^{D-1}.  Returns a0 + Q and a1.
+func power2round(a uint32) (a0plusQ, a1 uint32) {
+	// We effectively compute a0 = a mod± 2^d
+	//                    and a1 = (a - a0) / 2^d.
+	a0 := a & ((1 << D) - 1) // a mod 2^d
+
+	// a0 is one of  0, 1, ..., 2^{d-1}-1, 2^{d-1}, 2^{d-1}+1, ..., 2^d-1
+	a0 -= (1 << (D - 1)) + 1
+	// now a0 is     -2^{d-1}-1, -2^{d-1}, ..., -2, -1, 0, ..., 2^{d-1}-2
+	// Next, we add 2^D to those a0 that are negative (seen as int32).
+	a0 += uint32(int32(a0)>>31) & (1 << D)
+	// now a0 is     2^{d-1}-1, 2^{d-1}, ..., 2^d-2, 2^d-1, 0, ..., 2^{d-1}-2
+	a0 -= (1 << (D - 1)) - 1
+	// now a0 id     0, 1, 2, ..., 2^{d-1}-1, 2^{d-1}-1, -2^{d-1}-1, ...
+	// which is what we want.
+	a0plusQ = Q + a0
+	a1 = (a - a0) >> D
+	return
+}
+
+// Splits 0 ≤ a < Q into a0 and a1 with a = a1*α + a0 with -α/2 < a0 ≤ α/2,
+// except for when we would have a1 = (Q-1)/α in which case a1=0 is taken
+// and -α<2 ≤ a0 < 0.  Returns a0 + Q.  Note 0 ≤ a1 ≤ 15.
+func decompose(a uint32) (a0plusQ, a1 uint32) {
+	// Finds 0 ≤ t < 1.5α with t = a mod α.  (Recall α=2^19 - 2^9.)
+	t := int32(a & 0x7ffff)
+	t += (int32(a) >> 19) << 9
+
+	// If 0 ≤ t < α, then the following computes a mod± α with the same
+	// argument as in power2round().  If α ≤ t < 1.5α, then the following
+	// subtracts α, which thus also computes a mod± α.
+	t -= Alpha/2 + 1
+	t += (t >> 31) & Alpha
+	t -= Alpha/2 - 1
+
+	a1 = a - uint32(t)
+
+	// We want to divide α out of a1 (to get the proper value of a1).
+	// As our values are relatively small and α=2^19-2^9, we can simply
+	// divide by 2^19 and add one.  There is one corner case we have to deal
+	// with: if a1=0 then 0/α=0≠1=0/2^19+1, so we need to get rid of the +1.
+	u := ((a1 - 1) >> 31) & 1 // u=1 if a1=0
+	a1 = (a1 >> 19) + 1
+	a1 -= u // correct for the case a1=0
+
+	a0plusQ = Q + uint32(t)
+
+	// Now deal with the corner case of the definition, if a1=(Q-1)/α,
+	// then we use a1=0.  Note (Q-1)/α=2^4.
+	a0plusQ -= a1 >> 4 // to compensate, we only have to move the -1.
+	a1 &= 15           // set a0=0 if a1=16
+	return
+}
+
+// TODO document
+func makeHint(z0, r1 uint32) uint32 {
+	if z0 <= Gamma2 || z0 > Q-Gamma2 || (z0 == Q-Gamma2 && r1 == 0) {
+		return 0
+	}
+	return 1
+}
+
+// TODO document
+func useHint(r uint32, hint uint32) uint32 {
+	r0plusQ, r1 := decompose(r)
+	if hint == 0 {
+		return r1
+	}
+	if r0plusQ > Q {
+		return (r1 + 1) & 15
+	}
+	return (r1 - 1) & 15
+}

sign/dilithium/internal/field_test.go

@@ -0,0 +1,69 @@
+package internal
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func TestModQ(t *testing.T) {
+	for i := 0; i < 1000; i++ {
+		x := uint32(rand.Int31())
+		y := modQ(x)
+		if y > Q {
+			t.Fatalf("modQ(%d) > Q", x)
+		}
+		if y != x%Q {
+			t.Fatalf("modQ(%d) != %d (mod Q)", x, x)
+		}
+	}
+}
+
+func TestReduceLe2Q(t *testing.T) {
+	for i := 0; i < 1000; i++ {
+		x := uint32(rand.Int31())
+		y := reduceLe2Q(x)
+		if y > 2*Q {
+			t.Fatalf("recuce_le2q(%d) > 2Q", x)
+		}
+		if y%Q != x%Q {
+			t.Fatalf("recuce_le2q(%d) != %d (mod Q)", x, x)
+		}
+	}
+}
+
+func TestPower2Round(t *testing.T) {
+	for a := uint32(0); a < Q; a++ {
+		a0PlusQ, a1 := power2round(a)
+		a0 := int32(a0PlusQ) - int32(Q)
+		if int32(a) != a0+int32((1<<D)*a1) {
+			t.Fatalf("power2round(%v) doesn't recombine", a)
+		}
+		if (-(1 << (D - 1)) >= a0) || (a0 > 1<<(D-1)) {
+			t.Fatalf("power2round(%v): a0 out of bounds", a)
+		}
+		if a1 > (1 << (23 - 14)) { // 23 ~ 2log Q. XXX
+			t.Fatalf("power2round(%v): a1 out of bounds", a)
+		}
+	}
+}
+
+func TestDecompose(t *testing.T) {
+	for a := uint32(0); a < Q; a++ {
+		a0PlusQ, a1 := decompose(a)
+		a0 := int32(a0PlusQ) - int32(Q)
+		recombined := a0 + int32(Alpha*a1)
+		if a1 == 0 && recombined < 0 {
+			recombined += Q
+			if -(Alpha/2) > a0 || a0 >= 0 {
+				t.Fatalf("decompose(%v): a0 out of bounds", a)
+			}
+		} else {
+			if (-(Alpha / 2) >= a0) || (a0 > Alpha/2) {
+				t.Fatalf("decompose(%v): a0 out of bounds", a)
+			}
+		}
+		if int32(a) != recombined {
+			t.Fatalf("decompose(%v) doesn't recombine %v %v", a, a0, a1)
+		}
+	}
+}