diff --git a/README.md b/README.md
index ab9fd5d65..8851227d2 100644
--- a/README.md
+++ b/README.md
@@ -40,6 +40,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
 
diff --git a/sign/dilithium/README.md b/sign/dilithium/README.md
new file mode 100644
index 000000000..91ad4ae8c
--- /dev/null
+++ b/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.
diff --git a/sign/dilithium/mode1/internal/dilithium.go b/sign/dilithium/mode1/internal/dilithium.go
new file mode 120000
index 000000000..f4eea29aa
--- /dev/null
+++ b/sign/dilithium/mode1/internal/dilithium.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/dilithium_test.go b/sign/dilithium/mode1/internal/dilithium_test.go
new file mode 120000
index 000000000..baa8508b5
--- /dev/null
+++ b/sign/dilithium/mode1/internal/dilithium_test.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/field.go b/sign/dilithium/mode1/internal/field.go
new file mode 120000
index 000000000..2a2a32e6d
--- /dev/null
+++ b/sign/dilithium/mode1/internal/field.go
@@ -0,0 +1 @@
+../../mode3/internal/field.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/field_test.go b/sign/dilithium/mode1/internal/field_test.go
new file mode 120000
index 000000000..3e9043a03
--- /dev/null
+++ b/sign/dilithium/mode1/internal/field_test.go
@@ -0,0 +1 @@
+../../mode3/internal/field_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/kat_test.go b/sign/dilithium/mode1/internal/kat_test.go
new file mode 120000
index 000000000..54cd221f1
--- /dev/null
+++ b/sign/dilithium/mode1/internal/kat_test.go
@@ -0,0 +1 @@
+../../mode3/internal/kat_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/mat.go b/sign/dilithium/mode1/internal/mat.go
new file mode 120000
index 000000000..0601aa926
--- /dev/null
+++ b/sign/dilithium/mode1/internal/mat.go
@@ -0,0 +1 @@
+../../mode3/internal/mat.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/ntt.go b/sign/dilithium/mode1/internal/ntt.go
new file mode 120000
index 000000000..9b9aa00b8
--- /dev/null
+++ b/sign/dilithium/mode1/internal/ntt.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/ntt_test.go b/sign/dilithium/mode1/internal/ntt_test.go
new file mode 120000
index 000000000..02efd77bc
--- /dev/null
+++ b/sign/dilithium/mode1/internal/ntt_test.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/pack.go b/sign/dilithium/mode1/internal/pack.go
new file mode 120000
index 000000000..578d3dc68
--- /dev/null
+++ b/sign/dilithium/mode1/internal/pack.go
@@ -0,0 +1 @@
+../../mode3/internal/pack.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/pack_test.go b/sign/dilithium/mode1/internal/pack_test.go
new file mode 120000
index 000000000..5571475bb
--- /dev/null
+++ b/sign/dilithium/mode1/internal/pack_test.go
@@ -0,0 +1 @@
+../../mode3/internal/pack_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/params.go b/sign/dilithium/mode1/internal/params.go
new file mode 100644
index 000000000..e16606bc3
--- /dev/null
+++ b/sign/dilithium/mode1/internal/params.go
@@ -0,0 +1,14 @@
+package internal
+
+const (
+	Name           = "Dilithium1"
+	PublicKeySize  = 896
+	PrivateKeySize = 2096
+	SignatureSize  = 1387
+	K              = 3
+	L              = 2
+	Eta            = 7
+	DoubleEtaBits  = 4
+	Beta           = 375
+	Omega          = 64
+)
diff --git a/sign/dilithium/mode1/internal/params_test.go b/sign/dilithium/mode1/internal/params_test.go
new file mode 100644
index 000000000..9f4688516
--- /dev/null
+++ b/sign/dilithium/mode1/internal/params_test.go
@@ -0,0 +1,85 @@
+package internal
+
+import (
+	"testing"
+
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+// Tests specific to the current mode
+
+func hash(in []byte) [16]byte {
+	var ret [16]byte
+	h := shake.NewShake256()
+	h.Write(in)
+	h.Read(ret[:])
+	return ret
+}
+
+func TestNewKeyFromSeed(t *testing.T) {
+	var seed [SeedSize]byte
+	var pkp [PublicKeySize]byte
+	var skp [PrivateKeySize]byte
+	pk, sk := NewKeyFromSeed(&seed)
+	pk.Pack(&pkp)
+	sk.Pack(&skp)
+
+	// Generated with reference implementation.
+	ehpk := [16]byte{
+		131, 97, 105, 81, 185, 131, 18, 169,
+		126, 161, 14, 18, 183, 182, 150, 117,
+	}
+	ehsk := [16]byte{
+		175, 71, 14, 18, 165, 125, 0, 192, 76,
+		74, 43, 89, 152, 244, 28, 113,
+	}
+
+	if hash(pkp[:]) != ehpk {
+		t.Fatalf("pk not ok")
+	}
+	if hash(skp[:]) != ehsk {
+		t.Fatalf("sk not ok")
+	}
+}
+
+func TestVectorDeriveUniformLeqEta(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	p2 := Poly{
+		3, 0, 2, 3, 8380412, 8380414, 2, 2, 8380415, 8380416,
+		8380416, 1, 7, 8380411, 4, 7, 2, 8380414, 8380412, 0, 3,
+		8380413, 4, 8380410, 8380414, 8380410, 8380413, 8380412,
+		8380410, 8380410, 7, 3, 8380412, 3, 8380411, 8380414, 7,
+		8380415, 1, 8380416, 7, 8380414, 2, 4, 7, 8380416, 4, 5,
+		5, 2, 1, 8380411, 7, 3, 4, 8380416, 8380415, 2, 2, 8380410,
+		7, 5, 0, 8380411, 8380416, 7, 8380412, 5, 8380411, 5, 4,
+		1, 4, 2, 8380414, 8380416, 3, 4, 8380414, 8380410, 8380414,
+		7, 8380415, 8380412, 8380411, 6, 8380415, 8380410, 8380413,
+		0, 8380416, 8380412, 8380411, 6, 8380415, 8380410, 8380410,
+		8380413, 7, 8380415, 3, 1, 0, 8380414, 4, 8380410, 8380412,
+		5, 6, 0, 1, 8380410, 8380413, 8380410, 6, 8380411, 0,
+		8380415, 3, 8380411, 2, 3, 6, 8380410, 8380410, 0, 0, 6,
+		8380414, 0, 6, 8380411, 8380412, 8380412, 5, 7, 2, 8380413,
+		8380410, 0, 0, 4, 0, 5, 8380416, 1, 5, 7, 7, 8380410, 2,
+		0, 8380412, 6, 5, 4, 8380411, 0, 0, 7, 8380414, 8380416,
+		8380410, 4, 0, 8380411, 8380410, 8380413, 8380414, 3,
+		8380413, 8380414, 8380415, 7, 8380411, 8380414, 8380410,
+		5, 6, 8380410, 6, 8380415, 4, 0, 0, 8380416, 8380415,
+		8380415, 8380415, 4, 8380410, 8380411, 8380415, 8380412,
+		6, 1, 4, 1, 2, 8380411, 0, 8380416, 4, 8380414, 8380416,
+		6, 5, 1, 3, 8380414, 4, 8380414, 8380414, 3, 2, 2, 8380413,
+		2, 8380412, 0, 8380416, 8380415, 8380414, 0, 3, 8380413,
+		3, 8380410, 7, 8380410, 8380414, 8380412, 1, 4, 2, 8380415,
+		2, 8380410, 8380410, 8380410, 4, 8380413, 4, 4, 0, 3,
+		8380410, 5, 8380413, 8380414, 8380415, 4, 1, 8380410, 7,
+		8380413,
+	}
+	for i := 0; i < 32; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniformLeqEta(&seed, 30000)
+	p.Normalize()
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
diff --git a/sign/dilithium/mode1/internal/poly.go b/sign/dilithium/mode1/internal/poly.go
new file mode 120000
index 000000000..67843bdd7
--- /dev/null
+++ b/sign/dilithium/mode1/internal/poly.go
@@ -0,0 +1 @@
+../../mode3/internal/poly.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/sample.go b/sign/dilithium/mode1/internal/sample.go
new file mode 120000
index 000000000..7db77123c
--- /dev/null
+++ b/sign/dilithium/mode1/internal/sample.go
@@ -0,0 +1 @@
+../../mode3/internal/sample.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/sample_test.go b/sign/dilithium/mode1/internal/sample_test.go
new file mode 120000
index 000000000..a7d78fe33
--- /dev/null
+++ b/sign/dilithium/mode1/internal/sample_test.go
@@ -0,0 +1 @@
+../../mode3/internal/sample_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode1/internal/vec.go b/sign/dilithium/mode1/internal/vec.go
new file mode 120000
index 000000000..82e7f9b35
--- /dev/null
+++ b/sign/dilithium/mode1/internal/vec.go
@@ -0,0 +1 @@
+../../mode3/internal/vec.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/dilithium.go b/sign/dilithium/mode2/internal/dilithium.go
new file mode 120000
index 000000000..f4eea29aa
--- /dev/null
+++ b/sign/dilithium/mode2/internal/dilithium.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/dilithium_test.go b/sign/dilithium/mode2/internal/dilithium_test.go
new file mode 120000
index 000000000..baa8508b5
--- /dev/null
+++ b/sign/dilithium/mode2/internal/dilithium_test.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/field.go b/sign/dilithium/mode2/internal/field.go
new file mode 120000
index 000000000..2a2a32e6d
--- /dev/null
+++ b/sign/dilithium/mode2/internal/field.go
@@ -0,0 +1 @@
+../../mode3/internal/field.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/field_test.go b/sign/dilithium/mode2/internal/field_test.go
new file mode 120000
index 000000000..3e9043a03
--- /dev/null
+++ b/sign/dilithium/mode2/internal/field_test.go
@@ -0,0 +1 @@
+../../mode3/internal/field_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/kat_test.go b/sign/dilithium/mode2/internal/kat_test.go
new file mode 120000
index 000000000..54cd221f1
--- /dev/null
+++ b/sign/dilithium/mode2/internal/kat_test.go
@@ -0,0 +1 @@
+../../mode3/internal/kat_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/mat.go b/sign/dilithium/mode2/internal/mat.go
new file mode 120000
index 000000000..0601aa926
--- /dev/null
+++ b/sign/dilithium/mode2/internal/mat.go
@@ -0,0 +1 @@
+../../mode3/internal/mat.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/ntt.go b/sign/dilithium/mode2/internal/ntt.go
new file mode 120000
index 000000000..9b9aa00b8
--- /dev/null
+++ b/sign/dilithium/mode2/internal/ntt.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/ntt_test.go b/sign/dilithium/mode2/internal/ntt_test.go
new file mode 120000
index 000000000..02efd77bc
--- /dev/null
+++ b/sign/dilithium/mode2/internal/ntt_test.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/pack.go b/sign/dilithium/mode2/internal/pack.go
new file mode 120000
index 000000000..578d3dc68
--- /dev/null
+++ b/sign/dilithium/mode2/internal/pack.go
@@ -0,0 +1 @@
+../../mode3/internal/pack.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/pack_test.go b/sign/dilithium/mode2/internal/pack_test.go
new file mode 120000
index 000000000..5571475bb
--- /dev/null
+++ b/sign/dilithium/mode2/internal/pack_test.go
@@ -0,0 +1 @@
+../../mode3/internal/pack_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/params.go b/sign/dilithium/mode2/internal/params.go
new file mode 100644
index 000000000..bbf530159
--- /dev/null
+++ b/sign/dilithium/mode2/internal/params.go
@@ -0,0 +1,14 @@
+package internal
+
+const (
+	Name           = "Dilithium2"
+	PublicKeySize  = 1184
+	PrivateKeySize = 2800
+	SignatureSize  = 2044
+	K              = 4
+	L              = 3
+	Eta            = 6
+	DoubleEtaBits  = 4
+	Beta           = 325
+	Omega          = 80
+)
diff --git a/sign/dilithium/mode2/internal/params_test.go b/sign/dilithium/mode2/internal/params_test.go
new file mode 100644
index 000000000..64a32e113
--- /dev/null
+++ b/sign/dilithium/mode2/internal/params_test.go
@@ -0,0 +1,85 @@
+package internal
+
+import (
+	"testing"
+
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+// Tests specific to the current mode
+
+func hash(in []byte) [16]byte {
+	var ret [16]byte
+	h := shake.NewShake256()
+	h.Write(in)
+	h.Read(ret[:])
+	return ret
+}
+
+func TestNewKeyFromSeed(t *testing.T) {
+	var seed [SeedSize]byte
+	var pkp [PublicKeySize]byte
+	var skp [PrivateKeySize]byte
+	pk, sk := NewKeyFromSeed(&seed)
+	pk.Pack(&pkp)
+	sk.Pack(&skp)
+
+	// Generated with reference implementation.
+	ehpk := [16]byte{
+		56, 231, 51, 157, 0, 230, 67, 72, 203,
+		47, 150, 94, 207, 158, 227, 139,
+	}
+	ehsk := [16]byte{
+		72, 222, 195, 214, 136, 51, 13, 252,
+		104, 249, 191, 66, 119, 251, 146, 225,
+	}
+
+	if hash(pkp[:]) != ehpk {
+		t.Fatalf("pk not ok")
+	}
+	if hash(skp[:]) != ehsk {
+		t.Fatalf("sk not ok")
+	}
+}
+
+func TestVectorDeriveUniformLeqEta(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	p2 := Poly{
+		2, 8380416, 1, 2, 8380411, 8380413, 1, 1, 8380414, 8380415,
+		8380415, 0, 6, 3, 6, 1, 8380413, 8380411, 8380416, 2,
+		8380412, 3, 8380413, 8380412, 8380411, 6, 2, 8380411, 2,
+		8380413, 6, 8380414, 0, 8380415, 6, 8380413, 1, 3, 6,
+		8380415, 3, 4, 4, 1, 0, 6, 2, 3, 8380415, 8380414, 1, 1,
+		6, 4, 8380416, 8380415, 6, 8380411, 4, 4, 3, 0, 3, 1,
+		8380413, 8380415, 2, 3, 8380413, 8380413, 6, 8380414,
+		8380411, 5, 8380414, 8380412, 8380416, 8380415, 8380411,
+		5, 8380414, 8380412, 6, 8380414, 2, 0, 8380416, 8380413,
+		3, 8380411, 4, 5, 8380416, 0, 8380412, 5, 8380416, 8380414,
+		2, 1, 2, 5, 8380416, 8380416, 5, 8380413, 8380416, 5,
+		8380411, 8380411, 4, 6, 1, 8380412, 8380416, 8380416, 3,
+		8380416, 4, 8380415, 0, 4, 6, 6, 1, 8380416, 8380411, 5,
+		4, 3, 8380416, 8380416, 6, 8380413, 8380415, 3, 8380416,
+		8380412, 8380413, 2, 8380412, 8380413, 8380414, 6, 8380413,
+		4, 5, 5, 8380414, 3, 8380416, 8380416, 8380415, 8380414,
+		8380414, 8380414, 3, 8380414, 8380411, 5, 0, 3, 0, 1,
+		8380416, 8380415, 3, 8380413, 8380415, 5, 4, 0, 2, 8380413,
+		3, 8380413, 8380413, 2, 1, 1, 8380412, 1, 8380411, 8380416,
+		8380415, 8380414, 8380413, 8380416, 2, 8380412, 2, 6,
+		8380413, 8380411, 0, 3, 1, 8380414, 1, 3, 8380412, 3, 3,
+		8380416, 2, 4, 8380412, 8380413, 8380414, 3, 0, 6, 8380412,
+		8380415, 8380416, 5, 5, 0, 6, 8380415, 8380413, 8380414,
+		5, 6, 3, 6, 3, 5, 5, 8380415, 6, 4, 8380416, 3, 2, 6,
+		8380414, 2, 8380415, 6, 8380412, 8380413, 8380415, 1, 2,
+		0, 8380411, 6, 8380413, 6, 8380414, 8380416, 8380415,
+		8380414, 1, 8380412,
+	}
+	for i := 0; i < 32; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniformLeqEta(&seed, 30000)
+	p.Normalize()
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
diff --git a/sign/dilithium/mode2/internal/poly.go b/sign/dilithium/mode2/internal/poly.go
new file mode 120000
index 000000000..67843bdd7
--- /dev/null
+++ b/sign/dilithium/mode2/internal/poly.go
@@ -0,0 +1 @@
+../../mode3/internal/poly.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/sample.go b/sign/dilithium/mode2/internal/sample.go
new file mode 120000
index 000000000..7db77123c
--- /dev/null
+++ b/sign/dilithium/mode2/internal/sample.go
@@ -0,0 +1 @@
+../../mode3/internal/sample.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/sample_test.go b/sign/dilithium/mode2/internal/sample_test.go
new file mode 120000
index 000000000..a7d78fe33
--- /dev/null
+++ b/sign/dilithium/mode2/internal/sample_test.go
@@ -0,0 +1 @@
+../../mode3/internal/sample_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode2/internal/vec.go b/sign/dilithium/mode2/internal/vec.go
new file mode 120000
index 000000000..82e7f9b35
--- /dev/null
+++ b/sign/dilithium/mode2/internal/vec.go
@@ -0,0 +1 @@
+../../mode3/internal/vec.go
\ No newline at end of file
diff --git a/sign/dilithium/mode3/internal/dilithium.go b/sign/dilithium/mode3/internal/dilithium.go
new file mode 100644
index 000000000..0fb4ad243
--- /dev/null
+++ b/sign/dilithium/mode3/internal/dilithium.go
@@ -0,0 +1,382 @@
+package internal
+
+import (
+	cryptoRand "crypto/rand"
+	"io"
+
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+const (
+	SeedSize = 32
+	N        = 256
+	Q        = 8380417
+	QBits    = 23
+	D        = 14
+	Gamma1   = ((Q - 1) / 16)
+	Gamma2   = Gamma1 / 2
+	Alpha    = 2 * Gamma2
+
+	// Size of a packed polynomial of norm ≤η.
+	// (Note that the  formula is not valid in general.)
+	PolyLeqEtaSize = (N * DoubleEtaBits) / 8
+
+	// Size of T1 packed.  (Note that the formula is not valid in general.)
+	PolyT1Size = (N * (QBits - D)) / 8
+
+	// Size of T0 packed.  (Note that the formula is not valid in general.)
+	PolyT0Size = (N * D) / 8
+
+	// Size of a packed polynomial of norm <γ1.
+	PolyLeGamma1Size = (N * (QBits - 3)) / 8
+
+	// Size of a packed polynomial whose coeffients are in [0,16).
+	PolyLe16Size = N / 2
+)
+
+// PublicKey is the type of Dilithium public keys.
+type PublicKey struct {
+	rho [32]byte
+	t1  VecK
+
+	// Cached values
+	t1p [PolyT1Size * K]byte
+}
+
+// PublicKey is the type of Dilithium private keys.
+type PrivateKey struct {
+	rho [32]byte
+	key [32]byte
+	s1  VecL
+	s2  VecK
+	t0  VecK
+	tr  [48]byte
+}
+
+type unpackedSignature struct {
+	z    VecL
+	hint VecK
+	c    Poly
+}
+
+// Packs the signature into buf.
+func (sig *unpackedSignature) Pack(buf []byte) {
+	sig.z.PackLeGamma1(buf[:])
+	sig.hint.PackHint(buf[L*PolyLeGamma1Size:])
+	sig.c.PackB60(buf[L*PolyLeGamma1Size+Omega+K:])
+}
+
+// Sets sig to the signature encoded in the buffer.
+//
+// Returns whether buf contains a properly packed signature.
+func (sig *unpackedSignature) Unpack(buf []byte) bool {
+	if len(buf) < SignatureSize {
+		return false
+	}
+	sig.z.UnpackLeGamma1(buf[:])
+	if sig.z.Exceeds(Gamma1 - Beta) {
+		return false
+	}
+	if !sig.hint.UnpackHint(buf[L*PolyLeGamma1Size:]) {
+		return false
+	}
+	sig.c.UnpackB60(buf[L*PolyLeGamma1Size+Omega+K:])
+	return true
+}
+
+// Packs the public key into buf.
+func (pk *PublicKey) Pack(buf *[PublicKeySize]byte) {
+	copy(buf[:32], pk.rho[:])
+	copy(buf[32:], pk.t1p[:])
+}
+
+// Sets pk to the public key encoded in buf.
+func (pk *PublicKey) Unpack(buf *[PublicKeySize]byte) {
+	copy(pk.rho[:], buf[:32])
+	copy(pk.t1p[:], buf[32:])
+	pk.t1.UnpackT1(pk.t1p[:])
+}
+
+// Packs the private key into buf.
+func (sk *PrivateKey) Pack(buf *[PrivateKeySize]byte) {
+	copy(buf[:32], sk.rho[:])
+	copy(buf[32:64], sk.key[:])
+	copy(buf[64:112], sk.tr[:])
+	offset := 112
+	sk.s1.PackLeqEta(buf[offset:])
+	offset += PolyLeqEtaSize * L
+	sk.s2.PackLeqEta(buf[offset:])
+	offset += PolyLeqEtaSize * K
+	sk.t0.PackT0(buf[offset:])
+}
+
+// Sets sk to the private key encoded in buf.
+func (sk *PrivateKey) Unpack(buf *[PrivateKeySize]byte) {
+	copy(sk.rho[:], buf[:32])
+	copy(sk.key[:], buf[32:64])
+	copy(sk.tr[:], buf[64:112])
+	offset := 112
+	sk.s1.UnpackLeqEta(buf[offset:])
+	offset += PolyLeqEtaSize * L
+	sk.s2.UnpackLeqEta(buf[offset:])
+	offset += PolyLeqEtaSize * K
+	sk.t0.UnpackT0(buf[offset:])
+}
+
+// TODO implement Signer
+// TODO cache A
+// TODO cache tr
+
+// GenerateKey generates a public/private key pair using entropy from rand.
+// If rand is nil, crypto/rand.Reader will be used.
+func GenerateKey(rand io.Reader) (*PublicKey, *PrivateKey, error) {
+	var seed [32]byte
+	if rand == nil {
+		rand = cryptoRand.Reader
+	}
+	_, err := io.ReadFull(rand, seed[:])
+	if err != nil {
+		return nil, nil, err
+	}
+	pk, sk := NewKeyFromSeed(&seed)
+	return pk, sk, nil
+}
+
+// NewKeyFromExpandedSeed derives a public/private key pair using the
+// given expanded seed.
+//
+// Use NewKeyFromSeed instead of this function.  This function is only exposed
+// to generate the NIST KAT test vectors.
+func NewKeyFromExpandedSeed(seed *[96]byte) (*PublicKey, *PrivateKey) {
+	var pk PublicKey
+	var sk PrivateKey
+	var sSeed [32]byte
+	var t VecK
+
+	var A Mat
+
+	copy(pk.rho[:], seed[:32])
+	copy(sSeed[:], seed[32:64])
+	copy(sk.key[:], seed[64:])
+	copy(sk.rho[:], pk.rho[:])
+
+	A.Derive(&pk.rho)
+
+	for i := uint16(0); i < L; i++ {
+		sk.s1[i].DeriveUniformLeqEta(&sSeed, i)
+	}
+
+	for i := uint16(0); i < K; i++ {
+		sk.s2[i].DeriveUniformLeqEta(&sSeed, i+L)
+	}
+
+	// Set t to A s1 + s2
+	s1h := sk.s1
+	s1h.NTT()
+	for i := 0; i < K; i++ {
+		t[i].DotHat(&A[i], &s1h)
+		t[i].ReduceLe2Q()
+		t[i].InvNTT()
+	}
+	t.Add(&t, &sk.s2)
+	t.Normalize()
+
+	// Compute t0, t1 = Power2Round(t)
+	t.Power2Round(&sk.t0, &pk.t1)
+
+	// Finish public key
+	pk.t1.PackT1(pk.t1p[:])
+
+	var packedPk [PublicKeySize]byte
+	pk.Pack(&packedPk)
+	h := shake.NewShake256()
+	h.Write(packedPk[:])
+	h.Read(sk.tr[:])
+
+	return &pk, &sk
+}
+
+// NewKeyFromSeed derives a public/private key pair using the given seed.
+func NewKeyFromSeed(seed *[SeedSize]byte) (*PublicKey, *PrivateKey) {
+	var buf [96]byte
+	h := shake.NewShake128()
+	h.Write(seed[:])
+	h.Read(buf[:])
+	return NewKeyFromExpandedSeed(&buf)
+}
+
+// Verify checks whether sig is a signature on msg by pk.
+func Verify(pk *PublicKey, msg []byte, signature []byte) bool {
+	var sig unpackedSignature
+	var A Mat
+	var tr, mu [48]byte
+	var zh VecL
+	var Az, Az2dct1, w1 VecK
+	var ch, cp Poly
+
+	// Note that Unpack() checked whether ‖z‖_∞ < γ_1 - β
+	// and ensured that there at most ω ones in pk.hint.
+	if !sig.Unpack(signature) {
+		return false
+	}
+
+	A.Derive(&pk.rho)
+
+	// tr = CRH(ρ ‖ t1)
+	h := shake.NewShake256()
+	h.Write(pk.rho[:])
+	h.Write(pk.t1p[:])
+	h.Read(tr[:])
+
+	// μ = CRH(tr ‖ msg)
+	h.Reset()
+	h.Write(tr[:])
+	h.Write(msg)
+	h.Read(mu[:])
+
+	// Compute Az
+	zh = sig.z
+	zh.NTT()
+
+	for i := 0; i < K; i++ {
+		Az[i].DotHat(&A[i], &zh)
+	}
+
+	// Next, we compute Az - 2^d·c·t1.
+	// Note that the coefficients of t1 are bounded by 256 = 2^9,
+	// so the coefficients of Az2dct1 will bounded by 2^{9+d} = 2^23 < 2q,
+	// which is small enough for NTT().
+	Az2dct1.MulBy2toD(&pk.t1)
+	Az2dct1.NTT()
+	ch = sig.c
+	ch.NTT()
+	for i := 0; i < K; i++ {
+		Az2dct1[i].MulHat(&Az2dct1[i], &ch)
+	}
+	Az2dct1.Sub(&Az, &Az2dct1)
+	Az2dct1.ReduceLe2Q()
+	Az2dct1.InvNTT()
+	Az2dct1.NormalizeAssumingLe2Q()
+
+	// w1 = UseHint(pk.hint, Az - 2^d·c·t1).
+	w1.UseHint(&Az2dct1, &sig.hint)
+
+	// c' = H(μ, w1)
+	cp.DeriveUniformB60(&mu, &w1)
+	if sig.c != cp {
+		return false
+	}
+
+	return true
+}
+
+// Sign signs the given message and writes the signature into signature
+func SignTo(sk *PrivateKey, msg []byte, signature []byte) {
+	var A Mat
+	var mu, rhop [48]byte
+	var s1h, y, yh VecL
+	var s2h, t0h, w, w0, w1, w0mcs2, ct0, w0mcs2pct0 VecK
+	var ch Poly
+	var yNonce uint16
+	var sig unpackedSignature
+
+	A.Derive(&sk.rho)
+
+	//  μ = CRH(tr ‖ msg)
+	h := shake.NewShake256()
+	h.Write(sk.tr[:])
+	h.Write(msg)
+	h.Read(mu[:])
+
+	// ρ' = CRH(μ ‖ key)
+	h.Reset()
+	h.Write(sk.key[:])
+	h.Write(mu[:])
+	h.Read(rhop[:])
+
+	// Precompute NTT(s1), NTT(s2) ad NTT(t0)
+	s1h = sk.s1
+	s1h.NTT()
+	s2h = sk.s2
+	s2h.NTT()
+	t0h = sk.t0
+	t0h.NTT()
+
+	// Main rejection loop
+	for {
+		// y = ExpandMask(ρ', key)
+		for i := 0; i < L; i++ {
+			y[i].DeriveUniformLeGamma1(&rhop, yNonce)
+			yNonce++
+		}
+
+		// Set w to A y
+		yh = y
+		yh.NTT()
+		for i := 0; i < K; i++ {
+			w[i].DotHat(&A[i], &yh)
+			w[i].ReduceLe2Q()
+			w[i].InvNTT()
+		}
+
+		// Decompose w into w0 and w1
+		w.NormalizeAssumingLe2Q()
+		w.Decompose(&w0, &w1)
+
+		// c = H(μ, w1)
+		sig.c.DeriveUniformB60(&mu, &w1)
+		ch = sig.c
+		ch.NTT()
+
+		// TODO check reference implementation reduction of checks to these
+		//      three cases.
+		// Ensure ‖ w0 - c·s2 ‖_∞ < γ2 - β.  See Lemma 2 in the spec.
+		for i := 0; i < K; i++ {
+			w0mcs2[i].MulHat(&ch, &s2h[i])
+			w0mcs2[i].InvNTT()
+		}
+		w0mcs2.Sub(&w0, &w0mcs2)
+		w0mcs2.Normalize()
+
+		if w0mcs2.Exceeds(Gamma2 - Beta) {
+			continue
+		}
+
+		// z = y + c·s1
+		for i := 0; i < L; i++ {
+			sig.z[i].MulHat(&ch, &s1h[i])
+			sig.z[i].InvNTT()
+		}
+		sig.z.Add(&sig.z, &y)
+		sig.z.Normalize()
+
+		// Ensure  ‖z‖_∞ < γ1 - β
+		if sig.z.Exceeds(Gamma1 - Beta) {
+			continue
+		}
+
+		// Compute c·t0
+		for i := 0; i < K; i++ {
+			ct0[i].MulHat(&ch, &t0h[i])
+			ct0[i].InvNTT()
+		}
+		ct0.NormalizeAssumingLe2Q()
+
+		// Ensure ‖c·t0‖_∞ < γ2.
+		if ct0.Exceeds(Gamma2) {
+			continue
+		}
+
+		// Create the hint
+		w0mcs2pct0.Add(&w0mcs2, &ct0)
+		w0mcs2pct0.NormalizeAssumingLe2Q()
+		hintPop := sig.hint.MakeHint(&w0mcs2pct0, &w1)
+		if hintPop > Omega {
+			continue
+		}
+
+		break
+	}
+
+	sig.Pack(signature[:])
+}
diff --git a/sign/dilithium/mode3/internal/dilithium_test.go b/sign/dilithium/mode3/internal/dilithium_test.go
new file mode 100644
index 000000000..fe3a4500d
--- /dev/null
+++ b/sign/dilithium/mode3/internal/dilithium_test.go
@@ -0,0 +1,89 @@
+package internal
+
+import (
+	"encoding/binary"
+	"testing"
+)
+
+// Checks whether p is normalized.  Only used in tests.
+func (p *Poly) Normalized() bool {
+	p2 := *p
+	p2.Normalize()
+	return p2 == *p
+}
+
+func BenchmarkSign(b *testing.B) {
+	var seed [32]byte
+	var msg [8]byte
+	var sig [SignatureSize]byte
+	_, sk := NewKeyFromSeed(&seed)
+	for i := 0; i < b.N; i++ {
+		binary.LittleEndian.PutUint64(msg[:], uint64(i))
+		SignTo(sk, msg[:], sig[:])
+	}
+}
+
+func BenchmarkGenerateKey(b *testing.B) {
+	var seed [32]byte
+	for i := 0; i < b.N; i++ {
+		NewKeyFromSeed(&seed)
+	}
+}
+
+func TestSign(t *testing.T) {
+	var seed [SeedSize]byte
+	var sig [SignatureSize]byte
+	msg := []byte{
+		86, 101, 114, 105, 108, 121, 44, 32, 118, 101, 114,
+		105, 108, 121, 44, 32, 73, 32, 115, 97, 121, 32,
+		117, 110, 116, 111, 32, 121, 111, 117, 44, 32, 208,
+		181, 120, 99, 101, 112, 116, 32, 97, 32, 99, 111,
+		114, 110, 32, 111, 102, 32, 119, 104, 101, 97, 116,
+		32, 102, 97, 108, 108, 32, 105, 110, 116, 111, 32,
+		116, 104, 101, 32, 103, 114, 111, 117, 110, 100,
+		32, 97, 110, 100, 32, 100, 105, 101, 44, 32, 105,
+		116, 32, 97, 98, 105, 100, 101, 116, 104, 32, 97,
+		108, 111, 110, 101, 58, 32, 98, 117, 116, 32, 105,
+		102, 32, 105, 116, 32, 100, 105, 101, 44, 32, 105,
+		116, 32, 98, 114, 105, 110, 103, 101, 116, 104, 32,
+		102, 111, 114, 116, 104, 32, 109, 117, 99, 104, 32,
+		102, 114, 117, 105, 116, 46,
+	}
+
+	pk, sk := NewKeyFromSeed(&seed)
+	SignTo(sk, msg, sig[:])
+	if !Verify(pk, msg, sig[:]) {
+		t.Fatal()
+	}
+}
+
+func TestSignThenVerifyAndPkSkPacking(t *testing.T) {
+	var seed [SeedSize]byte
+	var sig [SignatureSize]byte
+	var msg [8]byte
+	var pkb [PublicKeySize]byte
+	var skb [PrivateKeySize]byte
+	var pk2 PublicKey
+	var sk2 PrivateKey
+	for i := uint64(0); i < 100; i++ {
+		binary.LittleEndian.PutUint64(seed[:], i)
+		pk, sk := NewKeyFromSeed(&seed)
+		for j := uint64(0); j < 10; j++ {
+			binary.LittleEndian.PutUint64(msg[:], j)
+			SignTo(sk, msg[:], sig[:])
+			if !Verify(pk, msg[:], sig[:]) {
+				t.Fatal()
+			}
+		}
+		pk.Pack(&pkb)
+		pk2.Unpack(&pkb)
+		if *pk != pk2 {
+			t.Fatal()
+		}
+		sk.Pack(&skb)
+		sk2.Unpack(&skb)
+		if *sk != sk2 {
+			t.Fatal()
+		}
+	}
+}
diff --git a/sign/dilithium/mode3/internal/field.go b/sign/dilithium/mode3/internal/field.go
new file mode 100644
index 000000000..6d2c42c79
--- /dev/null
+++ b/sign/dilithium/mode3/internal/field.go
@@ -0,0 +1,105 @@
+package internal
+
+// Returns y with y < 2q and y = x mod q.
+func reduce_le2q(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 mod_q(x uint32) uint32 {
+	return le2q_mod_q(reduce_le2q(x))
+}
+
+// For x R ≤ q 2^32, find y ≤ 2q with y = x mod q.
+func mont_reduce_le2q(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 le2q_mod_q(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
+	// substracts α, 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 make_hint(z0, r1 uint32) uint32 {
+	if z0 <= Gamma2 || z0 > Q-Gamma2 || (z0 == Q-Gamma2 && r1 == 0) {
+		return 0
+	}
+	return 1
+}
+
+// TODO document
+func use_hint(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
+}
diff --git a/sign/dilithium/mode3/internal/field_test.go b/sign/dilithium/mode3/internal/field_test.go
new file mode 100644
index 000000000..0eda445a7
--- /dev/null
+++ b/sign/dilithium/mode3/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 := mod_q(x)
+		if y > Q {
+			t.Fatalf("mod_q(%d) > Q", x)
+		}
+		if y != x%Q {
+			t.Fatalf("mod_q(%d) != %d (mod Q)", x, x)
+		}
+	}
+}
+
+func TestReduceLe2Q(t *testing.T) {
+	for i := 0; i < 1000; i++ {
+		x := uint32(rand.Int31())
+		y := reduce_le2q(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)
+		}
+	}
+}
diff --git a/sign/dilithium/mode3/internal/kat_test.go b/sign/dilithium/mode3/internal/kat_test.go
new file mode 100644
index 000000000..5a0621d35
--- /dev/null
+++ b/sign/dilithium/mode3/internal/kat_test.go
@@ -0,0 +1,118 @@
+package internal
+
+// Code to generate the NIST "PQCsignKAT" test vectors.
+// See PQCsignKAT_sign.c and randombytes.c in the reference implementation.
+
+import (
+	"crypto/aes"
+	"crypto/sha256"
+	"fmt"
+	"testing"
+)
+
+// From SHA256SUMS in the reference implementation
+var ExpectedHashes = map[string]string{
+	"Dilithium1":     "dd83f8584fded0398547827edff081969335c32069f3e4a9dbd865fd5c2ecd2b",
+	"Dilithium2":     "532f4a7a416bba96b607395a6d07fc0eaab1f1f968e49758d2a97c718de832e7",
+	"Dilithium3":     "37a16744627f2566180a547d022f03a36d22c50080303027179751070e626c72",
+	"Dilithium4":     "4c2e6d7c8675e9345e3ab7036a4e9fb786549d242462ba9b68f58db94e84147a",
+	"Dilithium1-AES": "68fabe91565c9a664d2461c7510ac32419eadfac0566dc3e9141d276bb98e11a",
+	"Dilithium2-AES": "08865a608edcdb5723769c583b37c17c9ff8cae578f1d88df7e173ed06dd23fa",
+	"Dilithium3-AES": "f3c5fcceafa9fb2462721f272791a26c9a123b3a07fad7e07dfec232085fdd7f",
+	"Dilithium4-AES": "8de4e2ac2032f714263aa0d045275ec62b6f192f8828cfe82b63ec0b0b32deb6",
+}
+
+// See NIST's PQCgenKAT.c.
+type DRBG struct {
+	key [32]byte
+	v   [16]byte
+}
+
+func (g *DRBG) incV() {
+	for j := 15; j >= 0; j-- {
+		if g.v[j] == 255 {
+			g.v[j] = 0
+		} else {
+			g.v[j]++
+			break
+		}
+	}
+}
+
+// AES256_CTR_DRBG_Update(pd, &g.key, &g.v)
+func (g *DRBG) update(pd *[48]byte) {
+	var buf [48]byte
+	b, _ := aes.NewCipher(g.key[:])
+	for i := 0; i < 3; i++ {
+		g.incV()
+		b.Encrypt(buf[i*16:(i+1)*16], g.v[:])
+	}
+	if pd != nil {
+		for i := 0; i < 48; i++ {
+			buf[i] ^= pd[i]
+		}
+	}
+	copy(g.key[:], buf[:32])
+	copy(g.v[:], buf[32:])
+}
+
+// randombyte_init(seed, NULL, 256)
+func NewDRBG(seed *[48]byte) (g DRBG) {
+	g.update(seed)
+	return
+}
+
+// randombytes
+func (g *DRBG) Fill(x []byte) {
+	var block [16]byte
+
+	b, _ := aes.NewCipher(g.key[:])
+	for len(x) > 0 {
+		g.incV()
+		b.Encrypt(block[:], g.v[:])
+		if len(x) < 16 {
+			copy(x[:], block[:len(x)])
+			break
+		}
+		copy(x[:], block[:])
+		x = x[16:]
+	}
+	g.update(nil)
+}
+
+func TestPQCgenKATSign(t *testing.T) {
+	var seed [48]byte
+	var eseed [96]byte
+	var ppk [PublicKeySize]byte
+	var psk [PrivateKeySize]byte
+	var sig [SignatureSize]byte
+	for i := 0; i < 48; i++ {
+		seed[i] = byte(i)
+	}
+	f := sha256.New()
+	g := NewDRBG(&seed)
+	fmt.Fprintf(f, "# %s\n\n", Name)
+	for i := 0; i < 100; i++ {
+		mlen := 33 * (i + 1)
+		g.Fill(seed[:])
+		msg := make([]byte, mlen)
+		g.Fill(msg[:])
+		fmt.Fprintf(f, "count = %d\n", i)
+		fmt.Fprintf(f, "seed = %X\n", seed)
+		fmt.Fprintf(f, "mlen = %d\n", mlen)
+		fmt.Fprintf(f, "msg = %X\n", msg)
+		g2 := NewDRBG(&seed)
+		g2.Fill(eseed[:])
+		pk, sk := NewKeyFromExpandedSeed(&eseed)
+		pk.Pack(&ppk)
+		sk.Pack(&psk)
+		fmt.Fprintf(f, "pk = %X\n", ppk)
+		fmt.Fprintf(f, "sk = %X\n", psk)
+		fmt.Fprintf(f, "smlen = %d\n", mlen+SignatureSize)
+		SignTo(sk, msg[:], sig[:])
+		fmt.Fprintf(f, "sm = %X%X\n\n", sig, msg)
+	}
+	if fmt.Sprintf("%x", f.Sum(nil)) != ExpectedHashes[Name] {
+		t.Fatal()
+	}
+}
diff --git a/sign/dilithium/mode3/internal/mat.go b/sign/dilithium/mode3/internal/mat.go
new file mode 100644
index 000000000..e38a78327
--- /dev/null
+++ b/sign/dilithium/mode3/internal/mat.go
@@ -0,0 +1,29 @@
+package internal
+
+// A k by l matrix of polynomials.
+type Mat [K]VecL
+
+// Expands the given seed to a complete matrix.
+//
+// This function is called ExpandA in the specification.
+func (m *Mat) Derive(seed *[32]byte) {
+	for i := uint16(0); i < K; i++ {
+		for j := uint16(0); j < L; j++ {
+			m[i][j].DeriveUniform(seed, (i<<8)+j)
+		}
+	}
+}
+
+// Set p to the inner product of a and b using pointwise multiplication.
+//
+// Assumes a and b are in Montgomery form and their coeffients are
+// pairwise sufficiently small to add, see Poly.MulHat().  Resulting
+// coefficients are bounded by 2Lq.
+func (p *Poly) DotHat(a, b *VecL) {
+	var t Poly
+	*p = Poly{} // zero p
+	for i := 0; i < L; i++ {
+		t.MulHat(&a[i], &b[i])
+		p.Add(&t, p)
+	}
+}
diff --git a/sign/dilithium/mode3/internal/ntt.go b/sign/dilithium/mode3/internal/ntt.go
new file mode 100644
index 000000000..04b68c328
--- /dev/null
+++ b/sign/dilithium/mode3/internal/ntt.go
@@ -0,0 +1,190 @@
+package internal
+
+// Precomputed powers of the root of unity in Montgomery representation used
+// for the NTT:
+//
+//     Zetas[i] = zeta^brv(i) R mod q,
+//
+// where zeta = 1753, brv(i) is the bitreversal of a 8-bit number
+// and R=2^32 mod q.
+var Zetas [N]uint32 = [N]uint32{
+	4193792, 25847, 5771523, 7861508, 237124, 7602457, 7504169,
+	466468, 1826347, 2353451, 8021166, 6288512, 3119733, 5495562,
+	3111497, 2680103, 2725464, 1024112, 7300517, 3585928, 7830929,
+	7260833, 2619752, 6271868, 6262231, 4520680, 6980856, 5102745,
+	1757237, 8360995, 4010497, 280005, 2706023, 95776, 3077325,
+	3530437, 6718724, 4788269, 5842901, 3915439, 4519302, 5336701,
+	3574422, 5512770, 3539968, 8079950, 2348700, 7841118, 6681150,
+	6736599, 3505694, 4558682, 3507263, 6239768, 6779997, 3699596,
+	811944, 531354, 954230, 3881043, 3900724, 5823537, 2071892,
+	5582638, 4450022, 6851714, 4702672, 5339162, 6927966, 3475950,
+	2176455, 6795196, 7122806, 1939314, 4296819, 7380215, 5190273,
+	5223087, 4747489, 126922, 3412210, 7396998, 2147896, 2715295,
+	5412772, 4686924, 7969390, 5903370, 7709315, 7151892, 8357436,
+	7072248, 7998430, 1349076, 1852771, 6949987, 5037034, 264944,
+	508951, 3097992, 44288, 7280319, 904516, 3958618, 4656075,
+	8371839, 1653064, 5130689, 2389356, 8169440, 759969, 7063561,
+	189548, 4827145, 3159746, 6529015, 5971092, 8202977, 1315589,
+	1341330, 1285669, 6795489, 7567685, 6940675, 5361315, 4499357,
+	4751448, 3839961, 2091667, 3407706, 2316500, 3817976, 5037939,
+	2244091, 5933984, 4817955, 266997, 2434439, 7144689, 3513181,
+	4860065, 4621053, 7183191, 5187039, 900702, 1859098, 909542,
+	819034, 495491, 6767243, 8337157, 7857917, 7725090, 5257975,
+	2031748, 3207046, 4823422, 7855319, 7611795, 4784579, 342297,
+	286988, 5942594, 4108315, 3437287, 5038140, 1735879, 203044,
+	2842341, 2691481, 5790267, 1265009, 4055324, 1247620, 2486353,
+	1595974, 4613401, 1250494, 2635921, 4832145, 5386378, 1869119,
+	1903435, 7329447, 7047359, 1237275, 5062207, 6950192, 7929317,
+	1312455, 3306115, 6417775, 7100756, 1917081, 5834105, 7005614,
+	1500165, 777191, 2235880, 3406031, 7838005, 5548557, 6709241,
+	6533464, 5796124, 4656147, 594136, 4603424, 6366809, 2432395,
+	2454455, 8215696, 1957272, 3369112, 185531, 7173032, 5196991,
+	162844, 1616392, 3014001, 810149, 1652634, 4686184, 6581310,
+	5341501, 3523897, 3866901, 269760, 2213111, 7404533, 1717735,
+	472078, 7953734, 1723600, 6577327, 1910376, 6712985, 7276084,
+	8119771, 4546524, 5441381, 6144432, 7959518, 6094090, 183443,
+	7403526, 1612842, 4834730, 7826001, 3919660, 8332111, 7018208,
+	3937738, 1400424, 7534263, 1976782,
+}
+
+// Precomputed powers of the inverse root of unity in Montgomery
+// representation used for the inverse NTT:
+//
+//     InvZetas[i] = zeta^{brv(255-i)-256} R mod q,
+//
+// where zeta = 1753, brv(i) is the bitreversal of a 8-bit number
+// and R=2^32 mod q.
+var InvZetas [N]uint32 = [N]uint32{
+	6403635, 846154, 6979993, 4442679, 1362209, 48306, 4460757,
+	554416, 3545687, 6767575, 976891, 8196974, 2286327, 420899,
+	2235985, 2939036, 3833893, 260646, 1104333, 1667432, 6470041,
+	1803090, 6656817, 426683, 7908339, 6662682, 975884, 6167306,
+	8110657, 4513516, 4856520, 3038916, 1799107, 3694233, 6727783,
+	7570268, 5366416, 6764025, 8217573, 3183426, 1207385, 8194886,
+	5011305, 6423145, 164721, 5925962, 5948022, 2013608, 3776993,
+	7786281, 3724270, 2584293, 1846953, 1671176, 2831860, 542412,
+	4974386, 6144537, 7603226, 6880252, 1374803, 2546312, 6463336,
+	1279661, 1962642, 5074302, 7067962, 451100, 1430225, 3318210,
+	7143142, 1333058, 1050970, 6476982, 6511298, 2994039, 3548272,
+	5744496, 7129923, 3767016, 6784443, 5894064, 7132797, 4325093,
+	7115408, 2590150, 5688936, 5538076, 8177373, 6644538, 3342277,
+	4943130, 4272102, 2437823, 8093429, 8038120, 3595838, 768622,
+	525098, 3556995, 5173371, 6348669, 3122442, 655327, 522500,
+	43260, 1613174, 7884926, 7561383, 7470875, 6521319, 7479715,
+	3193378, 1197226, 3759364, 3520352, 4867236, 1235728, 5945978,
+	8113420, 3562462, 2446433, 6136326, 3342478, 4562441, 6063917,
+	4972711, 6288750, 4540456, 3628969, 3881060, 3019102, 1439742,
+	812732, 1584928, 7094748, 7039087, 7064828, 177440, 2409325,
+	1851402, 5220671, 3553272, 8190869, 1316856, 7620448, 210977,
+	5991061, 3249728, 6727353, 8578, 3724342, 4421799, 7475901,
+	1100098, 8336129, 5282425, 7871466, 8115473, 3343383, 1430430,
+	6527646, 7031341, 381987, 1308169, 22981, 1228525, 671102,
+	2477047, 411027, 3693493, 2967645, 5665122, 6232521, 983419,
+	4968207, 8253495, 3632928, 3157330, 3190144, 1000202, 4083598,
+	6441103, 1257611, 1585221, 6203962, 4904467, 1452451, 3041255,
+	3677745, 1528703, 3930395, 2797779, 6308525, 2556880, 4479693,
+	4499374, 7426187, 7849063, 7568473, 4680821, 1600420, 2140649,
+	4873154, 3821735, 4874723, 1643818, 1699267, 539299, 6031717,
+	300467, 4840449, 2867647, 4805995, 3043716, 3861115, 4464978,
+	2537516, 3592148, 1661693, 4849980, 5303092, 8284641, 5674394,
+	8100412, 4369920, 19422, 6623180, 3277672, 1399561, 3859737,
+	2118186, 2108549, 5760665, 1119584, 549488, 4794489, 1079900,
+	7356305, 5654953, 5700314, 5268920, 2884855, 5260684, 2091905,
+	359251, 6026966, 6554070, 7913949, 876248, 777960, 8143293,
+	518909, 2608894, 8354570, 4186625,
+}
+
+// Execute an in-place forward NTT on as.
+//
+// Assumes the coefficients are in Montgomery representation and bounded
+// by 2*Q.  The resulting coefficients are again in Montgomery representation,
+// but are only bounded bt 18*Q.
+func (p *Poly) NTT() {
+	// Writing z := zeta for our root of unity zeta := 1753, note z^256=-1
+	// (otherwise the order of z wouldn't be 512) and so
+	//
+	//  x^256 + 1 = x^256 - z^256
+	//            = (x^128 - z^128)(x^128 + z^128)
+	//            = (x^64 - z^64)(x^64 + z^64)(x^64 + z^192)(x^64 - z^192)
+	//          ...
+	//            = (x-z)(x+z)(x - z^129)(x + z^129) ... (x - z^255)(x + z^255)
+	//
+	// Note that the powers of z that appear (from the second line) are
+	//  in binary
+	//
+	//  01000000 11000000
+	//  00100000 10100000 01100000 11100000
+	//  00010000 10010000 01010000 11010000 00110000 10110000 01110000 11110000
+	//     ...
+	//
+	// i.e. brv(2), brv(3), brv(4), ... and these powers of z are given by
+	// the Zetas array.
+	//
+	// The polynomials x ± z^i are irreducable and coprime, hence by the
+	// Chinese Remainder Theorem we know
+	//
+	//  R[x]/(x^256+1) ---> R[x] / (x-z) x ... x R[x] / (x+z^255)
+	//                      ~= ∏_i R
+	//
+	// given by
+	//
+	//  a |---> = ( a mod x-z, ..., a mod x+z^255 )
+	//          ~ ( a(z), a(-z), a(z^129), a(-z^129), ..., a(z^255), a(-z^255) )
+	//
+	// is an isomorphism, which is the forward NTT.  It can be computed
+	// efficiently by computing
+	//
+	//  a |---> ( a mod x^128 - z^128, a mod x^128 + z^128 )
+	//    |---> ( a mod x^64 - z^64, a mod x^64 + z^64,
+	//            a mod x^64 - z^192, a mod x^64 + z^192 )
+	//       et cetera
+
+	k := 0 // Index into Zetas
+
+	for l := uint(N / 2); l > 0; l >>= 1 {
+		// On the n-th iteration of the l-loop, the coefficients start off
+		// bounded by n*2*Q.
+		for offset := uint(0); offset < N-l; offset += 2 * l {
+			k++
+			zeta := uint64(Zetas[k])
+			for j := offset; j < offset+l; j++ {
+				t := mont_reduce_le2q(zeta * uint64(p[j+l]))
+				p[j+l] = p[j] + (2*Q - t) // Cooley--Tukey butterfly
+				p[j] += t
+			}
+		}
+	}
+}
+
+// Execute an in-place inverse NTT and multiply by Montgomery factor R
+//
+// Assumes the coefficients are in Montgomery representation and bounded
+// by 2*Q.  The resulting coefficients are again in Montgomery representation
+// and bounded by 2*Q.
+func (p *Poly) InvNTT() {
+	k := 0 // Index into InvZetas
+
+	// We basically do the opposite of NTT, but postpone dividing by 2 in the
+	// inverse of the Cooley--Tukey butterfly and accumelate that to a big
+	// division by 2^8 at the end.
+
+	for l := uint(1); l < N; l <<= 1 {
+		// On the n-th iteration of the l-loop, the coefficients start off
+		// bounded by 2^(n-1)*2*Q, so by 256*Q on the last.
+		for offset := uint(0); offset < N-l; offset += 2 * l {
+			zeta := uint64(InvZetas[k])
+			k++
+			for j := offset; j < offset+l; j++ {
+				t := p[j] // Gentleman--Sande butterfly
+				p[j] = t + p[j+l]
+				t += 256*Q - p[j+l]
+				p[j+l] = mont_reduce_le2q(zeta * uint64(t))
+			}
+		}
+	}
+
+	for j := uint(0); j < N; j++ {
+		// 41978 = (256)^-1 R^2
+		p[j] = mont_reduce_le2q(41978 * uint64(p[j]))
+	}
+}
diff --git a/sign/dilithium/mode3/internal/ntt_test.go b/sign/dilithium/mode3/internal/ntt_test.go
new file mode 100644
index 000000000..fc8f6e9c4
--- /dev/null
+++ b/sign/dilithium/mode3/internal/ntt_test.go
@@ -0,0 +1,54 @@
+package internal
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func (p *Poly) RandLe2Q() {
+	for i := uint(0); i < N; i++ {
+		p[i] = uint32(rand.Intn(int(2 * Q)))
+	}
+}
+
+func TestNTT(t *testing.T) {
+	for k := 0; k < 1000; k++ {
+		var p, q Poly
+		p.RandLe2Q()
+		q = p
+		q.Normalize()
+		p.NTT()
+		for i := uint(0); i < N; i++ {
+			if p[i] > 16*Q {
+				t.Fatalf("NTT(%v)[%d] = %d > 18*Q", q, i, p[i])
+			}
+		}
+		p.ReduceLe2Q()
+		p.InvNTT()
+		for i := uint(0); i < N; i++ {
+			if p[i] > 2*Q {
+				t.Fatalf("InvNTT(%v)[%d] > 2*Q", q, i)
+			}
+		}
+		p.Normalize()
+		for i := uint(0); i < N; i++ {
+			if p[i] != uint32((uint64(q[i])*uint64(1<<32))%Q) {
+				t.Fatalf("%v != %v", p, q)
+			}
+		}
+	}
+}
+
+func BenchmarkNTT(b *testing.B) {
+	var p Poly
+	for i := 0; i < b.N; i++ {
+		p.NTT()
+	}
+}
+
+func BenchmarkInvNTT(b *testing.B) {
+	var p Poly
+	for i := 0; i < b.N; i++ {
+		p.InvNTT()
+	}
+}
diff --git a/sign/dilithium/mode3/internal/pack.go b/sign/dilithium/mode3/internal/pack.go
new file mode 100644
index 000000000..acd39c61a
--- /dev/null
+++ b/sign/dilithium/mode3/internal/pack.go
@@ -0,0 +1,333 @@
+package internal
+
+// Writes p with norm less than or equal η into buf, which must be of
+// size PolyLeqEtaSize.
+//
+// Assumes coefficients of p are not normalized, but in [q-η,q+η].
+func (p *Poly) PackLeqEta(buf []byte) {
+	// TODO check if Go compiler eliminates branch
+	if DoubleEtaBits == 4 {
+		j := 0
+		for i := 0; i < PolyLeqEtaSize; i++ {
+			buf[i] = (byte(Q+Eta-p[j]) |
+				byte(Q+Eta-p[j+1])<<4)
+			j += 2
+		}
+	} else if DoubleEtaBits == 3 {
+		j := 0
+		for i := 0; i < PolyLeqEtaSize; i += 3 {
+			buf[i] = (byte(Q+Eta-p[j]) |
+				(byte(Q+Eta-p[j+1]) << 3) |
+				(byte(Q+Eta-p[j+2]) << 6))
+			buf[i+1] = ((byte(Q+Eta-p[j+2]) >> 2) |
+				(byte(Q+Eta-p[j+3]) << 1) |
+				(byte(Q+Eta-p[j+4]) << 4) |
+				(byte(Q+Eta-p[j+5]) << 7))
+			buf[i+2] = ((byte(Q+Eta-p[j+5]) >> 1) |
+				(byte(Q+Eta-p[j+6]) << 2) |
+				(byte(Q+Eta-p[j+7]) << 5))
+			j += 8
+		}
+
+	} else {
+		panic("eta not supported")
+	}
+}
+
+// Sets p to the polynomial of norm less than or equal η encoded in the
+// given buffer of size PolyLeqEtaSize.
+//
+// Output coefficients of p are not normalized, but in [q-η,q+η] provided
+// buf was created using PackLeqEta.
+//
+// Beware, for arbitrary buf the coefficients of p might en up in
+// the interval [q-2^b,q+2^b] where b is the least b with η≤2^b.
+func (p *Poly) UnpackLeqEta(buf []byte) {
+	// TODO check if branch is optimized out
+	if DoubleEtaBits == 4 {
+		j := 0
+		for i := 0; i < PolyLeqEtaSize; i++ {
+			p[j] = Q + Eta - uint32(buf[i]&15)
+			p[j+1] = Q + Eta - uint32(buf[i]>>4)
+			j += 2
+		}
+	} else if DoubleEtaBits == 3 {
+		j := 0
+		for i := 0; i < PolyLeqEtaSize; i += 3 {
+			p[j] = Q + Eta - uint32(buf[i]&7)
+			p[j+1] = Q + Eta - uint32((buf[i]>>3)&7)
+			p[j+2] = Q + Eta - uint32((buf[i]>>6)|((buf[i+1]<<2)&7))
+			p[j+3] = Q + Eta - uint32((buf[i+1]>>1)&7)
+			p[j+4] = Q + Eta - uint32((buf[i+1]>>4)&7)
+			p[j+5] = Q + Eta - uint32((buf[i+1]>>7)|((buf[i+2]<<1)&7))
+			p[j+6] = Q + Eta - uint32((buf[i+2]>>2)&7)
+			p[j+7] = Q + Eta - uint32((buf[i+2]>>5)&7)
+			j += 8
+		}
+	} else {
+		panic("eta not supported")
+	}
+}
+
+// Writes p whose coefficients are less than 512 into buf, which must be
+// of size at least PolyT1Size .
+//
+// Assumes coefficients of p are normalized.
+func (p *Poly) PackT1(buf []byte) {
+	j := 0
+	for i := 0; i < PolyT1Size; i += 9 {
+		buf[i] = byte(p[j])
+		buf[i+1] = byte(p[j]>>8) | byte(p[j+1]<<1)
+		buf[i+2] = byte(p[j+1]>>7) | byte(p[j+2]<<2)
+		buf[i+3] = byte(p[j+2]>>6) | byte(p[j+3]<<3)
+		buf[i+4] = byte(p[j+3]>>5) | byte(p[j+4]<<4)
+		buf[i+5] = byte(p[j+4]>>4) | byte(p[j+5]<<5)
+		buf[i+6] = byte(p[j+5]>>3) | byte(p[j+6]<<6)
+		buf[i+7] = byte(p[j+6]>>2) | byte(p[j+7]<<7)
+		buf[i+8] = byte(p[j+7] >> 1)
+		j += 8
+	}
+}
+
+// Sets p to the polynomial whose coefficients are less than 512 encoded
+// into buf (which must be of size PolyT1Size).
+//
+// p will be normalized.
+func (p *Poly) UnpackT1(buf []byte) {
+	j := 0
+	for i := 0; i < PolyT1Size; i += 9 {
+		p[j] = (uint32(buf[i]) | (uint32(buf[i+1]) << 8)) & 0x1ff
+		p[j+1] = (uint32(buf[i+1]>>1) | (uint32(buf[i+2]) << 7)) & 0x1ff
+		p[j+2] = (uint32(buf[i+2]>>2) | (uint32(buf[i+3]) << 6)) & 0x1ff
+		p[j+3] = (uint32(buf[i+3]>>3) | (uint32(buf[i+4]) << 5)) & 0x1ff
+		p[j+4] = (uint32(buf[i+4]>>4) | (uint32(buf[i+5]) << 4)) & 0x1ff
+		p[j+5] = (uint32(buf[i+5]>>5) | (uint32(buf[i+6]) << 3)) & 0x1ff
+		p[j+6] = (uint32(buf[i+6]>>6) | (uint32(buf[i+7]) << 2)) & 0x1ff
+		p[j+7] = (uint32(buf[i+7]>>7) | (uint32(buf[i+8]) << 1)) & 0x1ff
+		j += 8
+	}
+}
+
+// Writes p whose coefficients are in (-2^{d-1}, 2^{d-1}] into buf which
+// has to be of length at least PolyT0Size.
+//
+// Assumes that the coefficients are not normalized, but lie in the
+// range (q-2^{d-1}, q+2^{d-1}].
+func (p *Poly) PackT0(buf []byte) {
+	j := 0
+	for i := 0; i < PolyT0Size; i += 7 {
+		// TODO constant for 1<<(D-1)
+		p0 := Q + (1 << (D - 1)) - p[j]
+		p1 := Q + (1 << (D - 1)) - p[j+1]
+		p2 := Q + (1 << (D - 1)) - p[j+2]
+		p3 := Q + (1 << (D - 1)) - p[j+3]
+
+		buf[i] = byte(p0)
+		buf[i+1] = byte(p0>>8) | byte(p1<<6)
+		buf[i+2] = byte(p1 >> 2)
+		buf[i+3] = byte(p1>>10) | byte(p2<<4)
+		buf[i+4] = byte(p2 >> 4)
+		buf[i+5] = byte(p2>>12) | byte(p3<<2)
+		buf[i+6] = byte(p3 >> 6)
+		j += 4
+	}
+}
+
+// Sets p to the polynomial packed into buf by PackT0.
+//
+// The coefficients of p will not be normalized, but will lie
+// in (-2^{d-1}, 2^{d-1}].
+func (p *Poly) UnpackT0(buf []byte) {
+	j := 0
+	for i := 0; i < PolyT0Size; i += 7 {
+		p[j] = Q + (1 << (D - 1)) - ((uint32(buf[i]) |
+			(uint32(buf[i+1]) << 8)) & 0x3fff)
+		p[j+1] = Q + (1 << (D - 1)) - (((uint32(buf[i+1]) >> 6) |
+			(uint32(buf[i+2]) << 2) |
+			(uint32(buf[i+3]) << 10)) & 0x3fff)
+		p[j+2] = Q + (1 << (D - 1)) - (((uint32(buf[i+3]) >> 4) |
+			(uint32(buf[i+4]) << 4) |
+			(uint32(buf[i+5]) << 12)) & 0x3fff)
+		p[j+3] = Q + (1 << (D - 1)) - ((uint32(buf[i+5]) >> 2) |
+			(uint32(buf[i+6]) << 6))
+		j += 4
+	}
+}
+
+// Writes p whose coefficients are in norm less than γ1 into buf
+// which has to be of length PolyLeGamma1Size.
+//
+// Assumes p is normalized.
+func (p *Poly) PackLeGamma1(buf []byte) {
+	j := 0
+	for i := 0; i < PolyLeGamma1Size; i += 5 {
+		// Coefficients are in [0, γ1) ∪ (Q-γ1, Q)
+		p0 := Gamma1 - 1 - p[j]         // ... in [0, γ1) ∪ [γ1-1-Q, 2(γ1-1)-Q]
+		p0 += uint32(int32(p0)>>31) & Q // ... in [0, 2(γ1-1)]
+		p1 := Gamma1 - 1 - p[j+1]
+		p1 += uint32(int32(p1)>>31) & Q
+
+		buf[i] = byte(p0)
+		buf[i+1] = byte(p0 >> 8)
+		buf[i+2] = byte(p0>>16) | byte(p1<<4)
+		buf[i+3] = byte(p1 >> 4)
+		buf[i+4] = byte(p1 >> 12)
+		j += 2
+	}
+}
+
+// Sets p to the polynomial packed into buf by PackLeGamma1.
+//
+// p will be normalized.
+//
+// Beware, for arbitrary buf the coefficients of p might exceed γ1.
+func (p *Poly) UnpackLeGamma1(buf []byte) {
+	j := 0
+	for i := 0; i < PolyLeGamma1Size; i += 5 {
+		p0 := Gamma1 - 1 - ((uint32(buf[i]) |
+			(uint32(buf[i+1]) << 8) |
+			(uint32(buf[i+2]) << 16)) & 0xfffff)
+		p1 := Gamma1 - 1 - ((uint32(buf[i+2]) >> 4) |
+			(uint32(buf[i+3]) << 4) |
+			(uint32(buf[i+4]) << 12))
+		p0 += uint32(int32(p0)>>31) & Q
+		p1 += uint32(int32(p1)>>31) & Q
+		p[j] = p0
+		p[j+1] = p1
+		j += 2
+	}
+}
+
+// Writes p whose coefficients are in [0, 16) to buf, which must be of
+// length N/2.
+func (p *Poly) PackLe16(buf []byte) {
+	j := 0
+	for i := 0; i < PolyLe16Size; i++ {
+		buf[i] = byte(p[j]) | byte(p[j+1]<<4)
+		j += 2
+	}
+}
+
+// Writes v with coefficients in {0, 1} of which at most ω non-zero
+// to buf, which must have length ω+k.
+func (v *VecK) PackHint(buf []byte) {
+	// TODO optimize?
+
+	// The packed hint starts with the indices of the non-zero coefficients
+	// For instance:
+	//
+	//    (x^56 + x^100, x^255, 0, x^2 + x^23, x^1)
+	//
+	// Yields
+	//
+	//  56, 100, 255, 2, 23, 1
+	//
+	// Then we pad with zeroes until we have a list of ω items:
+	// //  56, 100, 255, 2, 23, 1, 0, 0, ..., 0
+	//
+	// Then we finish with a list of the switch-over-indices in this
+	// list between polynomials, so:
+	//
+	//  56, 100, 255, 2, 23, 1, 0, 0, ..., 0, 2, 3, 3, 5, 6
+
+	off := uint8(0)
+	for i := 0; i < K; i++ {
+		for j := uint16(0); j < N; j++ {
+			if v[i][j] != 0 {
+				buf[off] = uint8(j)
+				off += 1
+			}
+		}
+		buf[Omega+i] = off
+	}
+	for ; off < Omega; off++ {
+		buf[off] = 0
+	}
+}
+
+// Sets v to the vector encoded using VecK.PackHint()
+//
+// Returns whether unpacking was successful.
+func (v *VecK) UnpackHint(buf []byte) bool {
+	// A priori, there would be several reasonable ways to encode the same
+	// hint vector.  We take care to only allow only one encoding, to ensure
+	// "strong unforgeability".
+	//
+	// See PackHint() source for description of the encoding.
+	*v = VecK{}         // zero v
+	prevSOP := uint8(0) // previous switch-over-point
+	for i := 0; i < K; i++ {
+		SOP := buf[Omega+i]
+		if SOP < prevSOP || SOP > Omega {
+			return false // ensures switch-over-points are increasing
+		}
+		for j := prevSOP; j < SOP; j++ {
+			if j > prevSOP && buf[j] <= buf[j-1] {
+				return false // ensures indices are increasing (within a poly)
+			}
+			v[i][buf[j]] = 1
+		}
+		prevSOP = SOP
+	}
+	for j := prevSOP; j < Omega; j++ {
+		if buf[j] != 0 {
+			return false // ensures padding indices are zero
+		}
+	}
+
+	return true
+}
+
+// Writes p with 60 non-zero coefficients {-1,1} to buf, which must have
+// length 40.
+func (p *Poly) PackB60(buf []byte) {
+	// TODO optimize?
+
+	// We start with a mask of the non-zero positions of p (which is 32 bytes)
+	// and then append 60 packed bits, where a one indicates a negative
+	// coefficients.
+	var signs uint64
+	mask := uint64(1)
+	for i := 0; i < 32; i++ {
+		buf[i] = 0
+		for j := 0; j < 8; j++ {
+			if p[8*i+j] != 0 {
+				buf[i] |= 1 << uint(j)
+				if p[8*i+j] == Q-1 {
+					signs |= mask
+				}
+				mask <<= 1
+			}
+		}
+	}
+	for i := uint64(0); i < 8; i++ {
+		buf[i+32] = uint8(signs >> (8 * i))
+	}
+}
+
+// Sets p to the polynomial packed into buf iwth Poly.PackB60().
+//
+// Returns whether unpacking was successful.
+func (p *Poly) UnpackB60(buf []byte) bool {
+	*p = Poly{} // zero p
+	signs := (uint64(buf[32]) | (uint64(buf[33]) << 8) |
+		(uint64(buf[34]) << 16) | (uint64(buf[35]) << 24) |
+		(uint64(buf[36]) << 32) | (uint64(buf[37]) << 40) |
+		(uint64(buf[38]) << 48) | (uint64(buf[39]) << 56))
+	if signs>>60 != 0 {
+		return false // ensure unusud bits are zero for strong unforgeability
+	}
+
+	for i := 0; i < 32; i++ {
+		for j := 0; j < 8; j++ {
+			if (buf[i]>>uint(j))&1 == 1 {
+				p[8*i+j] = 1
+				// Note 1 ^ (1 | (Q-1)) = Q-1 and (-1)&x = x
+				p[8*i+j] ^= uint32(-(signs & 1)) & (1 | (Q - 1))
+				signs >>= 1
+			}
+		}
+	}
+
+	return true
+}
diff --git a/sign/dilithium/mode3/internal/pack_test.go b/sign/dilithium/mode3/internal/pack_test.go
new file mode 100644
index 000000000..09888bbf9
--- /dev/null
+++ b/sign/dilithium/mode3/internal/pack_test.go
@@ -0,0 +1,83 @@
+package internal
+
+import (
+	"testing"
+)
+
+func TestPolyPackLeqEta(t *testing.T) {
+	var p1, p2 Poly
+	var seed [32]byte
+	var buf [PolyLeqEtaSize]byte
+
+	for i := uint16(0); i < 100; i++ {
+		// Note that DeriveUniformLeqEta sets p to the right kind of
+		// unnormalized vector.
+		p1.DeriveUniformLeqEta(&seed, i)
+		for j := 0; j < PolyLeqEtaSize; j++ {
+			if p1[j] < Q-Eta || p1[j] > Q+Eta {
+				t.Fatalf("DerveUniformLeqEta out of bounds")
+			}
+		}
+		p1.PackLeqEta(buf[:])
+		p2.UnpackLeqEta(buf[:])
+		if p1 != p2 {
+			t.Fatalf("%v != %v", p1, p2)
+		}
+	}
+}
+
+func TestPolyPackT1(t *testing.T) {
+	var p1, p2 Poly
+	var seed [32]byte
+	var buf [PolyT1Size]byte
+
+	for i := uint16(0); i < 100; i++ {
+		p1.DeriveUniform(&seed, i)
+		p1.Normalize()
+		for j := 0; j < N; j++ {
+			p1[j] &= 0x1ff
+		}
+		p1.PackT1(buf[:])
+		p2.UnpackT1(buf[:])
+		if p1 != p2 {
+			t.Fatalf("%v != %v", p1, p2)
+		}
+	}
+}
+
+func TestPolyPackT0(t *testing.T) {
+	var p, p0, p1, p2 Poly
+	var seed [32]byte
+	var buf [PolyT0Size]byte
+
+	for i := uint16(0); i < 100; i++ {
+		p.DeriveUniform(&seed, i)
+		p.Normalize()
+		p.Power2Round(&p0, &p1)
+
+		p0.PackT0(buf[:])
+		p2.UnpackT0(buf[:])
+		if p0 != p2 {
+			t.Fatalf("%v != %v", p1, p2)
+		}
+	}
+}
+
+func TestPolyPackLeGamma1(t *testing.T) {
+	var p0, p1 Poly
+	var seed [48]byte
+	var buf [PolyLeGamma1Size]byte
+
+	for i := uint16(0); i < 100; i++ {
+		p0.DeriveUniformLeGamma1(&seed, i)
+		p0.Normalize()
+
+		p0.PackLeGamma1(buf[:])
+		p1.UnpackLeGamma1(buf[:])
+		if p0 != p1 {
+			t.Fatalf("%v != %v", p0, p1)
+		}
+	}
+}
+
+// TODO Test (Un)Hint
diff --git a/sign/dilithium/mode3/internal/params.go b/sign/dilithium/mode3/internal/params.go
new file mode 100644
index 000000000..880857f81
--- /dev/null
+++ b/sign/dilithium/mode3/internal/params.go
@@ -0,0 +1,14 @@
+package internal
+
+const (
+	Name           = "Dilithium3"
+	PublicKeySize  = 1472
+	PrivateKeySize = 3504
+	SignatureSize  = 2701
+	K              = 5
+	L              = 4
+	Eta            = 5
+	DoubleEtaBits  = 4
+	Beta           = 275
+	Omega          = 96
+)
diff --git a/sign/dilithium/mode3/internal/params_test.go b/sign/dilithium/mode3/internal/params_test.go
new file mode 100644
index 000000000..1c9069bce
--- /dev/null
+++ b/sign/dilithium/mode3/internal/params_test.go
@@ -0,0 +1,84 @@
+package internal
+
+import (
+	"testing"
+
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+// Tests specific to the current mode
+
+func hash(in []byte) [16]byte {
+	var ret [16]byte
+	h := shake.NewShake256()
+	h.Write(in)
+	h.Read(ret[:])
+	return ret
+}
+
+func TestNewKeyFromSeed(t *testing.T) {
+	var seed [SeedSize]byte
+	var pkp [PublicKeySize]byte
+	var skp [PrivateKeySize]byte
+	pk, sk := NewKeyFromSeed(&seed)
+	pk.Pack(&pkp)
+	sk.Pack(&skp)
+
+	// Generated with reference implementation.
+	ehpk := [16]byte{
+		183, 37, 211, 31, 183, 9, 102, 79, 133, 135,
+		226, 251, 106, 96, 254, 128,
+	}
+	ehsk := [16]byte{
+		164, 79, 207, 31, 67, 209, 36, 134, 92, 99,
+		203, 243, 129, 163, 183, 235,
+	}
+
+	if hash(pkp[:]) != ehpk {
+		t.Fatalf("pk not ok")
+	}
+	if hash(skp[:]) != ehsk {
+		t.Fatalf("sk not ok")
+	}
+}
+
+func TestVectorDeriveUniformLeqEta(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	p2 := Poly{
+		1, 8380415, 0, 1, 8380412, 0, 0, 8380413, 8380414, 8380414,
+		8380416, 5, 2, 5, 0, 8380412, 8380415, 1, 2, 8380412, 5,
+		1, 1, 8380412, 5, 8380413, 8380416, 8380414, 5, 8380412,
+		0, 2, 5, 8380414, 2, 3, 3, 0, 8380416, 5, 1, 2, 8380414,
+		8380413, 0, 0, 5, 3, 8380415, 8380414, 5, 3, 3, 2, 8380416,
+		2, 0, 8380412, 8380414, 1, 2, 8380412, 8380412, 5, 8380413,
+		4, 8380413, 8380415, 8380414, 4, 8380413, 5, 8380413, 1,
+		8380416, 8380415, 8380412, 2, 3, 4, 8380415, 8380416, 4,
+		8380415, 8380413, 1, 0, 1, 4, 8380415, 8380415, 4, 8380412,
+		8380415, 4, 3, 5, 0, 8380415, 8380415, 2, 8380415, 3,
+		8380414, 8380416, 3, 5, 5, 0, 8380415, 4, 3, 2, 8380415,
+		8380415, 5, 8380412, 8380414, 2, 8380415, 8380412, 1,
+		8380412, 8380413, 5, 8380412, 3, 4, 4, 8380413, 2, 8380415,
+		8380415, 8380414, 8380413, 8380413, 8380413, 2, 8380413,
+		4, 8380416, 2, 8380416, 0, 8380415, 8380414, 2, 8380412,
+		8380414, 4, 3, 8380416, 1, 8380412, 2, 8380412, 8380412,
+		1, 0, 0, 0, 8380415, 8380414, 8380413, 8380412, 8380415,
+		1, 1, 5, 8380412, 8380416, 2, 0, 8380413, 0, 2, 2, 2,
+		8380415, 1, 3, 8380412, 8380413, 2, 8380416, 5, 8380414,
+		8380415, 4, 4, 8380416, 5, 8380414, 8380412, 8380413, 4,
+		5, 2, 5, 2, 4, 4, 8380414, 5, 3, 8380415, 2, 1, 5, 8380413,
+		1, 8380414, 5, 8380412, 8380414, 0, 1, 8380416, 5, 8380412,
+		5, 8380413, 8380415, 8380414, 8380413, 0, 5, 1, 8380412,
+		8380415, 8380412, 8380416, 5, 2, 0, 3, 0, 4, 0, 8380416,
+		0, 8380413, 2, 0, 8380416, 8380416, 8380412, 3, 8380412,
+		1, 2, 8380414, 3, 4, 3, 0,
+	}
+	for i := 0; i < 32; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniformLeqEta(&seed, 30000)
+	p.Normalize()
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
diff --git a/sign/dilithium/mode3/internal/poly.go b/sign/dilithium/mode3/internal/poly.go
new file mode 100644
index 000000000..061a931f9
--- /dev/null
+++ b/sign/dilithium/mode3/internal/poly.go
@@ -0,0 +1,131 @@
+package internal
+
+// An element of our base ring R which are polynomials over Z_q modulo
+// the equation X^N = -1, where q=2^23 - 2^13 + 1 and N=256.
+//
+// Coefficients aren't always reduced.  See Freeze()
+type Poly [N]uint32
+
+// Reduces each of the coefficients to <2q.
+func (p *Poly) ReduceLe2Q() {
+	for i := uint(0); i < N; i++ {
+		p[i] = reduce_le2q(p[i])
+	}
+}
+
+// Reduce each of the coefficients to <q.
+func (p *Poly) Normalize() {
+	for i := uint(0); i < N; i++ {
+		p[i] = mod_q(p[i])
+	}
+}
+
+// Normalize the coefficients in this polynomial assuming they are already
+// bounded by 2q.
+func (p *Poly) NormalizeAssumingLe2Q() {
+	for i := 0; i < N; i++ {
+		p[i] = le2q_mod_q(p[i])
+	}
+}
+
+// Sets p to a + b.  Does not normalize polynomials.
+func (p *Poly) Add(a, b *Poly) {
+	for i := uint(0); i < N; i++ {
+		p[i] = a[i] + b[i]
+	}
+}
+
+// Sets p to a - b.
+//
+// Warning: assumes coefficients of b are less than 2q.
+func (p *Poly) Sub(a, b *Poly) {
+	for i := uint(0); i < N; i++ {
+		p[i] = a[i] + (2*Q - b[i])
+	}
+}
+
+// Checks whether the "supnorm" (see sec 2.1 of the spec) of p is equal
+// or greater than the given bound.
+//
+// Requires the coefficients of p to be normalized.
+func (p *Poly) Exceeds(bound uint32) bool {
+	// Note that we are allowed to leak which coefficients breaks the bound,
+	// but not its sign.
+	for i := 0; i < N; i++ {
+		// The central. reps. of {0,       1, ..., (Q-1)/2,  (Q+1)/2, ..., Q-1}
+		// are given by          {0,       1, ..., (Q-1)/2, -(Q-1)/2, ...,  -1}
+		// so their norms are    {0,       1, ..., (Q-1)/2l, (Q-1)/2, ...,   1}.
+		// We'll compute them in a different way though.
+
+		// Sets x to             {(Q-1)/2, (Q-3)/2, ..., 0, -1, ..., -(Q-1)/2}
+		x := int32((Q-1)/2) - int32(p[i])
+		// Sets x to             {(Q-1)/2, (Q-3)/2, ..., 0, 0, ...,  (Q-3)/2}
+		x ^= (x >> 31)
+		// Sets x to             {0,       1, ...,  (Q-1)/2, (Q-1)/2, ..., 1}
+		x = int32((Q-1)/2) - x
+		if uint32(x) >= bound {
+			return true
+		}
+	}
+	return false
+}
+
+// Splits each of the coefficients using decompose.
+//
+// Requires p to be normalized.
+func (p *Poly) Decompose(p0PlusQ, p1 *Poly) {
+	for i := 0; i < N; i++ {
+		p0PlusQ[i], p1[i] = decompose(p[i])
+	}
+}
+
+// Splits p into p1 and p0 such that [i]p1 * 2^D + [i]p0 = [i]p
+// with -2^{D-1} < [i]p0 ≤ 2^{D-1}.  Returns p0 + Q and p1.
+//
+// Requires the coefficients of p to be normalized.
+func (p *Poly) Power2Round(p0PlusQ, p1 *Poly) {
+	for i := 0; i < N; i++ {
+		p0PlusQ[i], p1[i] = power2round(p[i])
+	}
+}
+
+// Sets p to the hint polynomail for low part p0 and high part p1.
+//
+// Returns the number of ones in the hint vector.
+func (p *Poly) MakeHint(p0, p1 *Poly) (pop uint32) {
+	for i := 0; i < N; i++ {
+		h := make_hint(p0[i], p1[i])
+		pop += h
+		p[i] = h
+	}
+	return
+}
+
+// Computes corrections to the high bits of the polynomial q according
+// to the hints in h and sets p to the corrected high bits.  Returns p.
+func (p *Poly) UseHint(q, hint *Poly) *Poly {
+	for i := 0; i < N; i++ {
+		p[i] = use_hint(q[i], hint[i])
+	}
+	return p
+}
+
+// Sets p to the polynomial whose coefficients are the pointwise multiplication
+// of those of a and b.  The coefficients of p are bounded by 2q.
+//
+// Assumes a and b are in Montgomery form and that the pointwise product
+// of each coefficient is below 2^32 q.
+func (p *Poly) MulHat(a, b *Poly) {
+	for i := 0; i < N; i++ {
+		p[i] = mont_reduce_le2q(uint64(a[i]) * uint64(b[i]))
+	}
+}
+
+// Sets p to 2^d q without reducing.
+//
+// So it requires the coefficients of p  to be less than 2^{32-D}.
+func (p *Poly) MulBy2toD(q *Poly) {
+	for i := 0; i < N; i++ {
+		p[i] = q[i] << D
+	}
+}
diff --git a/sign/dilithium/mode3/internal/sample.go b/sign/dilithium/mode3/internal/sample.go
new file mode 100644
index 000000000..9e61a0d63
--- /dev/null
+++ b/sign/dilithium/mode3/internal/sample.go
@@ -0,0 +1,186 @@
+package internal
+
+import (
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+// Sample p uniformly from the given seed and nonce.
+//
+// p will be normalized.
+func (p *Poly) DeriveUniform(seed *[32]byte, nonce uint16) {
+	// TODO support AES
+	var iv [32 + 2]byte // 32 byte seed + uint16 nonce
+	var buf [168]byte   // 168 SHAKE-128 rate
+	var i int
+
+	h := shake.NewShake128()
+	copy(iv[:32], seed[:])
+	iv[32] = uint8(nonce)
+	iv[33] = uint8(nonce >> 8)
+	h.Write(iv[:])
+
+	for i < N {
+		h.Read(buf[:])
+
+		// Note that 3 divides into 168, so we use up buf completely.
+		for j := 0; j < 168 && i < N; j += 3 {
+			t := (uint32(buf[j]) | (uint32(buf[j+1]) << 8) |
+				(uint32(buf[j+2]) << 16)) & 0x7fffff
+
+			// We use rejection sampling
+			if t < Q {
+				p[i] = t
+				i++
+			}
+		}
+	}
+}
+
+// Sample p uniformly with coefficients of norm less than or equal η,
+// using the given seed and nonce.
+//
+// p will not be normalized, but will have coefficients in [q-η,q+η].
+func (p *Poly) DeriveUniformLeqEta(seed *[32]byte, nonce uint16) {
+	// Assumes 2 < η < 8.
+	// TODO support AES
+	var iv [32 + 2]byte // 32 byte seed + uint16 nonce
+	var buf [168]byte   // 168 SHAKE-128 rate
+	var i int
+
+	h := shake.NewShake128()
+	copy(iv[:32], seed[:])
+	iv[32] = uint8(nonce)
+	iv[33] = uint8(nonce >> 8)
+	h.Write(iv[:])
+
+	for i < N {
+		h.Read(buf[:])
+
+		// We use rejection sampling
+		for j := 0; j < 168 && i < N; j++ {
+			// TODO check whether the Go compiler gets rid of this conditional
+			//      (probably not.)
+			var t1, t2 uint32
+			if Eta <= 3 {
+				t1 = uint32(buf[j]) & 7
+				t2 = uint32(buf[j]) >> 5
+			} else {
+				t1 = uint32(buf[j]) & 15
+				t2 = uint32(buf[j]) >> 4
+			}
+			if t1 <= 2*Eta {
+				p[i] = Q + Eta - t1
+				i++
+			}
+			if t2 <= 2*Eta && i < N {
+				p[i] = Q + Eta - t2
+				i++
+			}
+		}
+	}
+}
+
+// Sample p uniformly with coefficients of norm less than γ1 using the
+// given seed and nonce.
+//
+// p will not be normalized, but have coefficients in the
+// interval (q-γ1,q+γ1)
+func (p *Poly) DeriveUniformLeGamma1(seed *[48]byte, nonce uint16) {
+	// Assumes γ1 is less than 2^20.
+	// TODO support AES
+	var iv [48 + 2]byte // 48 byte seed + uint16 nonce
+	// SHAKE-256 rate is 136 and we take 5 bytes at a same time.  As 136 is
+	// 1 modulo 5, we are left with 1 byte after the first block, which we
+	// include in the next block.  So we need 4 bytes leeway in our buffer.
+	var buf [136 + 4]byte
+	var i int
+	bufOffset := 0 // where to put the next block
+
+	h := shake.NewShake256()
+	copy(iv[:48], seed[:])
+	iv[48] = uint8(nonce)
+	iv[49] = uint8(nonce >> 8)
+	h.Write(iv[:])
+
+	for i < N {
+		h.Read(buf[bufOffset : bufOffset+136])
+
+		// We use rejection sampling
+		for j := 0; j < 136-4 && i < N; j += 5 {
+			t1 := (uint32(buf[j]) | (uint32(buf[j+1]) << 8) |
+				(uint32(buf[j+2]) << 16)) & 0xfffff
+			t2 := ((uint32(buf[j+2]) >> 4) | (uint32(buf[j+3]) << 4) |
+				(uint32(buf[j+4]) << 12))
+
+			if t1 <= 2*Gamma1-2 {
+				p[i] = Q + Gamma1 - 1 - t1
+				i++
+			}
+			if t2 <= 2*Gamma1-2 && i < N {
+				p[i] = Q + Gamma1 - 1 - t2
+				i++
+			}
+		}
+
+		bufOffset++
+		if bufOffset == 5 {
+			bufOffset = 0
+		}
+
+		// Move remaining bytes at the end to the start.
+		for j := 0; j < bufOffset; j++ {
+			buf[j] = buf[135+j]
+		}
+	}
+}
+
+// Samples p uniformly with 60 non-zero coefficients in {q-1,1}.
+//
+// The polynomial p will be normalized.
+func (p *Poly) DeriveUniformB60(seed *[48]byte, w1 *VecK) {
+	var w1Packed [PolyLe16Size * K]byte
+	var buf [136]byte // SHAKE-256 rate is 136
+
+	w1.PackLe16(w1Packed[:])
+
+	h := shake.NewShake256()
+	h.Write(seed[:])
+	h.Write(w1Packed[:])
+	h.Read(buf[:])
+
+	// Essentially we generate a sequence of 60 ones or minus ones,
+	// prepend 196 zeroes and shuffle the concatenation using the
+	// usual algorithm (Fisher--Yates.)
+	signs := (uint64(buf[0]) | (uint64(buf[1]) << 8) | (uint64(buf[2]) << 16) |
+		(uint64(buf[3]) << 24) | (uint64(buf[4]) << 32) | (uint64(buf[5]) << 40) |
+		(uint64(buf[6]) << 48) | (uint64(buf[7]) << 56))
+	bufOff := 8 // offset into buf
+
+	*p = Poly{} // zero p
+	for i := uint16(256 - 60); i < 256; i++ {
+		var b uint16
+
+		// Find location of where to move the new coefficient to using
+		// rejection sampling.
+		for {
+			if bufOff >= 136 {
+				h.Read(buf[:])
+				bufOff = 0
+			}
+
+			b = uint16(buf[bufOff])
+			bufOff++
+
+			if b <= i {
+				break
+			}
+		}
+
+		p[i] = p[b]
+		p[b] = 1
+		// Takes least significant bit of signs and uses it for the sign.
+		// Note 1 ^ (1 | (Q-1)) = Q-1.
+		p[b] ^= uint32((-(signs & 1)) & (1 | (Q - 1)))
+		signs >>= 1
+	}
+}
diff --git a/sign/dilithium/mode3/internal/sample_test.go b/sign/dilithium/mode3/internal/sample_test.go
new file mode 100644
index 000000000..5a59152db
--- /dev/null
+++ b/sign/dilithium/mode3/internal/sample_test.go
@@ -0,0 +1,173 @@
+package internal
+
+import (
+	"encoding/binary"
+	"testing"
+)
+
+func TestVectorDeriveUniform(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	p2 := Poly{
+		2901364, 562527, 5258502, 3885002, 4190126, 4460268, 6884052,
+		3514511, 5383040, 213206, 2155865, 5179607, 3551954, 2312357,
+		6066350, 8126097, 1179080, 4787182, 6552182, 6713644,
+		1561067, 7626063, 7859743, 5052321, 7032876, 7815031, 157938,
+		1865184, 490802, 5717642, 3451902, 7000218, 3743250, 1677431,
+		1875427, 5596150, 671623, 3819041, 6247594, 1014875, 4933545,
+		7122446, 6682963, 3388398, 3335295, 943002, 1145083, 3113071,
+		105967, 1916675, 7474561, 1107006, 700548, 2147909, 1603855,
+		5049181, 437882, 6118899, 5656914, 6731065, 3066622, 865453,
+		5427634, 981549, 4650873, 861291, 4003872, 5104220, 6171453,
+		3723302, 7426315, 6137283, 4874820, 6052561, 53441, 5032874,
+		5614778, 2248550, 1756499, 8280764, 8263880, 7600081,
+		5118374, 795344, 7543392, 6869925, 1841187, 4181568, 584562,
+		7483939, 4938664, 6863397, 5126354, 5218129, 6236086,
+		4149293, 379169, 4368487, 7490569, 3409215, 1580463, 3081737,
+		1278732, 7109719, 7371700, 2097931, 399836, 1700274, 7188595,
+		6830029, 1548850, 6593138, 6849097, 1518037, 2859442,
+		7772265, 7325153, 3281191, 7856131, 4995056, 4684325,
+		1351194, 8223904, 6817307, 2484146, 131782, 397032, 7436778,
+		7973479, 3171829, 5624626, 3540123, 7150120, 8313283,
+		3604714, 1043574, 117692, 7797783, 7909392, 903315, 7335342,
+		7501562, 5826142, 2709813, 8245473, 2369045, 2782257,
+		5762833, 6474114, 6862031, 424522, 594248, 2626630, 7659983,
+		5642869, 4075194, 1592129, 245547, 5271031, 3205046, 982375,
+		267873, 1286496, 7230481, 3208972, 7485411, 676111, 4944500,
+		2959742, 5934456, 1414847, 6067948, 1709895, 4648315, 126008,
+		8258986, 2183134, 2302072, 4674924, 4306056, 7465311,
+		6500270, 4247428, 4016815, 4973426, 294287, 2456847, 3289700,
+		2732169, 1159447, 5569724, 140001, 3237977, 8007761, 5874533,
+		255652, 3119586, 2102434, 6248250, 8152822, 8006066, 7708625,
+		6997719, 6260212, 6186962, 6636650, 7836834, 7998017,
+		2061516, 1197591, 1706544, 733027, 2392907, 2700000, 8254598,
+		4488002, 160495, 2985325, 2036837, 2703633, 6406550, 3579947,
+		6195178, 5552390, 6804584, 6305468, 5731980, 6095195,
+		3323409, 1322661, 6690942, 3374630, 5615167, 479044, 3136054,
+		4380418, 2833144, 7829577, 1770522, 6056687, 240415, 14780,
+		3740517, 5224226, 3547288, 2083124, 4699398, 3654239,
+		5624978, 585593, 3655369, 2281739, 3338565, 1908093, 7784706,
+		4352830,
+	}
+	for i := 0; i < 32; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniform(&seed, 30000)
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
+
+func TestDeriveUniform(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	for i := 0; i < 100; i++ {
+		binary.LittleEndian.PutUint64(seed[:], uint64(i))
+		p.DeriveUniform(&seed, uint16(i))
+		if !p.Normalized() {
+			t.Fatal()
+		}
+	}
+}
+
+func TestDeriveUniformLeqEta(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	for i := 0; i < 100; i++ {
+		binary.LittleEndian.PutUint64(seed[:], uint64(i))
+		p.DeriveUniformLeqEta(&seed, uint16(i))
+		for j := 0; j < N; j++ {
+			if p[j] < Q-Eta || p[j] > Q+Eta {
+				t.Fatal()
+			}
+		}
+	}
+}
+
+func TestDeriveUniformLeGamma1(t *testing.T) {
+	var p Poly
+	var seed [48]byte
+	for i := 0; i < 100; i++ {
+		binary.LittleEndian.PutUint64(seed[:], uint64(i))
+		p.DeriveUniformLeGamma1(&seed, uint16(i))
+		for j := 0; j < N; j++ {
+			if p[j] < Q-Gamma1 || p[j] > Q+Gamma1 {
+				t.Fatal()
+			}
+		}
+	}
+}
+
+func TestVectorDeriveUniformLeGamma1(t *testing.T) {
+	var p Poly
+	var seed [48]byte
+	p2 := Poly{
+		338083, 7978692, 8044913, 373628, 427855, 79725, 91018,
+		349821, 501552, 7955127, 8316400, 290708, 216142, 8199666,
+		8040144, 109426, 8177916, 8200218, 8125680, 358131, 160961,
+		497383, 25361, 156297, 8033745, 7897189, 48397, 498732,
+		464556, 7862704, 8308667, 236080, 91240, 8328377, 326190,
+		509979, 8313264, 8106493, 8210965, 8328036, 172602, 8108765,
+		8192963, 8361660, 8026473, 7932022, 322006, 8305874, 8254440,
+		7866474, 373371, 8023413, 8221878, 486124, 36080, 8324512,
+		8042056, 7984472, 8048111, 7910387, 8205382, 8259636, 188609,
+		8316587, 59306, 119803, 8067108, 8155455, 8153450, 7983908,
+		8222256, 21521, 114297, 8069037, 151621, 8014482, 8052856,
+		376107, 8004652, 175001, 8079461, 8351123, 8021484, 144547,
+		7908116, 8278100, 8136941, 142399, 8026843, 8081852, 124334,
+		242796, 266768, 7919478, 7954016, 28927, 8329064, 514031,
+		423911, 27907, 8142788, 8078298, 273978, 382723, 8148646,
+		186476, 8030712, 8067268, 356250, 145817, 60045, 122764,
+		56856, 8225416, 136437, 199652, 8343127, 8109765, 7936848,
+		446966, 8351681, 288663, 409663, 512988, 8350788, 8191864,
+		8366223, 281267, 7921696, 8213978, 442484, 67457, 8030602,
+		238514, 230458, 8301866, 8359700, 150320, 143893, 461735,
+		225443, 8027502, 151113, 365244, 7911438, 82498, 405398,
+		8207009, 8108255, 367485, 514660, 8294055, 8168958, 127725,
+		402955, 8051625, 7859029, 7980052, 321819, 7949587, 125778,
+		8287078, 131972, 499609, 256795, 8180323, 8269393, 5878,
+		8145473, 8238676, 383855, 415547, 424071, 241989, 8165743,
+		8207329, 149608, 8315331, 7901850, 8114275, 360650, 516061,
+		255090, 8277977, 270877, 8125200, 479248, 7991711, 8028595,
+		73426, 8215429, 208217, 153872, 429336, 229856, 461204,
+		236682, 7930158, 8298847, 228327, 8009399, 8111520, 345025,
+		386495, 93454, 8336429, 8161305, 7980811, 286795, 162808,
+		224476, 7972825, 85118, 287488, 8029791, 119071, 371, 518524,
+		473496, 451205, 127000, 19233, 211519, 2442, 7950959, 481888,
+		8041598, 8281176, 437202, 7912610, 8080153, 8237500, 7926828,
+		8009421, 204880, 62495, 8192310, 8314388, 98616, 182368,
+		323894, 59901, 481049, 8139275, 7872144, 254106, 376257,
+		93339, 301342, 366536, 438920, 84773, 461471, 8125755,
+		7930085, 405116,
+	}
+	for i := 0; i < 48; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniformLeGamma1(&seed, 30000)
+	p.Normalize()
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
+
+func TestDeriveUniformB60(t *testing.T) {
+	var p Poly
+	var w1 VecK
+	var seed [48]byte
+	for i := 0; i < 100; i++ {
+		binary.LittleEndian.PutUint64(seed[:], uint64(i))
+		p.DeriveUniformB60(&seed, &w1)
+		nonzero := 0
+		for j := 0; j < N; j++ {
+			if p[j] != 0 {
+				if p[j] != 1 && p[j] != Q-1 {
+					t.Fatal()
+				}
+				nonzero++
+			}
+		}
+		if nonzero != 60 {
+			t.Fatal()
+		}
+	}
+}
diff --git a/sign/dilithium/mode3/internal/vec.go b/sign/dilithium/mode3/internal/vec.go
new file mode 100644
index 000000000..7a71a5a2c
--- /dev/null
+++ b/sign/dilithium/mode3/internal/vec.go
@@ -0,0 +1,273 @@
+package internal
+
+// A vector of L polynomials
+type VecL [L]Poly
+
+// A vector of L polynomials
+type VecK [K]Poly
+
+// Normalize the polynomials in this vector
+func (v *VecL) Normalize() {
+	for i := 0; i < L; i++ {
+		v[i].Normalize()
+	}
+}
+
+// Normalize the polynomials in this vector assuming their coefficients
+// are already bounded by 2q.
+func (v *VecL) NormalizeAssumingLe2Q() {
+	for i := 0; i < L; i++ {
+		v[i].NormalizeAssumingLe2Q()
+	}
+}
+
+// Sets v to w + u.  Does not normalize.
+func (v *VecL) Add(w, u *VecL) {
+	for i := 0; i < L; i++ {
+		v[i].Add(&w[i], &u[i])
+	}
+}
+
+// Applies NTT componentwise. See Poly.NTT() for details.
+func (v *VecL) NTT() {
+	for i := 0; i < L; i++ {
+		v[i].NTT()
+	}
+}
+
+// Checks whether any of the coefficients exceeds the given bound in supnorm
+//
+// Requires the vector to be normalized.
+func (v *VecL) Exceeds(bound uint32) bool {
+	for i := 0; i < L; i++ {
+		if v[i].Exceeds(bound) {
+			return true
+		}
+	}
+	return false
+}
+
+// Applies Poly.Power2Round componentwise.
+//
+// Requires the vector to be normalized.
+func (v *VecL) Power2Round(v0PlusQ, v1 *VecL) {
+	for i := 0; i < L; i++ {
+		v[i].Power2Round(&v0PlusQ[i], &v1[i])
+	}
+}
+
+// Applies Poly.Decompose componentwise.
+//
+// Requires the vector to be normalized.
+func (v *VecL) Decmopose(v0PlusQ, v1 *VecL) {
+	for i := 0; i < L; i++ {
+		v[i].Decompose(&v0PlusQ[i], &v1[i])
+	}
+}
+
+// Sequentially packs each polynomial using Poly.PackLeqEta().
+func (v *VecL) PackLeqEta(buf []byte) {
+	offset := 0
+	for i := 0; i < L; i++ {
+		v[i].PackLeqEta(buf[offset:])
+		offset += PolyLeqEtaSize
+	}
+}
+
+// Sets v to the polynomials packed in buf using VecL.PackLeqEta().
+func (v *VecL) UnpackLeqEta(buf []byte) {
+	offset := 0
+	for i := 0; i < L; i++ {
+		v[i].UnpackLeqEta(buf[offset:])
+		offset += PolyLeqEtaSize
+	}
+}
+
+// Sequentially packs each polynomial using Poly.PackLeGamma1().
+func (v *VecL) PackLeGamma1(buf []byte) {
+	offset := 0
+	for i := 0; i < L; i++ {
+		v[i].PackLeGamma1(buf[offset:])
+		offset += PolyLeGamma1Size
+	}
+}
+
+// Sets v to the polynomials packed in buf using VecL.PackLeGamma1().
+func (v *VecL) UnpackLeGamma1(buf []byte) {
+	offset := 0
+	for i := 0; i < L; i++ {
+		v[i].UnpackLeGamma1(buf[offset:])
+		offset += PolyLeGamma1Size
+	}
+}
+
+// Normalize the polynomials in this vector
+func (v *VecK) Normalize() {
+	for i := 0; i < K; i++ {
+		v[i].Normalize()
+	}
+}
+
+// Normalize the polynomials in this vector assuming their coefficients
+// are already bounded by 2q.
+func (v *VecK) NormalizeAssumingLe2Q() {
+	for i := 0; i < K; i++ {
+		v[i].NormalizeAssumingLe2Q()
+	}
+}
+
+// Sets v to w + u.  Does not normalize.
+func (v *VecK) Add(w, u *VecK) {
+	for i := 0; i < K; i++ {
+		v[i].Add(&w[i], &u[i])
+	}
+}
+
+// Checks whether any of the coefficients exceeds the given bound in supnorm
+//
+// Requires the vector to be normalized.
+func (v *VecK) Exceeds(bound uint32) bool {
+	for i := 0; i < K; i++ {
+		if v[i].Exceeds(bound) {
+			return true
+		}
+	}
+	return false
+}
+
+// Applies Poly.Power2Round componentwise.
+//
+// Requires the vector to be normalized.
+func (v *VecK) Power2Round(v0PlusQ, v1 *VecK) {
+	for i := 0; i < K; i++ {
+		v[i].Power2Round(&v0PlusQ[i], &v1[i])
+	}
+}
+
+// Applies Poly.Decompose componentwise.
+//
+// Requires the vector to be normalized.
+func (v *VecK) Decompose(v0PlusQ, v1 *VecK) {
+	for i := 0; i < K; i++ {
+		v[i].Decompose(&v0PlusQ[i], &v1[i])
+	}
+}
+
+// Sets p to the hint polynomail for low part p0 and high part p1.
+//
+// Returns the number of ones in the hint vector.
+func (v *VecK) MakeHint(v0, v1 *VecK) (pop uint32) {
+	for i := 0; i < K; i++ {
+		pop += v[i].MakeHint(&v0[i], &v1[i])
+	}
+	return
+}
+
+// Computes corrections to the high bits of the polynomials in the vector
+// w using the hints in h and sets v to the corrected high bits.  Returns v.
+func (v *VecK) UseHint(q, hint *VecK) *VecK {
+	for i := 0; i < K; i++ {
+		v[i].UseHint(&q[i], &hint[i])
+	}
+	return v
+}
+
+// Sequentially packs each polynomial using Poly.PackT1().
+func (v *VecK) PackT1(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].PackT1(buf[offset:])
+		offset += PolyT1Size
+	}
+}
+
+// Sets v to the vector packed into buf by PackT1().
+func (v *VecK) UnpackT1(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].UnpackT1(buf[offset:])
+		offset += PolyT1Size
+	}
+}
+
+// Sequentially packs each polynomial using Poly.PackT0().
+func (v *VecK) PackT0(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].PackT0(buf[offset:])
+		offset += PolyT0Size
+	}
+}
+
+// Sets v to the vector packed into buf by PackT0().
+func (v *VecK) UnpackT0(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].UnpackT0(buf[offset:])
+		offset += PolyT0Size
+	}
+}
+
+// Sequentially packs each polynomial using Poly.PackLeqEta().
+func (v *VecK) PackLeqEta(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].PackLeqEta(buf[offset:])
+		offset += PolyLeqEtaSize
+	}
+}
+
+// Sets v to the polynomials packed in buf using VecK.PackLeqEta().
+func (v *VecK) UnpackLeqEta(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].UnpackLeqEta(buf[offset:])
+		offset += PolyLeqEtaSize
+	}
+}
+
+// Applies NTT componentwise. See Poly.NTT() for details.
+func (v *VecK) NTT() {
+	for i := 0; i < K; i++ {
+		v[i].NTT()
+	}
+}
+
+// Sequentially packs each polynomial using Poly.PackLe16().
+func (v *VecK) PackLe16(buf []byte) {
+	offset := 0
+	for i := 0; i < K; i++ {
+		v[i].PackLe16(buf[offset:])
+		offset += PolyLe16Size
+	}
+}
+
+// Sets v to a - b.
+//
+// Warning: assumes coefficients of the polynomials of  b are less than 2q.
+func (v *VecK) Sub(a, b *VecK) {
+	for i := 0; i < K; i++ {
+		v[i].Sub(&a[i], &b[i])
+	}
+}
+
+// Sets v to 2^d w without reducing.
+func (v *VecK) MulBy2toD(w *VecK) {
+	for i := 0; i < K; i++ {
+		v[i].MulBy2toD(&w[i])
+	}
+}
+
+// Applies InvNTT componentwise. See Poly.InvNTT() for details.
+func (v *VecK) InvNTT() {
+	for i := 0; i < K; i++ {
+		v[i].InvNTT()
+	}
+}
+
+// Applies Poly.ReduceLe2Q() componentwise.
+func (v *VecK) ReduceLe2Q() {
+	for i := 0; i < K; i++ {
+		v[i].ReduceLe2Q()
+	}
+}
diff --git a/sign/dilithium/mode4/internal/dilithium.go b/sign/dilithium/mode4/internal/dilithium.go
new file mode 120000
index 000000000..f4eea29aa
--- /dev/null
+++ b/sign/dilithium/mode4/internal/dilithium.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/dilithium_test.go b/sign/dilithium/mode4/internal/dilithium_test.go
new file mode 120000
index 000000000..baa8508b5
--- /dev/null
+++ b/sign/dilithium/mode4/internal/dilithium_test.go
@@ -0,0 +1 @@
+../../mode3/internal/dilithium_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/field.go b/sign/dilithium/mode4/internal/field.go
new file mode 120000
index 000000000..2a2a32e6d
--- /dev/null
+++ b/sign/dilithium/mode4/internal/field.go
@@ -0,0 +1 @@
+../../mode3/internal/field.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/field_test.go b/sign/dilithium/mode4/internal/field_test.go
new file mode 120000
index 000000000..3e9043a03
--- /dev/null
+++ b/sign/dilithium/mode4/internal/field_test.go
@@ -0,0 +1 @@
+../../mode3/internal/field_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/kat_test.go b/sign/dilithium/mode4/internal/kat_test.go
new file mode 120000
index 000000000..54cd221f1
--- /dev/null
+++ b/sign/dilithium/mode4/internal/kat_test.go
@@ -0,0 +1 @@
+../../mode3/internal/kat_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/mat.go b/sign/dilithium/mode4/internal/mat.go
new file mode 120000
index 000000000..0601aa926
--- /dev/null
+++ b/sign/dilithium/mode4/internal/mat.go
@@ -0,0 +1 @@
+../../mode3/internal/mat.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/ntt.go b/sign/dilithium/mode4/internal/ntt.go
new file mode 120000
index 000000000..9b9aa00b8
--- /dev/null
+++ b/sign/dilithium/mode4/internal/ntt.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/ntt_test.go b/sign/dilithium/mode4/internal/ntt_test.go
new file mode 120000
index 000000000..02efd77bc
--- /dev/null
+++ b/sign/dilithium/mode4/internal/ntt_test.go
@@ -0,0 +1 @@
+../../mode3/internal/ntt_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/pack.go b/sign/dilithium/mode4/internal/pack.go
new file mode 120000
index 000000000..578d3dc68
--- /dev/null
+++ b/sign/dilithium/mode4/internal/pack.go
@@ -0,0 +1 @@
+../../mode3/internal/pack.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/pack_test.go b/sign/dilithium/mode4/internal/pack_test.go
new file mode 120000
index 000000000..5571475bb
--- /dev/null
+++ b/sign/dilithium/mode4/internal/pack_test.go
@@ -0,0 +1 @@
+../../mode3/internal/pack_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/params.go b/sign/dilithium/mode4/internal/params.go
new file mode 100644
index 000000000..b9b141064
--- /dev/null
+++ b/sign/dilithium/mode4/internal/params.go
@@ -0,0 +1,14 @@
+package internal
+
+const (
+	Name           = "Dilithium4"
+	PublicKeySize  = 1760
+	PrivateKeySize = 3856
+	SignatureSize  = 3366
+	K              = 6
+	L              = 5
+	Eta            = 3
+	DoubleEtaBits  = 3
+	Beta           = 175
+	Omega          = 120
+)
diff --git a/sign/dilithium/mode4/internal/params_test.go b/sign/dilithium/mode4/internal/params_test.go
new file mode 100644
index 000000000..efcb0eb0b
--- /dev/null
+++ b/sign/dilithium/mode4/internal/params_test.go
@@ -0,0 +1,83 @@
+package internal
+
+import (
+	"testing"
+
+	"github.com/cloudflare/circl/internal/shake"
+)
+
+// Tests specific to the current mode
+
+func hash(in []byte) [16]byte {
+	var ret [16]byte
+	h := shake.NewShake256()
+	h.Write(in)
+	h.Read(ret[:])
+	return ret
+}
+
+func TestNewKeyFromSeed(t *testing.T) {
+	var seed [SeedSize]byte
+	var pkp [PublicKeySize]byte
+	var skp [PrivateKeySize]byte
+	pk, sk := NewKeyFromSeed(&seed)
+	pk.Pack(&pkp)
+	sk.Pack(&skp)
+
+	// Generated with reference implementation.
+	ehpk := [16]byte{
+		231, 153, 127, 199, 26, 103, 150, 5,
+		109, 70, 51, 164, 7, 105, 196, 149,
+	}
+	ehsk := [16]byte{
+		224, 84, 49, 155, 186, 189, 46, 21,
+		108, 86, 232, 238, 146, 60, 42, 142,
+	}
+
+	if hash(pkp[:]) != ehpk {
+		t.Fatalf("pk not ok")
+	}
+	if hash(skp[:]) != ehsk {
+		t.Fatalf("sk not ok")
+	}
+}
+
+func TestVectorDeriveUniformLeqEta(t *testing.T) {
+	var p Poly
+	var seed [32]byte
+	p2 := Poly{
+		8380416, 0, 8380415, 1, 8380416, 8380415, 1, 8380415,
+		8380416, 3, 8380416, 8380414, 3, 8380415, 2, 3, 1, 1,
+		8380414, 1, 0, 2, 8380414, 8380415, 8380415, 8380416,
+		8380414, 3, 8380416, 8380414, 8380416, 8380414, 1, 3, 2,
+		0, 8380416, 3, 8380415, 8380415, 2, 3, 8380416, 0, 2, 1,
+		1, 8380414, 8380414, 3, 1, 0, 8380416, 2, 1, 8380415, 3,
+		2, 8380414, 3, 3, 8380414, 1, 8380414, 1, 2, 8380414, 2,
+		8380415, 8380415, 3, 1, 0, 8380415, 8380414, 8380415, 3,
+		8380416, 8380416, 8380414, 2, 8380416, 8380414, 0, 0, 3,
+		8380414, 8380415, 2, 8380416, 8380414, 0, 3, 2, 1, 8380414,
+		0, 1, 0, 8380416, 2, 2, 0, 8380414, 0, 2, 8380414, 8380416,
+		1, 8380415, 1, 8380416, 3, 8380414, 0, 2, 8380415, 3,
+		8380415, 8380414, 8380416, 2, 3, 1, 0, 0, 2, 1, 3, 0, 1,
+		3, 3, 8380415, 0, 8380416, 3, 1, 2, 8380415, 0, 3, 1,
+		8380416, 8380414, 2, 8380414, 8380414, 0, 8380415, 8380416,
+		8380415, 1, 8380416, 3, 8380414, 1, 1, 3, 8380414, 3, 2,
+		2, 0, 3, 8380416, 2, 8380416, 0, 8380415, 8380416, 8380416,
+		3, 8380414, 0, 0, 8380415, 8380414, 8380416, 2, 1, 8380416,
+		2, 2, 8380414, 1, 1, 2, 1, 8380415, 8380416, 1, 1, 0, 1,
+		8380416, 0, 8380416, 2, 8380415, 1, 0, 1, 8380414, 3,
+		8380414, 8380415, 8380416, 0, 0, 1, 2, 8380415, 8380414,
+		0, 8380415, 0, 2, 0, 8380416, 1, 8380415, 1, 8380416, 0,
+		0, 8380414, 3, 0, 3, 0, 2, 3, 0, 3, 8380416, 8380415,
+		8380415, 8380416, 2, 3, 0, 3, 2, 8380415, 3, 2, 8380416,
+		8380415, 3, 1, 0, 8380414, 2, 8380416, 3,
+	}
+	for i := 0; i < 32; i++ {
+		seed[i] = byte(i)
+	}
+	p.DeriveUniformLeqEta(&seed, 30000)
+	p.Normalize()
+	if p != p2 {
+		t.Fatalf("%v != %v", p, p2)
+	}
+}
diff --git a/sign/dilithium/mode4/internal/poly.go b/sign/dilithium/mode4/internal/poly.go
new file mode 120000
index 000000000..67843bdd7
--- /dev/null
+++ b/sign/dilithium/mode4/internal/poly.go
@@ -0,0 +1 @@
+../../mode3/internal/poly.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/sample.go b/sign/dilithium/mode4/internal/sample.go
new file mode 120000
index 000000000..7db77123c
--- /dev/null
+++ b/sign/dilithium/mode4/internal/sample.go
@@ -0,0 +1 @@
+../../mode3/internal/sample.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/sample_test.go b/sign/dilithium/mode4/internal/sample_test.go
new file mode 120000
index 000000000..a7d78fe33
--- /dev/null
+++ b/sign/dilithium/mode4/internal/sample_test.go
@@ -0,0 +1 @@
+../../mode3/internal/sample_test.go
\ No newline at end of file
diff --git a/sign/dilithium/mode4/internal/vec.go b/sign/dilithium/mode4/internal/vec.go
new file mode 120000
index 000000000..82e7f9b35
--- /dev/null
+++ b/sign/dilithium/mode4/internal/vec.go
@@ -0,0 +1 @@
+../../mode3/internal/vec.go
\ No newline at end of file