ntt27.mx raw
1 package hamcrypto
2
3 // Radix-3 NTT over Z_271 for trinary Hamadryad.
4 // n=27 = 3^3, p=271. psi=188 (primitive 54th root of unity).
5
6 const GnarlP = 271
7 const GnarlN = 27
8
9 // mod271 reduces x into [0, 270] using Barrett reduction.
10 func mod271(x uint32) (r uint16) {
11 q := (x * 483) >> 17
12 rv := x - q*271
13 if rv >= 271 {
14 rv -= 271
15 }
16 return uint16(rv)
17 }
18
19 // mod271s reduces a signed int into [0, 270].
20 func mod271s(x int32) (r uint16) {
21 xv := x % GnarlP
22 if xv < 0 {
23 xv += GnarlP
24 }
25 return uint16(xv)
26 }
27
28 // powMod271 computes base^exp mod 271.
29 func powMod271(base, exp int32) (r uint16) {
30 result := int32(1)
31 b := base % GnarlP
32 if b < 0 {
33 b += GnarlP
34 }
35 e := exp
36 for e > 0 {
37 if e&1 == 1 {
38 result = (result * b) % GnarlP
39 }
40 b = (b * b) % GnarlP
41 e >>= 1
42 }
43 return uint16(result)
44 }
45
46 // invMod271 computes the modular inverse of a mod 271.
47 func invMod271(a uint16) (r uint16) {
48 return powMod271(int32(a), GnarlP-2)
49 }
50
51 // digitRev3 reverses the base-3 digits of x with the given number of digits.
52 func digitRev3(x, digits int32) (r int32) {
53 r = 0
54 xv := x
55 for i := int32(0); i < digits; i++ {
56 r = r*3 + xv%3
57 xv /= 3
58 }
59 return r
60 }
61
62 // Pre-computed tables for the radix-3 NTT.
63 var psiPows27 [54]uint16
64 var psiInvPows27 [54]uint16
65 var invN27 uint16
66 var ntt27TablesReady bool
67
68 // initNTT27Tables populates the radix-3 NTT lookup tables.
69 func initNTT27Tables() {
70 if ntt27TablesReady {
71 return
72 }
73 const psi = 188
74
75 psiPows27[0] = 1
76 for i := int32(1); i < 54; i++ {
77 psiPows27[i] = mod271(uint32(psiPows27[i-1]) * psi)
78 }
79
80 psiInv := invMod271(psi)
81 psiInvPows27[0] = 1
82 for i := int32(1); i < 54; i++ {
83 psiInvPows27[i] = mod271(uint32(psiInvPows27[i-1]) * uint32(psiInv))
84 }
85
86 invN27 = invMod271(uint16(GnarlN))
87 ntt27TablesReady = true
88 }
89
90 // ntt27 computes the forward negacyclic NTT of a length-27 polynomial over Z_271.
91 func ntt27(a *[GnarlN]uint16) {
92 // Pre-multiply by psi^i for negacyclic.
93 a[1] = mod271(uint32(a[1]) * 188)
94 a[2] = mod271(uint32(a[2]) * 114)
95 a[3] = mod271(uint32(a[3]) * 23)
96 a[4] = mod271(uint32(a[4]) * 259)
97 a[5] = mod271(uint32(a[5]) * 183)
98 a[6] = mod271(uint32(a[6]) * 258)
99 a[7] = mod271(uint32(a[7]) * 266)
100 a[8] = mod271(uint32(a[8]) * 144)
101 a[9] = mod271(uint32(a[9]) * 243)
102 a[10] = mod271(uint32(a[10]) * 156)
103 a[11] = mod271(uint32(a[11]) * 60)
104 a[12] = mod271(uint32(a[12]) * 169)
105 a[13] = mod271(uint32(a[13]) * 65)
106 a[14] = mod271(uint32(a[14]) * 25)
107 a[15] = mod271(uint32(a[15]) * 93)
108 a[16] = mod271(uint32(a[16]) * 140)
109 a[17] = mod271(uint32(a[17]) * 33)
110 a[18] = mod271(uint32(a[18]) * 242)
111 a[19] = mod271(uint32(a[19]) * 239)
112 a[20] = mod271(uint32(a[20]) * 217)
113 a[21] = mod271(uint32(a[21]) * 146)
114 a[22] = mod271(uint32(a[22]) * 77)
115 a[23] = mod271(uint32(a[23]) * 113)
116 a[24] = mod271(uint32(a[24]) * 106)
117 a[25] = mod271(uint32(a[25]) * 145)
118 a[26] = mod271(uint32(a[26]) * 160)
119
120 // Digit-reversal permutation (9 swaps).
121 a[1], a[9] = a[9], a[1]
122 a[2], a[18] = a[18], a[2]
123 a[4], a[12] = a[12], a[4]
124 a[5], a[21] = a[21], a[5]
125 a[7], a[15] = a[15], a[7]
126 a[8], a[24] = a[24], a[8]
127 a[11], a[19] = a[19], a[11]
128 a[14], a[22] = a[22], a[14]
129 a[17], a[25] = a[25], a[17]
130
131 // Stage 0: stride=1, all tw1=tw2=1
132 nttBfly1(a, 0, 1, 2)
133 nttBfly1(a, 3, 4, 5)
134 nttBfly1(a, 6, 7, 8)
135 nttBfly1(a, 9, 10, 11)
136 nttBfly1(a, 12, 13, 14)
137 nttBfly1(a, 15, 16, 17)
138 nttBfly1(a, 18, 19, 20)
139 nttBfly1(a, 21, 22, 23)
140 nttBfly1(a, 24, 25, 26)
141
142 // Stage 1: stride=3
143 nttBfly1(a, 0, 3, 6)
144 nttBfly1(a, 9, 12, 15)
145 nttBfly1(a, 18, 21, 24)
146 nttBfly(a, 1, 4, 7, 258, 169)
147 nttBfly(a, 10, 13, 16, 258, 169)
148 nttBfly(a, 19, 22, 25, 258, 169)
149 nttBfly(a, 2, 5, 8, 169, 106)
150 nttBfly(a, 11, 14, 17, 169, 106)
151 nttBfly(a, 20, 23, 26, 169, 106)
152
153 // Stage 2: stride=9
154 nttBfly1(a, 0, 9, 18)
155 nttBfly(a, 1, 10, 19, 114, 259)
156 nttBfly(a, 2, 11, 20, 259, 144)
157 nttBfly(a, 3, 12, 21, 258, 169)
158 nttBfly(a, 4, 13, 22, 144, 140)
159 nttBfly(a, 5, 14, 23, 156, 217)
160 nttBfly(a, 6, 15, 24, 169, 106)
161 nttBfly(a, 7, 16, 25, 25, 83)
162 nttBfly(a, 8, 17, 26, 140, 88)
163 }
164
165 // nttBfly1 performs a radix-3 butterfly with trivial twiddles.
166 func nttBfly1(a *[27]uint16, i0, i1, i2 int32) {
167 v0 := uint32(a[i0])
168 v1 := uint32(a[i1])
169 v2 := uint32(a[i2])
170 a[i0] = mod271(v0 + v1 + v2)
171 a[i1] = mod271(v0 + uint32(mod271(v1*242)) + uint32(mod271(v2*28)))
172 a[i2] = mod271(v0 + uint32(mod271(v1*28)) + uint32(mod271(v2*242)))
173 }
174
175 // nttBfly performs a radix-3 butterfly with given twiddle factors.
176 func nttBfly(a *[27]uint16, i0, i1, i2 int32, tw1, tw2 uint32) {
177 v0 := uint32(a[i0])
178 a1tw := uint32(mod271(uint32(a[i1]) * tw1))
179 a2tw := uint32(mod271(uint32(a[i2]) * tw2))
180 a[i0] = mod271(v0 + a1tw + a2tw)
181 a[i1] = mod271(v0 + uint32(mod271(a1tw*242)) + uint32(mod271(a2tw*28)))
182 a[i2] = mod271(v0 + uint32(mod271(a1tw*28)) + uint32(mod271(a2tw*242)))
183 }
184
185 // intt27 computes the inverse negacyclic NTT, recovering coefficients.
186 func intt27(a *[GnarlN]uint16) {
187 // Stage 2 (DIF, top-down)
188 inttBfly1(a, 0, 9, 18)
189 inttBfly(a, 1, 10, 19, 126, 158)
190 inttBfly(a, 2, 11, 20, 158, 32)
191 inttBfly(a, 3, 12, 21, 125, 178)
192 inttBfly(a, 4, 13, 22, 32, 211)
193 inttBfly(a, 5, 14, 23, 238, 5)
194 inttBfly(a, 6, 15, 24, 178, 248)
195 inttBfly(a, 7, 16, 25, 206, 160)
196 inttBfly(a, 8, 17, 26, 211, 77)
197
198 // Stage 1
199 inttBfly1(a, 0, 3, 6)
200 inttBfly(a, 1, 4, 7, 125, 178)
201 inttBfly(a, 2, 5, 8, 178, 248)
202 inttBfly1(a, 9, 12, 15)
203 inttBfly(a, 10, 13, 16, 125, 178)
204 inttBfly(a, 11, 14, 17, 178, 248)
205 inttBfly1(a, 18, 21, 24)
206 inttBfly(a, 19, 22, 25, 125, 178)
207 inttBfly(a, 20, 23, 26, 178, 248)
208
209 // Stage 0
210 inttBfly1(a, 0, 1, 2)
211 inttBfly1(a, 3, 4, 5)
212 inttBfly1(a, 6, 7, 8)
213 inttBfly1(a, 9, 10, 11)
214 inttBfly1(a, 12, 13, 14)
215 inttBfly1(a, 15, 16, 17)
216 inttBfly1(a, 18, 19, 20)
217 inttBfly1(a, 21, 22, 23)
218 inttBfly1(a, 24, 25, 26)
219
220 // Digit-reversal permutation.
221 a[1], a[9] = a[9], a[1]
222 a[2], a[18] = a[18], a[2]
223 a[4], a[12] = a[12], a[4]
224 a[5], a[21] = a[21], a[5]
225 a[7], a[15] = a[15], a[7]
226 a[8], a[24] = a[24], a[8]
227 a[11], a[19] = a[19], a[11]
228 a[14], a[22] = a[22], a[14]
229 a[17], a[25] = a[25], a[17]
230
231 // Post-multiply: fused invN27 * psiInvPows27[i].
232 a[0] = mod271(uint32(a[0]) * 261)
233 a[1] = mod271(uint32(a[1]) * 245)
234 a[2] = mod271(uint32(a[2]) * 95)
235 a[3] = mod271(uint32(a[3]) * 247)
236 a[4] = mod271(uint32(a[4]) * 46)
237 a[5] = mod271(uint32(a[5]) * 228)
238 a[6] = mod271(uint32(a[6]) * 105)
239 a[7] = mod271(uint32(a[7]) * 2)
240 a[8] = mod271(uint32(a[8]) * 222)
241 a[9] = mod271(uint32(a[9]) * 252)
242 a[10] = mod271(uint32(a[10]) * 59)
243 a[11] = mod271(uint32(a[11]) * 45)
244 a[12] = mod271(uint32(a[12]) * 117)
245 a[13] = mod271(uint32(a[13]) * 250)
246 a[14] = mod271(uint32(a[14]) * 108)
247 a[15] = mod271(uint32(a[15]) * 64)
248 a[16] = mod271(uint32(a[16]) * 58)
249 a[17] = mod271(uint32(a[17]) * 205)
250 a[18] = mod271(uint32(a[18]) * 262)
251 a[19] = mod271(uint32(a[19]) * 85)
252 a[20] = mod271(uint32(a[20]) * 221)
253 a[21] = mod271(uint32(a[21]) * 141)
254 a[22] = mod271(uint32(a[22]) * 204)
255 a[23] = mod271(uint32(a[23]) * 151)
256 a[24] = mod271(uint32(a[24]) * 230)
257 a[25] = mod271(uint32(a[25]) * 56)
258 a[26] = mod271(uint32(a[26]) * 254)
259 }
260
261 // inttBfly1 performs an inverse radix-3 DIF butterfly with trivial twiddles.
262 func inttBfly1(a *[27]uint16, i0, i1, i2 int32) {
263 v0 := uint32(a[i0])
264 v1 := uint32(a[i1])
265 v2 := uint32(a[i2])
266 a[i0] = mod271(v0 + v1 + v2)
267 a[i1] = mod271(v0 + uint32(mod271(v1*28)) + uint32(mod271(v2*242)))
268 a[i2] = mod271(v0 + uint32(mod271(v1*242)) + uint32(mod271(v2*28)))
269 }
270
271 // inttBfly performs an inverse radix-3 DIF butterfly with given twiddle factors.
272 func inttBfly(a *[27]uint16, i0, i1, i2 int32, tw1inv, tw2inv uint32) {
273 v0 := uint32(a[i0])
274 v1 := uint32(a[i1])
275 v2 := uint32(a[i2])
276 a[i0] = mod271(v0 + v1 + v2)
277 b1 := mod271(v0 + uint32(mod271(v1*28)) + uint32(mod271(v2*242)))
278 b2 := mod271(v0 + uint32(mod271(v1*242)) + uint32(mod271(v2*28)))
279 a[i1] = mod271(uint32(b1) * tw1inv)
280 a[i2] = mod271(uint32(b2) * tw2inv)
281 }
282