poulpy/spqlios/lib/spqlios/arithmetic/vec_rnx_svp_ref.c
spqlios/lib/spqlios/arithmetic/vec_rnx_svp_ref.cRaw
1#include <string.h>
2
3#include "../coeffs/coeffs_arithmetic.h"
4#include "vec_rnx_arithmetic_private.h"
5
6EXPORT uint64_t fft64_bytes_of_rnx_svp_ppol(const MOD_RNX* module) { return module->n * sizeof(double); }
7
8EXPORT RNX_SVP_PPOL* new_rnx_svp_ppol(const MOD_RNX* module) { return spqlios_alloc(bytes_of_rnx_svp_ppol(module)); }
9
10EXPORT void delete_rnx_svp_ppol(RNX_SVP_PPOL* ppol) { spqlios_free(ppol); }
11
12/** @brief prepares a svp polynomial */
13EXPORT void fft64_rnx_svp_prepare_ref(const MOD_RNX* module, // N
14 RNX_SVP_PPOL* ppol, // output
15 const double* pol // a
16) {
17 double* const dppol = (double*)ppol;
18 rnx_divide_by_m_ref(module->n, module->m, dppol, pol);
19 reim_fft(module->precomp.fft64.p_fft, dppol);
20}
21
22EXPORT void fft64_rnx_svp_apply_ref( //
23 const MOD_RNX* module, // N
24 double* res, uint64_t res_size, uint64_t res_sl, // output
25 const RNX_SVP_PPOL* ppol, // prepared pol
26 const double* a, uint64_t a_size, uint64_t a_sl // a
27) {
28 const uint64_t nn = module->n;
29 double* const dppol = (double*)ppol;
30
31 const uint64_t auto_end_idx = res_size < a_size ? res_size : a_size;
32 for (uint64_t i = 0; i < auto_end_idx; ++i) {
33 const double* a_ptr = a + i * a_sl;
34 double* const res_ptr = res + i * res_sl;
35 // copy the polynomial to res, apply fft in place, call fftvec
36 // _mul, apply ifft in place.
37 memcpy(res_ptr, a_ptr, nn * sizeof(double));
38 reim_fft(module->precomp.fft64.p_fft, (double*)res_ptr);
39 reim_fftvec_mul(module->precomp.fft64.p_fftvec_mul, res_ptr, res_ptr, dppol);
40 reim_ifft(module->precomp.fft64.p_ifft, res_ptr);
41 }
42
43 // then extend with zeros
44 for (uint64_t i = auto_end_idx; i < res_size; ++i) {
45 memset(res + i * res_sl, 0, nn * sizeof(double));
46 }
47}