1#include <memory.h>
 2
 3#include "zn_arithmetic_private.h"
 4
 5EXPORT TNDBL_APPROXDECOMP_GADGET* new_tndbl_approxdecomp_gadget(const MOD_Z* module,  //
 6                                                                uint64_t k, uint64_t ell) {
 7  if (k * ell > 50) {
 8    return spqlios_error("approx decomposition requested is too precise for doubles");
 9  }
10  if (k < 1) {
11    return spqlios_error("approx decomposition supports k>=1");
12  }
13  TNDBL_APPROXDECOMP_GADGET* res = malloc(sizeof(TNDBL_APPROXDECOMP_GADGET));
14  memset(res, 0, sizeof(TNDBL_APPROXDECOMP_GADGET));
15  res->k = k;
16  res->ell = ell;
17  double add_cst = INT64_C(3) << (51 - k * ell);
18  for (uint64_t i = 0; i < ell; ++i) {
19    add_cst += pow(2., -(double)(i * k + 1));
20  }
21  res->add_cst = add_cst;
22  res->and_mask = (UINT64_C(1) << k) - 1;
23  res->sub_cst = UINT64_C(1) << (k - 1);
24  for (uint64_t i = 0; i < ell; ++i) res->rshifts[i] = (ell - 1 - i) * k;
25  return res;
26}
27EXPORT void delete_tndbl_approxdecomp_gadget(TNDBL_APPROXDECOMP_GADGET* ptr) { free(ptr); }
28
29EXPORT int default_init_tndbl_approxdecomp_gadget(const MOD_Z* module,             //
30                                                  TNDBL_APPROXDECOMP_GADGET* res,  //
31                                                  uint64_t k, uint64_t ell) {
32  return 0;
33}
34
35typedef union {
36  double dv;
37  uint64_t uv;
38} du_t;
39
40#define IMPL_ixx_approxdecomp_from_tndbl_ref(ITYPE)        \
41  if (res_size != a_size * gadget->ell) NOT_IMPLEMENTED(); \
42  const uint64_t ell = gadget->ell;                        \
43  const double add_cst = gadget->add_cst;                  \
44  const uint8_t* const rshifts = gadget->rshifts;          \
45  const ITYPE and_mask = gadget->and_mask;                 \
46  const ITYPE sub_cst = gadget->sub_cst;                   \
47  ITYPE* rr = res;                                         \
48  const double* aa = a;                                    \
49  const double* aaend = a + a_size;                        \
50  while (aa < aaend) {                                     \
51    du_t t = {.dv = *aa + add_cst};                        \
52    for (uint64_t i = 0; i < ell; ++i) {                   \
53      ITYPE v = (ITYPE)(t.uv >> rshifts[i]);               \
54      *rr = (v & and_mask) - sub_cst;                      \
55      ++rr;                                                \
56    }                                                      \
57    ++aa;                                                  \
58  }
59
60/** @brief sets res = gadget_decompose(a) (int8_t* output) */
61EXPORT void default_i8_approxdecomp_from_tndbl_ref(const MOD_Z* module,                      // N
62                                                   const TNDBL_APPROXDECOMP_GADGET* gadget,  // gadget
63                                                   int8_t* res, uint64_t res_size,  // res (in general, size ell.a_size)
64                                                   const double* a, uint64_t a_size  //
65){IMPL_ixx_approxdecomp_from_tndbl_ref(int8_t)}
66
67/** @brief sets res = gadget_decompose(a) (int16_t* output) */
68EXPORT void default_i16_approxdecomp_from_tndbl_ref(const MOD_Z* module,                      // N
69                                                    const TNDBL_APPROXDECOMP_GADGET* gadget,  // gadget
70                                                    int16_t* res, uint64_t res_size,          // res
71                                                    const double* a, uint64_t a_size          // a
72){IMPL_ixx_approxdecomp_from_tndbl_ref(int16_t)}
73
74/** @brief sets res = gadget_decompose(a) (int32_t* output) */
75EXPORT void default_i32_approxdecomp_from_tndbl_ref(const MOD_Z* module,                      // N
76                                                    const TNDBL_APPROXDECOMP_GADGET* gadget,  // gadget
77                                                    int32_t* res, uint64_t res_size,          // res
78                                                    const double* a, uint64_t a_size          // a
79) {
80  IMPL_ixx_approxdecomp_from_tndbl_ref(int32_t)
81}