1#include <errno.h>
  2#include <string.h>
  3
  4#include "../commons_private.h"
  5#include "cplx_fft_internal.h"
  6#include "cplx_fft_private.h"
  7
  8EXPORT void cplx_from_znx32_ref(const CPLX_FROM_ZNX32_PRECOMP* precomp, void* r, const int32_t* x) {
  9  const uint32_t m = precomp->m;
 10  const int32_t* inre = x;
 11  const int32_t* inim = x + m;
 12  CPLX* out = r;
 13  for (uint32_t i = 0; i < m; ++i) {
 14    out[i][0] = (double)inre[i];
 15    out[i][1] = (double)inim[i];
 16  }
 17}
 18
 19EXPORT void cplx_from_tnx32_ref(const CPLX_FROM_TNX32_PRECOMP* precomp, void* r, const int32_t* x) {
 20  static const double _2p32 = 1. / (INT64_C(1) << 32);
 21  const uint32_t m = precomp->m;
 22  const int32_t* inre = x;
 23  const int32_t* inim = x + m;
 24  CPLX* out = r;
 25  for (uint32_t i = 0; i < m; ++i) {
 26    out[i][0] = ((double)inre[i]) * _2p32;
 27    out[i][1] = ((double)inim[i]) * _2p32;
 28  }
 29}
 30
 31EXPORT void cplx_to_tnx32_ref(const CPLX_TO_TNX32_PRECOMP* precomp, int32_t* r, const void* x) {
 32  static const double _2p32 = (INT64_C(1) << 32);
 33  const uint32_t m = precomp->m;
 34  double factor = _2p32 / precomp->divisor;
 35  int32_t* outre = r;
 36  int32_t* outim = r + m;
 37  const CPLX* in = x;
 38  // Note: this formula will only work if abs(in) < 2^32
 39  for (uint32_t i = 0; i < m; ++i) {
 40    outre[i] = (int32_t)(int64_t)(rint(in[i][0] * factor));
 41    outim[i] = (int32_t)(int64_t)(rint(in[i][1] * factor));
 42  }
 43}
 44
 45void* init_cplx_from_znx32_precomp(CPLX_FROM_ZNX32_PRECOMP* res, uint32_t m) {
 46  res->m = m;
 47  if (CPU_SUPPORTS("avx2")) {
 48    if (m >= 8) {
 49      res->function = cplx_from_znx32_avx2_fma;
 50    } else {
 51      res->function = cplx_from_znx32_ref;
 52    }
 53  } else {
 54    res->function = cplx_from_znx32_ref;
 55  }
 56  return res;
 57}
 58
 59CPLX_FROM_ZNX32_PRECOMP* new_cplx_from_znx32_precomp(uint32_t m) {
 60  CPLX_FROM_ZNX32_PRECOMP* res = malloc(sizeof(CPLX_FROM_ZNX32_PRECOMP));
 61  if (!res) return spqlios_error(strerror(errno));
 62  return spqlios_keep_or_free(res, init_cplx_from_znx32_precomp(res, m));
 63}
 64
 65void* init_cplx_from_tnx32_precomp(CPLX_FROM_TNX32_PRECOMP* res, uint32_t m) {
 66  res->m = m;
 67  if (CPU_SUPPORTS("avx2")) {
 68    if (m >= 8) {
 69      res->function = cplx_from_tnx32_avx2_fma;
 70    } else {
 71      res->function = cplx_from_tnx32_ref;
 72    }
 73  } else {
 74    res->function = cplx_from_tnx32_ref;
 75  }
 76  return res;
 77}
 78
 79CPLX_FROM_TNX32_PRECOMP* new_cplx_from_tnx32_precomp(uint32_t m) {
 80  CPLX_FROM_TNX32_PRECOMP* res = malloc(sizeof(CPLX_FROM_TNX32_PRECOMP));
 81  if (!res) return spqlios_error(strerror(errno));
 82  return spqlios_keep_or_free(res, init_cplx_from_tnx32_precomp(res, m));
 83}
 84
 85void* init_cplx_to_tnx32_precomp(CPLX_TO_TNX32_PRECOMP* res, uint32_t m, double divisor, uint32_t log2overhead) {
 86  if (is_not_pow2_double(&divisor)) return spqlios_error("divisor must be a power of 2");
 87  if (m & (m - 1)) return spqlios_error("m must be a power of 2");
 88  if (log2overhead > 52) return spqlios_error("log2overhead is too large");
 89  res->m = m;
 90  res->divisor = divisor;
 91  if (CPU_SUPPORTS("avx2")) {
 92    if (log2overhead <= 18) {
 93      if (m >= 8) {
 94        res->function = cplx_to_tnx32_avx2_fma;
 95      } else {
 96        res->function = cplx_to_tnx32_ref;
 97      }
 98    } else {
 99      res->function = cplx_to_tnx32_ref;
100    }
101  } else {
102    res->function = cplx_to_tnx32_ref;
103  }
104  return res;
105}
106
107EXPORT CPLX_TO_TNX32_PRECOMP* new_cplx_to_tnx32_precomp(uint32_t m, double divisor, uint32_t log2overhead) {
108  CPLX_TO_TNX32_PRECOMP* res = malloc(sizeof(CPLX_TO_TNX32_PRECOMP));
109  if (!res) return spqlios_error(strerror(errno));
110  return spqlios_keep_or_free(res, init_cplx_to_tnx32_precomp(res, m, divisor, log2overhead));
111}
112
113/**
114 * @brief Simpler API for the znx32 to cplx conversion.
115 * For each dimension, the precomputed tables for this dimension are generated automatically the first time.
116 * It is advised to do one dry-run call per desired dimension before using in a multithread environment */
117EXPORT void cplx_from_znx32_simple(uint32_t m, void* r, const int32_t* x) {
118  // not checking for log2bound which is not relevant here
119  static CPLX_FROM_ZNX32_PRECOMP precomp[32];
120  CPLX_FROM_ZNX32_PRECOMP* p = precomp + log2m(m);
121  if (!p->function) {
122    if (!init_cplx_from_znx32_precomp(p, m)) abort();
123  }
124  p->function(p, r, x);
125}
126
127/**
128 * @brief Simpler API for the tnx32 to cplx conversion.
129 * For each dimension, the precomputed tables for this dimension are generated automatically the first time.
130 * It is advised to do one dry-run call per desired dimension before using in a multithread environment */
131EXPORT void cplx_from_tnx32_simple(uint32_t m, void* r, const int32_t* x) {
132  static CPLX_FROM_TNX32_PRECOMP precomp[32];
133  CPLX_FROM_TNX32_PRECOMP* p = precomp + log2m(m);
134  if (!p->function) {
135    if (!init_cplx_from_tnx32_precomp(p, m)) abort();
136  }
137  p->function(p, r, x);
138}
139/**
140 * @brief Simpler API for the cplx to tnx32 conversion.
141 * For each dimension, the precomputed tables for this dimension are generated automatically the first time.
142 * It is advised to do one dry-run call per desired dimension before using in a multithread environment */
143EXPORT void cplx_to_tnx32_simple(uint32_t m, double divisor, uint32_t log2overhead, int32_t* r, const void* x) {
144  struct LAST_CPLX_TO_TNX32_PRECOMP {
145    CPLX_TO_TNX32_PRECOMP p;
146    double last_divisor;
147    double last_log2over;
148  };
149  static __thread struct LAST_CPLX_TO_TNX32_PRECOMP precomp[32];
150  struct LAST_CPLX_TO_TNX32_PRECOMP* p = precomp + log2m(m);
151  if (!p->p.function || divisor != p->last_divisor || log2overhead != p->last_log2over) {
152    memset(p, 0, sizeof(*p));
153    if (!init_cplx_to_tnx32_precomp(&p->p, m, divisor, log2overhead)) abort();
154    p->last_divisor = divisor;
155    p->last_log2over = log2overhead;
156  }
157  p->p.function(&p->p, r, x);
158}