poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/q120/q120_arithmetic_ref.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/q120/q120_arithmetic_ref.cRaw
1#include <assert.h>
2#include <math.h>
3#include <stdint.h>
4
5#include "q120_arithmetic.h"
6#include "q120_arithmetic_private.h"
7#include "q120_common.h"
8
9#define MODQ(val, q) ((val) % (q))
10
11double comp_bit_size_red(const uint64_t h, const uint64_t qs[4]) {
12 assert(h < 128);
13 double h_pow2_bs = 0;
14 for (uint64_t k = 0; k < 4; ++k) {
15 double t = log2((double)MODQ((__uint128_t)1 << h, qs[k]));
16 if (t > h_pow2_bs) h_pow2_bs = t;
17 }
18 return h_pow2_bs;
19}
20
21double comp_bit_size_sum(const uint64_t n, const double* const bs) {
22 double s = 0;
23 for (uint64_t i = 0; i < n; ++i) {
24 s += pow(2, bs[i]);
25 }
26 return log2(s);
27}
28
29void vec_mat1col_product_baa_precomp(q120_mat1col_product_baa_precomp* precomp) {
30 uint64_t qs[4] = {Q1, Q2, Q3, Q4};
31
32 double min_res_bs = 1000;
33 uint64_t min_h = -1;
34
35 double ell_bs = log2((double)MAX_ELL);
36 for (uint64_t h = 1; h < 64; ++h) {
37 double h_pow2_bs = comp_bit_size_red(h, qs);
38
39 const double bs[] = {h + ell_bs, 64 - h + ell_bs + h_pow2_bs};
40 const double res_bs = comp_bit_size_sum(2, bs);
41
42 if (min_res_bs > res_bs) {
43 min_res_bs = res_bs;
44 min_h = h;
45 }
46 }
47
48 assert(min_res_bs < 64);
49 precomp->h = min_h;
50 for (uint64_t k = 0; k < 4; ++k) {
51 precomp->h_pow_red[k] = MODQ(UINT64_C(1) << precomp->h, qs[k]);
52 }
53#ifndef NDEBUG
54 precomp->res_bit_size = min_res_bs;
55#endif
56 // printf("AA %lu %lf\n", min_h, min_res_bs);
57}
58
59EXPORT q120_mat1col_product_baa_precomp* q120_new_vec_mat1col_product_baa_precomp() {
60 q120_mat1col_product_baa_precomp* res = malloc(sizeof(q120_mat1col_product_baa_precomp));
61 vec_mat1col_product_baa_precomp(res);
62 return res;
63}
64
65EXPORT void q120_delete_vec_mat1col_product_baa_precomp(q120_mat1col_product_baa_precomp* addr) { free(addr); }
66
67EXPORT void q120_vec_mat1col_product_baa_ref(q120_mat1col_product_baa_precomp* precomp, const uint64_t ell,
68 q120b* const res, const q120a* const x, const q120a* const y) {
69 /**
70 * Algorithm:
71 * - res = acc1 + acc2 . ((2^H) % Q)
72 * - acc1 is the sum of H LSB of products x[i].y[i]
73 * - acc2 is the sum of 64-H MSB of products x[i]].y[i]
74 * - for l < 10k acc1 will have H + log2(10000) and acc2 64 - H + log2(10000) bits
75 * - final sum has max(H, 64 - H + bit_size((2^H) % Q)) + log2(10000) + 1 bits
76 */
77 const uint64_t H = precomp->h;
78 const uint64_t MASK = (UINT64_C(1) << H) - 1;
79
80 uint64_t acc1[4] = {0, 0, 0, 0}; // accumulate H least significant bits of product
81 uint64_t acc2[4] = {0, 0, 0, 0}; // accumulate 64 - H most significan bits of product
82
83 const uint64_t* const x_ptr = (uint64_t*)x;
84 const uint64_t* const y_ptr = (uint64_t*)y;
85
86 for (uint64_t i = 0; i < 4 * ell; i += 4) {
87 for (uint64_t j = 0; j < 4; ++j) {
88 uint64_t t = x_ptr[i + j] * y_ptr[i + j];
89 acc1[j] += t & MASK;
90 acc2[j] += t >> H;
91 }
92 }
93
94 uint64_t* const res_ptr = (uint64_t*)res;
95 for (uint64_t j = 0; j < 4; ++j) {
96 res_ptr[j] = acc1[j] + acc2[j] * precomp->h_pow_red[j];
97 assert(log2(res_ptr[j]) < precomp->res_bit_size);
98 }
99}
100
101void vec_mat1col_product_bbb_precomp(q120_mat1col_product_bbb_precomp* precomp) {
102 uint64_t qs[4] = {Q1, Q2, Q3, Q4};
103
104 double ell_bs = log2((double)MAX_ELL);
105 double min_res_bs = 1000;
106 uint64_t min_h = -1;
107
108 const double s1_bs = 32 + ell_bs;
109 const double s2_bs = 32 + ell_bs + log2(3);
110 const double s3_bs = 32 + ell_bs + log2(3);
111 const double s4_bs = 32 + ell_bs;
112 for (uint64_t h = 16; h < 32; ++h) {
113 const double s1l_bs = h;
114 const double s1h_bs = (s1_bs - h) + comp_bit_size_red(h, qs);
115 const double s2l_bs = h + comp_bit_size_red(32, qs);
116 const double s2h_bs = (s2_bs - h) + comp_bit_size_red(32 + h, qs);
117 const double s3l_bs = h + comp_bit_size_red(64, qs);
118 const double s3h_bs = (s3_bs - h) + comp_bit_size_red(64 + h, qs);
119 const double s4l_bs = h + comp_bit_size_red(96, qs);
120 const double s4h_bs = (s4_bs - h) + comp_bit_size_red(96 + h, qs);
121
122 const double bs[] = {s1l_bs, s1h_bs, s2l_bs, s2h_bs, s3l_bs, s3h_bs, s4l_bs, s4h_bs};
123 const double res_bs = comp_bit_size_sum(8, bs);
124
125 if (min_res_bs > res_bs) {
126 min_res_bs = res_bs;
127 min_h = h;
128 }
129 }
130
131 assert(min_res_bs < 64);
132 precomp->h = min_h;
133 for (uint64_t k = 0; k < 4; ++k) {
134 precomp->s1h_pow_red[k] = UINT64_C(1) << precomp->h; // 2^24
135 precomp->s2l_pow_red[k] = MODQ(UINT64_C(1) << 32, qs[k]); // 2^32
136 precomp->s2h_pow_red[k] = MODQ(precomp->s2l_pow_red[k] * precomp->s1h_pow_red[k], qs[k]); // 2^(32+24)
137 precomp->s3l_pow_red[k] = MODQ(precomp->s2l_pow_red[k] * precomp->s2l_pow_red[k], qs[k]); // 2^64 = 2^(32+32)
138 precomp->s3h_pow_red[k] = MODQ(precomp->s3l_pow_red[k] * precomp->s1h_pow_red[k], qs[k]); // 2^(64+24)
139 precomp->s4l_pow_red[k] = MODQ(precomp->s3l_pow_red[k] * precomp->s2l_pow_red[k], qs[k]); // 2^96 = 2^(64+32)
140 precomp->s4h_pow_red[k] = MODQ(precomp->s4l_pow_red[k] * precomp->s1h_pow_red[k], qs[k]); // 2^(96+24)
141 }
142// printf("AA %lu %lf\n", min_h, min_res_bs);
143#ifndef NDEBUG
144 precomp->res_bit_size = min_res_bs;
145#endif
146}
147
148EXPORT q120_mat1col_product_bbb_precomp* q120_new_vec_mat1col_product_bbb_precomp() {
149 q120_mat1col_product_bbb_precomp* res = malloc(sizeof(q120_mat1col_product_bbb_precomp));
150 vec_mat1col_product_bbb_precomp(res);
151 return res;
152}
153
154EXPORT void q120_delete_vec_mat1col_product_bbb_precomp(q120_mat1col_product_bbb_precomp* addr) { free(addr); }
155
156EXPORT void q120_vec_mat1col_product_bbb_ref(q120_mat1col_product_bbb_precomp* precomp, const uint64_t ell,
157 q120b* const res, const q120b* const x, const q120b* const y) {
158 /**
159 * Algorithm:
160 * 1. Split x_i and y_i in 2 32-bit parts and compute the cross-products:
161 * - x_i = xl_i + xh_i . 2^32
162 * - y_i = yl_i + yh_i . 2^32
163 * - A_i = xl_i . yl_i
164 * - B_i = xl_i . yh_i
165 * - C_i = xh_i . yl_i
166 * - D_i = xh_i . yh_i
167 * - we have x_i . y_i == A_i + (B_i + C_i) . 2^32 + D_i . 2^64
168 * 2. Split A_i, B_i, C_i and D_i into 2 32-bit parts
169 * - A_i = Al_i + Ah_i . 2^32
170 * - B_i = Bl_i + Bh_i . 2^32
171 * - C_i = Cl_i + Ch_i . 2^32
172 * - D_i = Dl_i + Dh_i . 2^32
173 * 3. Compute the sums:
174 * - S1 = \sum Al_i
175 * - S2 = \sum (Ah_i + Bl_i + Cl_i)
176 * - S3 = \sum (Bh_i + Ch_i + Dl_i)
177 * - S4 = \sum Dh_i
178 * - here S1, S4 have 32 + log2(ell) bits and S2, S3 have 32 + log2(ell) +
179 * log2(3) bits
180 * - for ell == 10000 S2, S3 have < 47 bits
181 * 4. Split S1, S2, S3 and S4 in 2 24-bit parts (24 = ceil(47/2))
182 * - S1 = S1l + S1h . 2^24
183 * - S2 = S2l + S2h . 2^24
184 * - S3 = S3l + S3h . 2^24
185 * - S4 = S4l + S4h . 2^24
186 * 5. Compute final result as:
187 * - \sum x_i . y_i = S1l + S1h . 2^24
188 * + S2l . 2^32 + S2h . 2^(32+24)
189 * + S3l . 2^64 + S3h . 2^(64 + 24)
190 * + S4l . 2^96 + S4l . 2^(96+24)
191 * - here the powers of 2 are reduced modulo the primes Q before
192 * multiplications
193 * - the result will be on 24 + 3 + bit size of primes Q
194 */
195 const uint64_t H1 = 32;
196 const uint64_t MASK1 = (UINT64_C(1) << H1) - 1;
197
198 uint64_t s1[4] = {0, 0, 0, 0};
199 uint64_t s2[4] = {0, 0, 0, 0};
200 uint64_t s3[4] = {0, 0, 0, 0};
201 uint64_t s4[4] = {0, 0, 0, 0};
202
203 const uint64_t* const x_ptr = (uint64_t*)x;
204 const uint64_t* const y_ptr = (uint64_t*)y;
205
206 for (uint64_t i = 0; i < 4 * ell; i += 4) {
207 for (uint64_t j = 0; j < 4; ++j) {
208 const uint64_t xl = x_ptr[i + j] & MASK1;
209 const uint64_t xh = x_ptr[i + j] >> H1;
210 const uint64_t yl = y_ptr[i + j] & MASK1;
211 const uint64_t yh = y_ptr[i + j] >> H1;
212
213 const uint64_t a = xl * yl;
214 const uint64_t al = a & MASK1;
215 const uint64_t ah = a >> H1;
216
217 const uint64_t b = xl * yh;
218 const uint64_t bl = b & MASK1;
219 const uint64_t bh = b >> H1;
220
221 const uint64_t c = xh * yl;
222 const uint64_t cl = c & MASK1;
223 const uint64_t ch = c >> H1;
224
225 const uint64_t d = xh * yh;
226 const uint64_t dl = d & MASK1;
227 const uint64_t dh = d >> H1;
228
229 s1[j] += al;
230 s2[j] += ah + bl + cl;
231 s3[j] += bh + ch + dl;
232 s4[j] += dh;
233 }
234 }
235
236 const uint64_t H2 = precomp->h;
237 const uint64_t MASK2 = (UINT64_C(1) << H2) - 1;
238
239 uint64_t* const res_ptr = (uint64_t*)res;
240 for (uint64_t j = 0; j < 4; ++j) {
241 const uint64_t s1l = s1[j] & MASK2;
242 const uint64_t s1h = s1[j] >> H2;
243 const uint64_t s2l = s2[j] & MASK2;
244 const uint64_t s2h = s2[j] >> H2;
245 const uint64_t s3l = s3[j] & MASK2;
246 const uint64_t s3h = s3[j] >> H2;
247 const uint64_t s4l = s4[j] & MASK2;
248 const uint64_t s4h = s4[j] >> H2;
249
250 uint64_t t = s1l;
251 t += s1h * precomp->s1h_pow_red[j];
252 t += s2l * precomp->s2l_pow_red[j];
253 t += s2h * precomp->s2h_pow_red[j];
254 t += s3l * precomp->s3l_pow_red[j];
255 t += s3h * precomp->s3h_pow_red[j];
256 t += s4l * precomp->s4l_pow_red[j];
257 t += s4h * precomp->s4h_pow_red[j];
258
259 res_ptr[j] = t;
260 assert(log2(res_ptr[j]) < precomp->res_bit_size);
261 }
262}
263
264void vec_mat1col_product_bbc_precomp(q120_mat1col_product_bbc_precomp* precomp) {
265 uint64_t qs[4] = {Q1, Q2, Q3, Q4};
266
267 double min_res_bs = 1000;
268 uint64_t min_h = -1;
269
270 double pow2_32_bs = comp_bit_size_red(32, qs);
271
272 double ell_bs = log2((double)MAX_ELL);
273 double s1_bs = 32 + ell_bs;
274 for (uint64_t h = 16; h < 32; ++h) {
275 double s2l_bs = pow2_32_bs + h;
276 double s2h_bs = s1_bs - h + comp_bit_size_red(32 + h, qs);
277
278 const double bs[] = {s1_bs, s2l_bs, s2h_bs};
279 const double res_bs = comp_bit_size_sum(3, bs);
280
281 if (min_res_bs > res_bs) {
282 min_res_bs = res_bs;
283 min_h = h;
284 }
285 }
286
287 assert(min_res_bs < 64);
288 precomp->h = min_h;
289 for (uint64_t k = 0; k < 4; ++k) {
290 precomp->s2l_pow_red[k] = MODQ(UINT64_C(1) << 32, qs[k]);
291 precomp->s2h_pow_red[k] = MODQ(UINT64_C(1) << (32 + precomp->h), qs[k]);
292 }
293#ifndef NDEBUG
294 precomp->res_bit_size = min_res_bs;
295#endif
296 // printf("AA %lu %lf\n", min_h, min_res_bs);
297}
298
299EXPORT q120_mat1col_product_bbc_precomp* q120_new_vec_mat1col_product_bbc_precomp() {
300 q120_mat1col_product_bbc_precomp* res = malloc(sizeof(q120_mat1col_product_bbc_precomp));
301 vec_mat1col_product_bbc_precomp(res);
302 return res;
303}
304
305EXPORT void q120_delete_vec_mat1col_product_bbc_precomp(q120_mat1col_product_bbc_precomp* addr) { free(addr); }
306
307EXPORT void q120_vec_mat1col_product_bbc_ref_old(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
308 q120b* const res, const q120b* const x, const q120c* const y) {
309 /**
310 * Algorithm:
311 * 0. We have
312 * - y0_i == y_i % Q and y1_i == (y_i . 2^32) % Q
313 * 1. Split x_i in 2 32-bit parts and compute the cross-products:
314 * - x_i = xl_i + xh_i . 2^32
315 * - A_i = xl_i . y1_i
316 * - B_i = xh_i . y2_i
317 * - we have x_i . y_i == A_i + B_i
318 * 2. Split A_i and B_i into 2 32-bit parts
319 * - A_i = Al_i + Ah_i . 2^32
320 * - B_i = Bl_i + Bh_i . 2^32
321 * 3. Compute the sums:
322 * - S1 = \sum Al_i + Bl_i
323 * - S2 = \sum Ah_i + Bh_i
324 * - here S1 and S2 have 32 + log2(ell) bits
325 * - for ell == 10000 S1, S2 have < 46 bits
326 * 4. Split S2 in 27-bit and 19-bit parts (27+19 == 46)
327 * - S2 = S2l + S2h . 2^27
328 * 5. Compute final result as:
329 * - \sum x_i . y_i = S1 + S2l . 2^32 + S2h . 2^(32+27)
330 * - here the powers of 2 are reduced modulo the primes Q before
331 * multiplications
332 * - the result will be on < 52 bits
333 */
334
335 const uint64_t H1 = 32;
336 const uint64_t MASK1 = (UINT64_C(1) << H1) - 1;
337
338 uint64_t s1[4] = {0, 0, 0, 0};
339 uint64_t s2[4] = {0, 0, 0, 0};
340
341 const uint64_t* const x_ptr = (uint64_t*)x;
342 const uint32_t* const y_ptr = (uint32_t*)y;
343
344 for (uint64_t i = 0; i < 4 * ell; i += 4) {
345 for (uint64_t j = 0; j < 4; ++j) {
346 const uint64_t xl = x_ptr[i + j] & MASK1;
347 const uint64_t xh = x_ptr[i + j] >> H1;
348 const uint64_t y0 = y_ptr[2 * (i + j)];
349 const uint64_t y1 = y_ptr[2 * (i + j) + 1];
350
351 const uint64_t a = xl * y0;
352 const uint64_t al = a & MASK1;
353 const uint64_t ah = a >> H1;
354
355 const uint64_t b = xh * y1;
356 const uint64_t bl = b & MASK1;
357 const uint64_t bh = b >> H1;
358
359 s1[j] += al + bl;
360 s2[j] += ah + bh;
361 }
362 }
363
364 const uint64_t H2 = precomp->h;
365 const uint64_t MASK2 = (UINT64_C(1) << H2) - 1;
366
367 uint64_t* const res_ptr = (uint64_t*)res;
368 for (uint64_t k = 0; k < 4; ++k) {
369 const uint64_t s2l = s2[k] & MASK2;
370 const uint64_t s2h = s2[k] >> H2;
371
372 uint64_t t = s1[k];
373 t += s2l * precomp->s2l_pow_red[k];
374 t += s2h * precomp->s2h_pow_red[k];
375
376 res_ptr[k] = t;
377 assert(log2(res_ptr[k]) < precomp->res_bit_size);
378 }
379}
380
381static __always_inline void accum_mul_q120_bc(uint64_t res[8], //
382 const uint32_t x_layb[8], const uint32_t y_layc[8]) {
383 for (uint64_t i = 0; i < 4; ++i) {
384 static const uint64_t MASK32 = 0xFFFFFFFFUL;
385 uint64_t x_lo = x_layb[2 * i];
386 uint64_t x_hi = x_layb[2 * i + 1];
387 uint64_t y_lo = y_layc[2 * i];
388 uint64_t y_hi = y_layc[2 * i + 1];
389 uint64_t xy_lo = x_lo * y_lo;
390 uint64_t xy_hi = x_hi * y_hi;
391 res[2 * i] += (xy_lo & MASK32) + (xy_hi & MASK32);
392 res[2 * i + 1] += (xy_lo >> 32) + (xy_hi >> 32);
393 }
394}
395
396static __always_inline void accum_to_q120b(uint64_t res[4], //
397 const uint64_t s[8], const q120_mat1col_product_bbc_precomp* precomp) {
398 const uint64_t H2 = precomp->h;
399 const uint64_t MASK2 = (UINT64_C(1) << H2) - 1;
400 for (uint64_t k = 0; k < 4; ++k) {
401 const uint64_t s2l = s[2 * k + 1] & MASK2;
402 const uint64_t s2h = s[2 * k + 1] >> H2;
403 uint64_t t = s[2 * k];
404 t += s2l * precomp->s2l_pow_red[k];
405 t += s2h * precomp->s2h_pow_red[k];
406 res[k] = t;
407 assert(log2(res[k]) < precomp->res_bit_size);
408 }
409}
410
411EXPORT void q120_vec_mat1col_product_bbc_ref(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
412 q120b* const res, const q120b* const x, const q120c* const y) {
413 uint64_t s[8] = {0, 0, 0, 0, 0, 0, 0, 0};
414
415 const uint32_t(*const x_ptr)[8] = (const uint32_t(*const)[8])x;
416 const uint32_t(*const y_ptr)[8] = (const uint32_t(*const)[8])y;
417
418 for (uint64_t i = 0; i < ell; i++) {
419 accum_mul_q120_bc(s, x_ptr[i], y_ptr[i]);
420 }
421 accum_to_q120b((uint64_t*)res, s, precomp);
422}
423
424EXPORT void q120x2_vec_mat1col_product_bbc_ref(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
425 q120b* const res, const q120b* const x, const q120c* const y) {
426 uint64_t s[2][16] = {0};
427
428 const uint32_t(*const x_ptr)[2][8] = (const uint32_t(*const)[2][8])x;
429 const uint32_t(*const y_ptr)[2][8] = (const uint32_t(*const)[2][8])y;
430 uint64_t(*re)[4] = (uint64_t(*)[4])res;
431
432 for (uint64_t i = 0; i < ell; i++) {
433 accum_mul_q120_bc(s[0], x_ptr[i][0], y_ptr[i][0]);
434 accum_mul_q120_bc(s[1], x_ptr[i][1], y_ptr[i][1]);
435 }
436 accum_to_q120b(re[0], s[0], precomp);
437 accum_to_q120b(re[1], s[1], precomp);
438}
439
440EXPORT void q120x2_vec_mat2cols_product_bbc_ref(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
441 q120b* const res, const q120b* const x, const q120c* const y) {
442 uint64_t s[4][16] = {0};
443
444 const uint32_t(*const x_ptr)[2][8] = (const uint32_t(*const)[2][8])x;
445 const uint32_t(*const y_ptr)[4][8] = (const uint32_t(*const)[4][8])y;
446 uint64_t(*re)[4] = (uint64_t(*)[4])res;
447
448 for (uint64_t i = 0; i < ell; i++) {
449 accum_mul_q120_bc(s[0], x_ptr[i][0], y_ptr[i][0]);
450 accum_mul_q120_bc(s[1], x_ptr[i][1], y_ptr[i][1]);
451 accum_mul_q120_bc(s[2], x_ptr[i][0], y_ptr[i][2]);
452 accum_mul_q120_bc(s[3], x_ptr[i][1], y_ptr[i][3]);
453 }
454 accum_to_q120b(re[0], s[0], precomp);
455 accum_to_q120b(re[1], s[1], precomp);
456 accum_to_q120b(re[2], s[2], precomp);
457 accum_to_q120b(re[3], s[3], precomp);
458}
459
460EXPORT void q120x2_extract_1blk_from_q120b_ref(uint64_t nn, uint64_t blk,
461 q120x2b* const dst, // 8 doubles
462 const q120b* const src // a q120b vector
463) {
464 const uint64_t* in = (uint64_t*)src;
465 uint64_t* out = (uint64_t*)dst;
466 for (uint64_t i = 0; i < 8; ++i) {
467 out[i] = in[8 * blk + i];
468 }
469}
470
471// function on layout c is the exact same as on layout b
472#ifdef __APPLE__
473#pragma weak q120x2_extract_1blk_from_q120c_ref = q120x2_extract_1blk_from_q120b_ref
474#else
475EXPORT void q120x2_extract_1blk_from_q120c_ref(uint64_t nn, uint64_t blk,
476 q120x2c* const dst, // 8 doubles
477 const q120c* const src // a q120c vector
478 ) __attribute__((alias("q120x2_extract_1blk_from_q120b_ref")));
479#endif
480
481EXPORT void q120x2_extract_1blk_from_contiguous_q120b_ref(
482 uint64_t nn, uint64_t nrows, uint64_t blk,
483 q120x2b* const dst, // nrows * 2 q120
484 const q120b* const src // a contiguous array of nrows q120b vectors
485) {
486 const uint64_t* in = (uint64_t*)src;
487 uint64_t* out = (uint64_t*)dst;
488 for (uint64_t row = 0; row < nrows; ++row) {
489 for (uint64_t i = 0; i < 8; ++i) {
490 out[i] = in[8 * blk + i];
491 }
492 in += 4 * nn;
493 out += 8;
494 }
495}
496
497EXPORT void q120x2b_save_1blk_to_q120b_ref(uint64_t nn, uint64_t blk,
498 q120b* dest, // 1 reim vector of length m
499 const q120x2b* src // 8 doubles
500) {
501 const uint64_t* in = (uint64_t*)src;
502 uint64_t* out = (uint64_t*)dest;
503 for (uint64_t i = 0; i < 8; ++i) {
504 out[8 * blk + i] = in[i];
505 }
506}