1#include <errno.h>
 2#include <string.h>
 3
 4#include "../commons_private.h"
 5#include "reim_fft_internal.h"
 6#include "reim_fft_private.h"
 7
 8EXPORT void reim_to_tnx_basic_ref(const REIM_TO_TNX_PRECOMP* tables, double* r, const double* x) {
 9  const uint64_t n = tables->m << 1;
10  double divisor = tables->divisor;
11  for (uint64_t i = 0; i < n; ++i) {
12    double ri = x[i] / divisor;
13    r[i] = ri - rint(ri);
14  }
15}
16
17typedef union {
18  double d;
19  uint64_t u;
20} dblui64_t;
21
22EXPORT void reim_to_tnx_ref(const REIM_TO_TNX_PRECOMP* tables, double* r, const double* x) {
23  const uint64_t n = tables->m << 1;
24  double add_cst = tables->add_cst;
25  uint64_t mask_and = tables->mask_and;
26  uint64_t mask_or = tables->mask_or;
27  double sub_cst = tables->sub_cst;
28  dblui64_t cur;
29  for (uint64_t i = 0; i < n; ++i) {
30    cur.d = x[i] + add_cst;
31    cur.u &= mask_and;
32    cur.u |= mask_or;
33    r[i] = cur.d - sub_cst;
34  }
35}
36
37void* init_reim_to_tnx_precomp(REIM_TO_TNX_PRECOMP* const res, uint32_t m, double divisor, uint32_t log2overhead) {
38  if (m & (m - 1)) return spqlios_error("m must be a power of 2");
39  if (is_not_pow2_double(&divisor)) return spqlios_error("divisor must be a power of 2");
40  if (log2overhead > 52) return spqlios_error("log2overhead is too large");
41  res->m = m;
42  res->divisor = divisor;
43  res->log2overhead = log2overhead;
44  //                 52 + 11 + 1
45  // ......1.......01(1)|expo|sign
46  // .......=========(1)|expo|sign  msbbits = log2ovh + 2 + 11 + 1
47  uint64_t nbits = 50 - log2overhead;
48  dblui64_t ovh_cst;
49  ovh_cst.d = 0.5 + (6 << log2overhead);
50  res->add_cst = ovh_cst.d * divisor;
51  res->mask_and = ((UINT64_C(1) << nbits) - 1);
52  res->mask_or = ovh_cst.u & ((UINT64_C(-1)) << nbits);
53  res->sub_cst = ovh_cst.d;
54  // TODO: check selection logic
55  if (CPU_SUPPORTS("avx2")) {
56    if (m >= 8) {
57      res->function = reim_to_tnx_avx;
58    } else {
59      res->function = reim_to_tnx_ref;
60    }
61  } else {
62    res->function = reim_to_tnx_ref;
63  }
64  return res;
65}
66
67EXPORT REIM_TO_TNX_PRECOMP* new_reim_to_tnx_precomp(uint32_t m, double divisor, uint32_t log2overhead) {
68  REIM_TO_TNX_PRECOMP* res = malloc(sizeof(*res));
69  if (!res) return spqlios_error(strerror(errno));
70  return spqlios_keep_or_free(res, init_reim_to_tnx_precomp(res, m, divisor, log2overhead));
71}
72
73EXPORT void reim_to_tnx(const REIM_TO_TNX_PRECOMP* tables, double* r, const double* x) {
74  tables->function(tables, r, x);
75}