1#include <memory.h>
 2
 3#include "vec_rnx_arithmetic_private.h"
 4
 5typedef union di {
 6  double dv;
 7  uint64_t uv;
 8} di_t;
 9
10/** @brief new gadget: delete with delete_tnxdbl_approxdecomp_gadget */
11EXPORT TNXDBL_APPROXDECOMP_GADGET* new_tnxdbl_approxdecomp_gadget(  //
12    const MOD_RNX* module,                                          // N
13    uint64_t k, uint64_t ell                                        // base 2^K and size
14) {
15  if (k * ell > 50) return spqlios_error("gadget requires a too large fp precision");
16  TNXDBL_APPROXDECOMP_GADGET* res = spqlios_alloc(sizeof(TNXDBL_APPROXDECOMP_GADGET));
17  res->k = k;
18  res->ell = ell;
19  // double add_cst;       // double(3.2^(51-ell.K) + 1/2.(sum 2^(-iK)) for i=[0,ell[)
20  union di add_cst;
21  add_cst.dv = UINT64_C(3) << (51 - ell * k);
22  for (uint64_t i = 0; i < ell; ++i) {
23    add_cst.uv |= UINT64_C(1) << ((i + 1) * k - 1);
24  }
25  res->add_cst = add_cst.dv;
26  // uint64_t and_mask;     // uint64(2^(K)-1)
27  res->and_mask = (UINT64_C(1) << k) - 1;
28  // uint64_t or_mask;      // double(2^52)
29  union di or_mask;
30  or_mask.dv = (UINT64_C(1) << 52);
31  res->or_mask = or_mask.uv;
32  // double sub_cst;       // double(2^52 + 2^(K-1))
33  res->sub_cst = ((UINT64_C(1) << 52) + (UINT64_C(1) << (k - 1)));
34  return res;
35}
36
37EXPORT void delete_tnxdbl_approxdecomp_gadget(TNXDBL_APPROXDECOMP_GADGET* gadget) { spqlios_free(gadget); }
38
39/** @brief sets res = gadget_decompose(a) */
40EXPORT void rnx_approxdecomp_from_tnxdbl_ref(         //
41    const MOD_RNX* module,                            // N
42    const TNXDBL_APPROXDECOMP_GADGET* gadget,         // output base 2^K
43    double* res, uint64_t res_size, uint64_t res_sl,  // res
44    const double* a                                   // a
45) {
46  const uint64_t nn = module->n;
47  const uint64_t k = gadget->k;
48  const uint64_t ell = gadget->ell;
49  const double add_cst = gadget->add_cst;
50  const uint64_t and_mask = gadget->and_mask;
51  const uint64_t or_mask = gadget->or_mask;
52  const double sub_cst = gadget->sub_cst;
53  const uint64_t msize = res_size <= ell ? res_size : ell;
54  const uint64_t first_rsh = (ell - msize) * k;
55  // gadget decompose column by column
56  if (msize > 0) {
57    double* const last_r = res + (msize - 1) * res_sl;
58    for (uint64_t j = 0; j < nn; ++j) {
59      double* rr = last_r + j;
60      di_t t = {.dv = a[j] + add_cst};
61      if (msize < ell) t.uv >>= first_rsh;
62      do {
63        di_t u;
64        u.uv = (t.uv & and_mask) | or_mask;
65        *rr = u.dv - sub_cst;
66        t.uv >>= k;
67        rr -= res_sl;
68      } while (rr >= res);
69    }
70  }
71  // zero-out the last slices (if any)
72  for (uint64_t i = msize; i < res_size; ++i) {
73    memset(res + i * res_sl, 0, nn * sizeof(double));
74  }
75}