miden-crypto/src/dsa/rpo_falcon512/falcon_c/rpo.c
src/dsa/rpo_falcon512/falcon_c/rpo.cRaw
1/*
2 * RPO implementation.
3 */
4
5#include <stdint.h>
6#include <string.h>
7#include <stdlib.h>
8
9/* ================================================================================================
10 * Modular Arithmetic
11 */
12
13#define P 0xFFFFFFFF00000001
14#define M 12289
15
16// From https://github.com/ncw/iprime/blob/master/mod_math_noasm.go
17uint64_t add_mod_p(uint64_t a, uint64_t b)
18{
19 a = P - a;
20 uint64_t res = b - a;
21 if (b < a)
22 res += P;
23 return res;
24}
25
26uint64_t sub_mod_p(uint64_t a, uint64_t b)
27{
28 uint64_t r = a - b;
29 if (a < b)
30 r += P;
31 return r;
32}
33
34uint64_t reduce_mod_p(uint64_t b, uint64_t a)
35{
36 uint32_t d = b >> 32,
37 c = b;
38 if (a >= P)
39 a -= P;
40 a = sub_mod_p(a, c);
41 a = sub_mod_p(a, d);
42 a = add_mod_p(a, ((uint64_t)c) << 32);
43 return a;
44}
45
46uint64_t mult_mod_p(uint64_t x, uint64_t y)
47{
48 uint32_t a = x,
49 b = x >> 32,
50 c = y,
51 d = y >> 32;
52
53 /* first synthesize the product using 32*32 -> 64 bit multiplies */
54 x = b * (uint64_t)c; /* b*c */
55 y = a * (uint64_t)d; /* a*d */
56 uint64_t e = a * (uint64_t)c, /* a*c */
57 f = b * (uint64_t)d, /* b*d */
58 t;
59
60 x += y; /* b*c + a*d */
61 /* carry? */
62 if (x < y)
63 f += 1LL << 32; /* carry into upper 32 bits - can't overflow */
64
65 t = x << 32;
66 e += t; /* a*c + LSW(b*c + a*d) */
67 /* carry? */
68 if (e < t)
69 f += 1; /* carry into upper 64 bits - can't overflow*/
70 t = x >> 32;
71 f += t; /* b*d + MSW(b*c + a*d) */
72 /* can't overflow */
73
74 /* now reduce: (b*d + MSW(b*c + a*d), a*c + LSW(b*c + a*d)) */
75 return reduce_mod_p(f, e);
76}
77
78/* ================================================================================================
79 * RPO128 Permutation
80 */
81
82static const uint64_t STATE_WIDTH = 12;
83static const uint64_t NUM_ROUNDS = 7;
84
85/*
86 * MDS matrix
87 */
88const uint64_t MDS[12][12] = {
89 { 7, 23, 8, 26, 13, 10, 9, 7, 6, 22, 21, 8 },
90 { 8, 7, 23, 8, 26, 13, 10, 9, 7, 6, 22, 21 },
91 { 21, 8, 7, 23, 8, 26, 13, 10, 9, 7, 6, 22 },
92 { 22, 21, 8, 7, 23, 8, 26, 13, 10, 9, 7, 6 },
93 { 6, 22, 21, 8, 7, 23, 8, 26, 13, 10, 9, 7 },
94 { 7, 6, 22, 21, 8, 7, 23, 8, 26, 13, 10, 9 },
95 { 9, 7, 6, 22, 21, 8, 7, 23, 8, 26, 13, 10 },
96 { 10, 9, 7, 6, 22, 21, 8, 7, 23, 8, 26, 13 },
97 { 13, 10, 9, 7, 6, 22, 21, 8, 7, 23, 8, 26 },
98 { 26, 13, 10, 9, 7, 6, 22, 21, 8, 7, 23, 8 },
99 { 8, 26, 13, 10, 9, 7, 6, 22, 21, 8, 7, 23 },
100 { 23, 8, 26, 13, 10, 9, 7, 6, 22, 21, 8, 7 },
101};
102
103/*
104 * Round constants.
105 */
106const uint64_t ARK1[7][12] = {
107 {
108 5789762306288267392ULL,
109 6522564764413701783ULL,
110 17809893479458208203ULL,
111 107145243989736508ULL,
112 6388978042437517382ULL,
113 15844067734406016715ULL,
114 9975000513555218239ULL,
115 3344984123768313364ULL,
116 9959189626657347191ULL,
117 12960773468763563665ULL,
118 9602914297752488475ULL,
119 16657542370200465908ULL,
120 },
121 {
122 12987190162843096997ULL,
123 653957632802705281ULL,
124 4441654670647621225ULL,
125 4038207883745915761ULL,
126 5613464648874830118ULL,
127 13222989726778338773ULL,
128 3037761201230264149ULL,
129 16683759727265180203ULL,
130 8337364536491240715ULL,
131 3227397518293416448ULL,
132 8110510111539674682ULL,
133 2872078294163232137ULL,
134 },
135 {
136 18072785500942327487ULL,
137 6200974112677013481ULL,
138 17682092219085884187ULL,
139 10599526828986756440ULL,
140 975003873302957338ULL,
141 8264241093196931281ULL,
142 10065763900435475170ULL,
143 2181131744534710197ULL,
144 6317303992309418647ULL,
145 1401440938888741532ULL,
146 8884468225181997494ULL,
147 13066900325715521532ULL,
148 },
149 {
150 5674685213610121970ULL,
151 5759084860419474071ULL,
152 13943282657648897737ULL,
153 1352748651966375394ULL,
154 17110913224029905221ULL,
155 1003883795902368422ULL,
156 4141870621881018291ULL,
157 8121410972417424656ULL,
158 14300518605864919529ULL,
159 13712227150607670181ULL,
160 17021852944633065291ULL,
161 6252096473787587650ULL,
162 },
163 {
164 4887609836208846458ULL,
165 3027115137917284492ULL,
166 9595098600469470675ULL,
167 10528569829048484079ULL,
168 7864689113198939815ULL,
169 17533723827845969040ULL,
170 5781638039037710951ULL,
171 17024078752430719006ULL,
172 109659393484013511ULL,
173 7158933660534805869ULL,
174 2955076958026921730ULL,
175 7433723648458773977ULL,
176 },
177 {
178 16308865189192447297ULL,
179 11977192855656444890ULL,
180 12532242556065780287ULL,
181 14594890931430968898ULL,
182 7291784239689209784ULL,
183 5514718540551361949ULL,
184 10025733853830934803ULL,
185 7293794580341021693ULL,
186 6728552937464861756ULL,
187 6332385040983343262ULL,
188 13277683694236792804ULL,
189 2600778905124452676ULL,
190 },
191 {
192 7123075680859040534ULL,
193 1034205548717903090ULL,
194 7717824418247931797ULL,
195 3019070937878604058ULL,
196 11403792746066867460ULL,
197 10280580802233112374ULL,
198 337153209462421218ULL,
199 13333398568519923717ULL,
200 3596153696935337464ULL,
201 8104208463525993784ULL,
202 14345062289456085693ULL,
203 17036731477169661256ULL,
204 }};
205
206const uint64_t ARK2[7][12] = {
207 {
208 6077062762357204287ULL,
209 15277620170502011191ULL,
210 5358738125714196705ULL,
211 14233283787297595718ULL,
212 13792579614346651365ULL,
213 11614812331536767105ULL,
214 14871063686742261166ULL,
215 10148237148793043499ULL,
216 4457428952329675767ULL,
217 15590786458219172475ULL,
218 10063319113072092615ULL,
219 14200078843431360086ULL,
220 },
221 {
222 6202948458916099932ULL,
223 17690140365333231091ULL,
224 3595001575307484651ULL,
225 373995945117666487ULL,
226 1235734395091296013ULL,
227 14172757457833931602ULL,
228 707573103686350224ULL,
229 15453217512188187135ULL,
230 219777875004506018ULL,
231 17876696346199469008ULL,
232 17731621626449383378ULL,
233 2897136237748376248ULL,
234 },
235 {
236 8023374565629191455ULL,
237 15013690343205953430ULL,
238 4485500052507912973ULL,
239 12489737547229155153ULL,
240 9500452585969030576ULL,
241 2054001340201038870ULL,
242 12420704059284934186ULL,
243 355990932618543755ULL,
244 9071225051243523860ULL,
245 12766199826003448536ULL,
246 9045979173463556963ULL,
247 12934431667190679898ULL,
248 },
249 {
250 18389244934624494276ULL,
251 16731736864863925227ULL,
252 4440209734760478192ULL,
253 17208448209698888938ULL,
254 8739495587021565984ULL,
255 17000774922218161967ULL,
256 13533282547195532087ULL,
257 525402848358706231ULL,
258 16987541523062161972ULL,
259 5466806524462797102ULL,
260 14512769585918244983ULL,
261 10973956031244051118ULL,
262 },
263 {
264 6982293561042362913ULL,
265 14065426295947720331ULL,
266 16451845770444974180ULL,
267 7139138592091306727ULL,
268 9012006439959783127ULL,
269 14619614108529063361ULL,
270 1394813199588124371ULL,
271 4635111139507788575ULL,
272 16217473952264203365ULL,
273 10782018226466330683ULL,
274 6844229992533662050ULL,
275 7446486531695178711ULL,
276 },
277 {
278 3736792340494631448ULL,
279 577852220195055341ULL,
280 6689998335515779805ULL,
281 13886063479078013492ULL,
282 14358505101923202168ULL,
283 7744142531772274164ULL,
284 16135070735728404443ULL,
285 12290902521256031137ULL,
286 12059913662657709804ULL,
287 16456018495793751911ULL,
288 4571485474751953524ULL,
289 17200392109565783176ULL,
290 },
291 {
292 17130398059294018733ULL,
293 519782857322261988ULL,
294 9625384390925085478ULL,
295 1664893052631119222ULL,
296 7629576092524553570ULL,
297 3485239601103661425ULL,
298 9755891797164033838ULL,
299 15218148195153269027ULL,
300 16460604813734957368ULL,
301 9643968136937729763ULL,
302 3611348709641382851ULL,
303 18256379591337759196ULL,
304 },
305};
306
307void apply_sbox(uint64_t *const state)
308{
309 for (uint64_t i = 0; i < STATE_WIDTH; i++)
310 {
311 uint64_t t2 = mult_mod_p(*(state + i), *(state + i));
312 uint64_t t4 = mult_mod_p(t2, t2);
313
314 *(state + i) = mult_mod_p(*(state + i), mult_mod_p(t2, t4));
315 }
316}
317
318void apply_mds(uint64_t *state)
319{
320 uint64_t res[STATE_WIDTH];
321 for (uint64_t i = 0; i < STATE_WIDTH; i++)
322 {
323 res[i] = 0;
324 }
325 for (uint64_t i = 0; i < STATE_WIDTH; i++)
326 {
327 for (uint64_t j = 0; j < STATE_WIDTH; j++)
328 {
329 res[i] = add_mod_p(res[i], mult_mod_p(MDS[i][j], *(state + j)));
330 }
331 }
332
333 for (uint64_t i = 0; i < STATE_WIDTH; i++)
334 {
335 *(state + i) = res[i];
336 }
337}
338
339void apply_constants(uint64_t *const state, const uint64_t *ark)
340{
341 for (uint64_t i = 0; i < STATE_WIDTH; i++)
342 {
343 *(state + i) = add_mod_p(*(state + i), *(ark + i));
344 }
345}
346
347void exp_acc(const uint64_t m, const uint64_t *base, const uint64_t *tail, uint64_t *const res)
348{
349 for (uint64_t i = 0; i < m; i++)
350 {
351 for (uint64_t j = 0; j < STATE_WIDTH; j++)
352 {
353 if (i == 0)
354 {
355 *(res + j) = mult_mod_p(*(base + j), *(base + j));
356 }
357 else
358 {
359 *(res + j) = mult_mod_p(*(res + j), *(res + j));
360 }
361 }
362 }
363
364 for (uint64_t i = 0; i < STATE_WIDTH; i++)
365 {
366 *(res + i) = mult_mod_p(*(res + i), *(tail + i));
367 }
368}
369
370void apply_inv_sbox(uint64_t *const state)
371{
372 uint64_t t1[STATE_WIDTH];
373 for (uint64_t i = 0; i < STATE_WIDTH; i++)
374 {
375 t1[i] = 0;
376 }
377 for (uint64_t i = 0; i < STATE_WIDTH; i++)
378 {
379 t1[i] = mult_mod_p(*(state + i), *(state + i));
380 }
381
382 uint64_t t2[STATE_WIDTH];
383 for (uint64_t i = 0; i < STATE_WIDTH; i++)
384 {
385 t2[i] = 0;
386 }
387 for (uint64_t i = 0; i < STATE_WIDTH; i++)
388 {
389 t2[i] = mult_mod_p(t1[i], t1[i]);
390 }
391
392 uint64_t t3[STATE_WIDTH];
393 for (uint64_t i = 0; i < STATE_WIDTH; i++)
394 {
395 t3[i] = 0;
396 }
397 exp_acc(3, t2, t2, t3);
398
399 uint64_t t4[STATE_WIDTH];
400 for (uint64_t i = 0; i < STATE_WIDTH; i++)
401 {
402 t4[i] = 0;
403 }
404 exp_acc(6, t3, t3, t4);
405
406 uint64_t tmp[STATE_WIDTH];
407 for (uint64_t i = 0; i < STATE_WIDTH; i++)
408 {
409 tmp[i] = 0;
410 }
411 exp_acc(12, t4, t4, tmp);
412
413 uint64_t t5[STATE_WIDTH];
414 for (uint64_t i = 0; i < STATE_WIDTH; i++)
415 {
416 t5[i] = 0;
417 }
418 exp_acc(6, tmp, t3, t5);
419
420 uint64_t t6[STATE_WIDTH];
421 for (uint64_t i = 0; i < STATE_WIDTH; i++)
422 {
423 t6[i] = 0;
424 }
425 exp_acc(31, t5, t5, t6);
426
427 for (uint64_t i = 0; i < STATE_WIDTH; i++)
428 {
429 uint64_t a = mult_mod_p(mult_mod_p(t6[i], t6[i]), t5[i]);
430 a = mult_mod_p(a, a);
431 a = mult_mod_p(a, a);
432 uint64_t b = mult_mod_p(mult_mod_p(t1[i], t2[i]), *(state + i));
433
434 *(state + i) = mult_mod_p(a, b);
435 }
436}
437
438void apply_round(uint64_t *const state, const uint64_t round)
439{
440 apply_mds(state);
441 apply_constants(state, ARK1[round]);
442 apply_sbox(state);
443
444 apply_mds(state);
445 apply_constants(state, ARK2[round]);
446 apply_inv_sbox(state);
447}
448
449static void apply_permutation(uint64_t *state)
450{
451 for (uint64_t i = 0; i < NUM_ROUNDS; i++)
452 {
453 apply_round(state, i);
454 }
455}
456
457/* ================================================================================================
458 * RPO128 implementation. This is supposed to substitute SHAKE256 in the hash-to-point algorithm.
459 */
460
461#include "rpo.h"
462
463void rpo128_init(rpo128_context *rc)
464{
465 rc->dptr = 32;
466
467 memset(rc->st.A, 0, sizeof rc->st.A);
468}
469
470void rpo128_absorb(rpo128_context *rc, const uint8_t *in, size_t len)
471{
472 size_t dptr;
473
474 dptr = (size_t)rc->dptr;
475 while (len > 0)
476 {
477 size_t clen, u;
478
479 /* 136 * 8 = 1088 bit for the rate portion in the case of SHAKE256
480 * For RPO, this is 64 * 8 = 512 bits
481 * The capacity for SHAKE256 is at the end while for RPO128 it is at the beginning
482 */
483 clen = 96 - dptr;
484 if (clen > len)
485 {
486 clen = len;
487 }
488
489 for (u = 0; u < clen; u++)
490 {
491 rc->st.dbuf[dptr + u] = in[u];
492 }
493
494 dptr += clen;
495 in += clen;
496 len -= clen;
497 if (dptr == 96)
498 {
499 apply_permutation(rc->st.A);
500 dptr = 32;
501 }
502 }
503 rc->dptr = dptr;
504}
505
506void rpo128_finalize(rpo128_context *rc)
507{
508 // Set dptr to the end of the buffer, so that first call to extract will call the permutation.
509 rc->dptr = 96;
510}
511
512void rpo128_squeeze(rpo128_context *rc, uint8_t *out, size_t len)
513{
514 size_t dptr;
515
516 dptr = (size_t)rc->dptr;
517 while (len > 0)
518 {
519 size_t clen;
520
521 if (dptr == 96)
522 {
523 apply_permutation(rc->st.A);
524 dptr = 32;
525 }
526 clen = 96 - dptr;
527 if (clen > len)
528 {
529 clen = len;
530 }
531 len -= clen;
532
533 memcpy(out, rc->st.dbuf + dptr, clen);
534 dptr += clen;
535 out += clen;
536 }
537 rc->dptr = dptr;
538}
539
540void rpo128_release(rpo128_context *rc)
541{
542 memset(rc->st.A, 0, sizeof rc->st.A);
543 rc->dptr = 32;
544}
545
546/* ================================================================================================
547 * Hash-to-Point algorithm implementation based on RPO128
548 */
549
550void PQCLEAN_FALCON512_CLEAN_hash_to_point_rpo(rpo128_context *rc, uint16_t *x, unsigned logn)
551{
552 /*
553 * This implementation avoids the rejection sampling step needed in the
554 * per-the-spec implementation. It uses a remark in https://falcon-sign.info/falcon.pdf
555 * page 31, which argues that the current variant is secure for the parameters set by NIST.
556 * Avoiding the rejection-sampling step leads to an implementation that is constant-time.
557 * TODO: Check that the current implementation is indeed constant-time.
558 */
559 size_t n;
560
561 n = (size_t)1 << logn;
562 while (n > 0)
563 {
564 uint8_t buf[8];
565 uint64_t w;
566
567 rpo128_squeeze(rc, (void *)buf, sizeof buf);
568 w = ((uint64_t)(buf[7]) << 56) |
569 ((uint64_t)(buf[6]) << 48) |
570 ((uint64_t)(buf[5]) << 40) |
571 ((uint64_t)(buf[4]) << 32) |
572 ((uint64_t)(buf[3]) << 24) |
573 ((uint64_t)(buf[2]) << 16) |
574 ((uint64_t)(buf[1]) << 8) |
575 ((uint64_t)(buf[0]));
576
577 w %= M;
578
579 *x++ = (uint16_t)w;
580 n--;
581 }
582}