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