Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_crypto/src/c/poly_mul.c

7.4 KiB, 1 run

created by r1870400018:216, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1#include "poly_mul.h"
2#include <stdint.h>
3#include <string.h>
4
5#define SCHB_N 16
6
7#define N_RES (SABER_N << 1)
8#define N_SB (SABER_N >> 2)
9#define N_SB_RES (2*N_SB-1)
10
11#define OVERFLOWING_MUL(X, Y) ((uint16_t)((uint32_t)(X) * (uint32_t)(Y)))
12
13#define KARATSUBA_N 64
14static void karatsuba_simple(const uint16_t *a_1, const uint16_t *b_1, uint16_t *result_final) {
15 uint16_t d01[KARATSUBA_N / 2 - 1];
16 uint16_t d0123[KARATSUBA_N / 2 - 1];
17 uint16_t d23[KARATSUBA_N / 2 - 1];
18 uint16_t result_d01[KARATSUBA_N - 1];
19
20 int32_t i, j;
21
22 memset(result_d01, 0, (KARATSUBA_N - 1)*sizeof(uint16_t));
23 memset(d01, 0, (KARATSUBA_N / 2 - 1)*sizeof(uint16_t));
24 memset(d0123, 0, (KARATSUBA_N / 2 - 1)*sizeof(uint16_t));
25 memset(d23, 0, (KARATSUBA_N / 2 - 1)*sizeof(uint16_t));
26 memset(result_final, 0, (2 * KARATSUBA_N - 1)*sizeof(uint16_t));
27
28 uint16_t acc1, acc2, acc3, acc4, acc5, acc6, acc7, acc8, acc9, acc10;
29
30
31 for (i = 0; i < KARATSUBA_N / 4; i++) {
32 acc1 = a_1[i]; //a0
33 acc2 = a_1[i + KARATSUBA_N / 4]; //a1
34 acc3 = a_1[i + 2 * KARATSUBA_N / 4]; //a2
35 acc4 = a_1[i + 3 * KARATSUBA_N / 4]; //a3
36 for (j = 0; j < KARATSUBA_N / 4; j++) {
37
38 acc5 = b_1[j]; //b0
39 acc6 = b_1[j + KARATSUBA_N / 4]; //b1
40
41 result_final[i + j + 0 * KARATSUBA_N / 4] =
42 result_final[i + j + 0 * KARATSUBA_N / 4] +
43 OVERFLOWING_MUL(acc1, acc5);
44 result_final[i + j + 2 * KARATSUBA_N / 4] =
45 result_final[i + j + 2 * KARATSUBA_N / 4] +
46 OVERFLOWING_MUL(acc2, acc6);
47
48 acc7 = acc5 + acc6; //b01
49 acc8 = acc1 + acc2; //a01
50 d01[i + j] = d01[i + j] + (uint16_t)(acc7 * (uint64_t)acc8);
51 //--------------------------------------------------------
52
53 acc7 = b_1[j + 2 * KARATSUBA_N / 4]; //b2
54 acc8 = b_1[j + 3 * KARATSUBA_N / 4]; //b3
55 result_final[i + j + 4 * KARATSUBA_N / 4] =
56 result_final[i + j + 4 * KARATSUBA_N / 4] +
57 OVERFLOWING_MUL(acc7, acc3);
58
59 result_final[i + j + 6 * KARATSUBA_N / 4] =
60 result_final[i + j + 6 * KARATSUBA_N / 4] +
61 OVERFLOWING_MUL(acc8, acc4);
62
63 acc9 = acc3 + acc4;
64 acc10 = acc7 + acc8;
65 d23[i + j] = d23[i + j] + OVERFLOWING_MUL(acc9, acc10);
66 //--------------------------------------------------------
67
68 acc5 = acc5 + acc7; //b02
69 acc7 = acc1 + acc3; //a02
70 result_d01[i + j + 0 * KARATSUBA_N / 4] =
71 result_d01[i + j + 0 * KARATSUBA_N / 4] +
72 OVERFLOWING_MUL(acc5, acc7);
73
74 acc6 = acc6 + acc8; //b13
75 acc8 = acc2 + acc4;
76 result_d01[i + j + 2 * KARATSUBA_N / 4] =
77 result_d01[i + j + 2 * KARATSUBA_N / 4] +
78 OVERFLOWING_MUL(acc6, acc8);
79
80 acc5 = acc5 + acc6;
81 acc7 = acc7 + acc8;
82 d0123[i + j] = d0123[i + j] + OVERFLOWING_MUL(acc5, acc7);
83 }
84 }
85
86 // 2nd last stage
87
88 for (i = 0; i < KARATSUBA_N / 2 - 1; i++) {
89 d0123[i] = d0123[i] - result_d01[i + 0 * KARATSUBA_N / 4] - result_d01[i + 2 * KARATSUBA_N / 4];
90 d01[i] = d01[i] - result_final[i + 0 * KARATSUBA_N / 4] - result_final[i + 2 * KARATSUBA_N / 4];
91 d23[i] = d23[i] - result_final[i + 4 * KARATSUBA_N / 4] - result_final[i + 6 * KARATSUBA_N / 4];
92 }
93
94 for (i = 0; i < KARATSUBA_N / 2 - 1; i++) {
95 result_d01[i + 1 * KARATSUBA_N / 4] = result_d01[i + 1 * KARATSUBA_N / 4] + d0123[i];
96 result_final[i + 1 * KARATSUBA_N / 4] = result_final[i + 1 * KARATSUBA_N / 4] + d01[i];
97 result_final[i + 5 * KARATSUBA_N / 4] = result_final[i + 5 * KARATSUBA_N / 4] + d23[i];
98 }
99
100 // Last stage
101 for (i = 0; i < KARATSUBA_N - 1; i++) {
102 result_d01[i] = result_d01[i] - result_final[i] - result_final[i + KARATSUBA_N];
103 }
104
105 for (i = 0; i < KARATSUBA_N - 1; i++) {
106 result_final[i + 1 * KARATSUBA_N / 2] = result_final[i + 1 * KARATSUBA_N / 2] + result_d01[i];
107 }
108
109}
110
111
112
113static void toom_cook_4way (const uint16_t *a1, const uint16_t *b1, uint16_t *result) {
114 uint16_t inv3 = 43691, inv9 = 36409, inv15 = 61167;
115
116 uint16_t aw1[N_SB], aw2[N_SB], aw3[N_SB], aw4[N_SB], aw5[N_SB], aw6[N_SB], aw7[N_SB];
117 uint16_t bw1[N_SB], bw2[N_SB], bw3[N_SB], bw4[N_SB], bw5[N_SB], bw6[N_SB], bw7[N_SB];
118 uint16_t w1[N_SB_RES] = {0}, w2[N_SB_RES] = {0}, w3[N_SB_RES] = {0}, w4[N_SB_RES] = {0},
119 w5[N_SB_RES] = {0}, w6[N_SB_RES] = {0}, w7[N_SB_RES] = {0};
120 uint16_t r0, r1, r2, r3, r4, r5, r6, r7;
121 uint16_t *A0, *A1, *A2, *A3, *B0, *B1, *B2, *B3;
122 A0 = (uint16_t *)a1;
123 A1 = (uint16_t *)&a1[N_SB];
124 A2 = (uint16_t *)&a1[2 * N_SB];
125 A3 = (uint16_t *)&a1[3 * N_SB];
126 B0 = (uint16_t *)b1;
127 B1 = (uint16_t *)&b1[N_SB];
128 B2 = (uint16_t *)&b1[2 * N_SB];
129 B3 = (uint16_t *)&b1[3 * N_SB];
130
131 uint16_t *C;
132 C = result;
133
134 int i, j;
135
136 // EVALUATION
137 for (j = 0; j < N_SB; ++j) {
138 r0 = A0[j];
139 r1 = A1[j];
140 r2 = A2[j];
141 r3 = A3[j];
142 r4 = r0 + r2;
143 r5 = r1 + r3;
144 r6 = r4 + r5;
145 r7 = r4 - r5;
146 aw3[j] = r6;
147 aw4[j] = r7;
148 r4 = ((r0 << 2) + r2) << 1;
149 r5 = (r1 << 2) + r3;
150 r6 = r4 + r5;
151 r7 = r4 - r5;
152 aw5[j] = r6;
153 aw6[j] = r7;
154 r4 = (r3 << 3) + (r2 << 2) + (r1 << 1) + r0;
155 aw2[j] = r4;
156 aw7[j] = r0;
157 aw1[j] = r3;
158 }
159 for (j = 0; j < N_SB; ++j) {
160 r0 = B0[j];
161 r1 = B1[j];
162 r2 = B2[j];
163 r3 = B3[j];
164 r4 = r0 + r2;
165 r5 = r1 + r3;
166 r6 = r4 + r5;
167 r7 = r4 - r5;
168 bw3[j] = r6;
169 bw4[j] = r7;
170 r4 = ((r0 << 2) + r2) << 1;
171 r5 = (r1 << 2) + r3;
172 r6 = r4 + r5;
173 r7 = r4 - r5;
174 bw5[j] = r6;
175 bw6[j] = r7;
176 r4 = (r3 << 3) + (r2 << 2) + (r1 << 1) + r0;
177 bw2[j] = r4;
178 bw7[j] = r0;
179 bw1[j] = r3;
180 }
181
182 // MULTIPLICATION
183
184 karatsuba_simple(aw1, bw1, w1);
185 karatsuba_simple(aw2, bw2, w2);
186 karatsuba_simple(aw3, bw3, w3);
187 karatsuba_simple(aw4, bw4, w4);
188 karatsuba_simple(aw5, bw5, w5);
189 karatsuba_simple(aw6, bw6, w6);
190 karatsuba_simple(aw7, bw7, w7);
191
192 // INTERPOLATION
193 for (i = 0; i < N_SB_RES; ++i) {
194 r0 = w1[i];
195 r1 = w2[i];
196 r2 = w3[i];
197 r3 = w4[i];
198 r4 = w5[i];
199 r5 = w6[i];
200 r6 = w7[i];
201
202 r1 = r1 + r4;
203 r5 = r5 - r4;
204 r3 = ((r3 - r2) >> 1);
205 r4 = r4 - r0;
206 r4 = r4 - (r6 << 6);
207 r4 = (r4 << 1) + r5;
208 r2 = r2 + r3;
209 r1 = r1 - (r2 << 6) - r2;
210 r2 = r2 - r6;
211 r2 = r2 - r0;
212 r1 = r1 + 45 * r2;
213 r4 = (uint16_t)(((r4 - (r2 << 3)) * (uint32_t)inv3) >> 3);
214 r5 = r5 + r1;
215 r1 = (uint16_t)(((r1 + (r3 << 4)) * (uint32_t)inv9) >> 1);
216 r3 = -(r3 + r1);
217 r5 = (uint16_t)(((30 * r1 - r5) * (uint32_t)inv15) >> 2);
218 r2 = r2 - r4;
219 r1 = r1 - r5;
220
221 C[i] += r6;
222 C[i + 64] += r5;
223 C[i + 128] += r4;
224 C[i + 192] += r3;
225 C[i + 256] += r2;
226 C[i + 320] += r1;
227 C[i + 384] += r0;
228 }
229}
230
231/* res += a*b */
232void poly_mul_acc(const uint16_t a[SABER_N], const uint16_t b[SABER_N], uint16_t res[SABER_N])
233{
234 uint16_t c[2 * SABER_N] = {0};
235 int i;
236
237 toom_cook_4way(a, b, c);
238
239 /* reduction */
240 for (i = SABER_N; i < 2 * SABER_N; i++)
241 {
242 res[i - SABER_N] += (c[i - SABER_N] - c[i]);
243 }
244
245}