ntt.go raw
1 package crypto
2
3 // Number-Theoretic Transform (NTT) over Z_257 for Hamadryad.
4 //
5 // The NTT is the finite-field analogue of the FFT. It transforms
6 // a polynomial from coefficient form to evaluation form, enabling
7 // O(n log n) polynomial multiplication via pointwise multiply.
8 //
9 // For Hamadryad: n=64, p=257.
10 //
11 // 257 is prime and 256 = 257-1 = 2^8, so Z_257* has order 256.
12 // g = 3 is a primitive root (generator) of Z_257*.
13 // psi = 3^(256/128) = 3^2 = 9 is a primitive 128th root of unity.
14 //
15 // For polynomials mod (x^64 + 1), we use the negacyclic NTT:
16 // pre-multiply coefficients by powers of psi (a 128th root),
17 // then apply a standard 64-point NTT with omega = psi^2 = 81
18 // (a primitive 64th root of unity).
19
20 // mod257 reduces x into [0, 256].
21 func mod257(x int) uint16 {
22 x %= HamP
23 if x < 0 {
24 x += HamP
25 }
26 return uint16(x)
27 }
28
29 // powMod computes base^exp mod 257.
30 func powMod(base, exp int) uint16 {
31 result := 1
32 b := base % HamP
33 if b < 0 {
34 b += HamP
35 }
36 for e := exp; e > 0; e >>= 1 {
37 if e&1 == 1 {
38 result = (result * b) % HamP
39 }
40 b = (b * b) % HamP
41 }
42 return uint16(result)
43 }
44
45 // invMod computes the modular inverse of a mod 257.
46 // Uses Fermat's little theorem: a^(-1) = a^(p-2) mod p.
47 func invMod(a uint16) uint16 {
48 return powMod(int(a), HamP-2)
49 }
50
51 // bitRev6 reverses the low 6 bits of x.
52 func bitRev6(x int) int {
53 r := 0
54 for range 6 {
55 r = (r << 1) | (x & 1)
56 x >>= 1
57 }
58 return r
59 }
60
61 // Pre-computed tables, filled by initNTTTables.
62 var (
63 // psiPows[i] = psi^i mod 257, for i=0..127.
64 // psi = 9 (primitive 128th root of unity).
65 psiPows [128]uint16
66
67 // psiInvPows[i] = psi^(-i) mod 257.
68 psiInvPows [128]uint16
69
70 // invN = 64^(-1) mod 257.
71 invN uint16
72
73 nttTablesReady bool
74 )
75
76 // initNTTTables populates the NTT lookup tables. Safe to call multiple
77 // times; only the first call has effect. Called explicitly from the
78 // hamadryad init to guarantee table availability before key generation
79 // (Go runs init functions in source file alphabetical order within a
80 // package, and "hamadryad.go" sorts before "ntt.go").
81 func initNTTTables() {
82 if nttTablesReady {
83 return
84 }
85 const psi = 9 // primitive 128th root of unity in Z_257
86
87 psiPows[0] = 1
88 for i := 1; i < 128; i++ {
89 psiPows[i] = mod257(int(psiPows[i-1]) * psi)
90 }
91
92 psiInv := invMod(psi)
93 psiInvPows[0] = 1
94 for i := 1; i < 128; i++ {
95 psiInvPows[i] = mod257(int(psiInvPows[i-1]) * int(psiInv))
96 }
97
98 invN = invMod(uint16(HamN))
99 nttTablesReady = true
100 }
101
102 func init() {
103 initNTTTables()
104 }
105
106 // ntt64 computes the forward negacyclic NTT of a length-64 polynomial over Z_257.
107 //
108 // The negacyclic NTT computes polynomial evaluation at the roots of (x^64 + 1),
109 // which are psi^(2k+1) for k=0..63. Implementation:
110 // 1. Pre-multiply: a[i] *= psi^i
111 // 2. Standard radix-2 DIT NTT with omega = psi^2 (primitive 64th root)
112 func ntt64(a *[HamN]uint16) {
113 const n = HamN
114
115 // Step 1: pre-multiply by psi^i.
116 for i := range n {
117 a[i] = mod257(int(a[i]) * int(psiPows[i]))
118 }
119
120 // Step 2: bit-reversal permutation.
121 for i := range n {
122 j := bitRev6(i)
123 if i < j {
124 a[i], a[j] = a[j], a[i]
125 }
126 }
127
128 // Step 3: Cooley-Tukey butterfly stages.
129 // omega = psi^2 = 9^2 = 81 (primitive 64th root of unity).
130 for length := 1; length < n; length <<= 1 {
131 // Twiddle factor base: omega^(n/(2*length)).
132 step := n / (2 * length)
133 for start := 0; start < n; start += 2 * length {
134 for j := range length {
135 // tw = omega^(j*step) = psi^(2*j*step)
136 tw := psiPows[(2*j*step)%128]
137 u := a[start+j]
138 v := mod257(int(a[start+j+length]) * int(tw))
139 a[start+j] = mod257(int(u) + int(v))
140 a[start+j+length] = mod257(int(u) - int(v))
141 }
142 }
143 }
144 }
145
146 // intt64 computes the inverse negacyclic NTT, recovering coefficients.
147 func intt64(a *[HamN]uint16) {
148 const n = HamN
149
150 // Step 1: bit-reversal permutation.
151 for i := range n {
152 j := bitRev6(i)
153 if i < j {
154 a[i], a[j] = a[j], a[i]
155 }
156 }
157
158 // Step 2: Gentleman-Sande butterfly (inverse DIT) with omega^(-1).
159 for length := 1; length < n; length <<= 1 {
160 step := n / (2 * length)
161 for start := 0; start < n; start += 2 * length {
162 for j := range length {
163 // tw = omega^(-(j*step)) = psi^(-2*j*step)
164 tw := psiInvPows[(2*j*step)%128]
165 u := a[start+j]
166 v := mod257(int(a[start+j+length]) * int(tw))
167 a[start+j] = mod257(int(u) + int(v))
168 a[start+j+length] = mod257(int(u) - int(v))
169 }
170 }
171 }
172
173 // Step 3: multiply by 1/n.
174 for i := range n {
175 a[i] = mod257(int(a[i]) * int(invN))
176 }
177
178 // Step 4: undo psi pre-multiplication: a[i] *= psi^(-i).
179 for i := range n {
180 a[i] = mod257(int(a[i]) * int(psiInvPows[i]))
181 }
182 }
183