1#include <memory.h>
  2
  3#include "../commons_private.h"
  4#include "cplx_fft.h"
  5#include "cplx_fft_internal.h"
  6#include "cplx_fft_private.h"
  7
  8/** @brief (a,b) <- (a+b,omegabar.(a-b)) */
  9void invctwiddle(CPLX a, CPLX b, const CPLX ombar) {
 10  double diffre = a[0] - b[0];
 11  double diffim = a[1] - b[1];
 12  a[0] = a[0] + b[0];
 13  a[1] = a[1] + b[1];
 14  b[0] = diffre * ombar[0] - diffim * ombar[1];
 15  b[1] = diffre * ombar[1] + diffim * ombar[0];
 16}
 17
 18/**  @brief (a,b) <- (a+b,-i.omegabar(a-b)) */
 19void invcitwiddle(CPLX a, CPLX b, const CPLX ombar) {
 20  double diffre = a[0] - b[0];
 21  double diffim = a[1] - b[1];
 22  a[0] = a[0] + b[0];
 23  a[1] = a[1] + b[1];
 24  //-i(x+iy)=-ix+y
 25  b[0] = diffre * ombar[1] + diffim * ombar[0];
 26  b[1] = -diffre * ombar[0] + diffim * ombar[1];
 27}
 28
 29/** @brief exp(-i.2pi.x) */
 30void cplx_set_e2pimx(CPLX res, double x) {
 31  res[0] = m_accurate_cos(2 * M_PI * x);
 32  res[1] = -m_accurate_sin(2 * M_PI * x);
 33}
 34/**
 35 * @brief this computes the sequence: 0,1/2,1/4,3/4,1/8,5/8,3/8,7/8,...
 36 * essentially: the bits of (i+1) in lsb order on the basis (1/2^k) mod 1*/
 37double fracrevbits(uint32_t i);
 38/** @brief fft modulo X^m-exp(i.2pi.entry+pwr) -- reference code */
 39void cplx_ifft_naive(const uint32_t m, const double entry_pwr, CPLX* data) {
 40  if (m == 1) return;
 41  const double pom = entry_pwr / 2.;
 42  const uint32_t h = m / 2;
 43  CPLX cpom;
 44  cplx_set_e2pimx(cpom, pom);
 45  // do the recursive calls
 46  cplx_ifft_naive(h, pom, data);
 47  cplx_ifft_naive(h, pom + 0.5, data + h);
 48  // apply the inverse twiddle factors
 49  for (uint64_t i = 0; i < h; ++i) {
 50    invctwiddle(data[i], data[i + h], cpom);
 51  }
 52}
 53
 54void cplx_ifft16_precomp(const double entry_pwr, CPLX** omg) {
 55  static const double j_pow = 1. / 8.;
 56  static const double k_pow = 1. / 16.;
 57  const double pom = entry_pwr / 2.;
 58  const double pom_2 = entry_pwr / 4.;
 59  const double pom_4 = entry_pwr / 8.;
 60  const double pom_8 = entry_pwr / 16.;
 61  cplx_set_e2pimx((*omg)[0], pom_8);
 62  cplx_set_e2pimx((*omg)[1], pom_8 + j_pow);
 63  cplx_set_e2pimx((*omg)[2], pom_8 + k_pow);
 64  cplx_set_e2pimx((*omg)[3], pom_8 + j_pow + k_pow);
 65  cplx_set_e2pimx((*omg)[4], pom_4);
 66  cplx_set_e2pimx((*omg)[5], pom_4 + j_pow);
 67  cplx_set_e2pimx((*omg)[6], pom_2);
 68  cplx_set_e2pimx((*omg)[7], pom);
 69  *omg += 8;
 70}
 71
 72/**
 73 * @brief iFFT modulo X^16-omega^2 (in registers)
 74 * @param data contains 16 complexes
 75 * @param omegabar 8 complexes in this order:
 76 *     gammabar,jb.gammabar,kb.gammabar,kbjb.gammabar,
 77 *     betabar,jb.betabar,alphabar,omegabar
 78 *     alpha = sqrt(omega), beta = sqrt(alpha), gamma = sqrt(beta)
 79 *     jb = sqrt(ib), kb=sqrt(jb)
 80 */
 81void cplx_ifft16_ref(void* data, const void* omegabar) {
 82  CPLX* d = data;
 83  const CPLX* om = omegabar;
 84  // fourth pass inverse
 85  invctwiddle(d[0], d[1], om[0]);
 86  invcitwiddle(d[2], d[3], om[0]);
 87  invctwiddle(d[4], d[5], om[1]);
 88  invcitwiddle(d[6], d[7], om[1]);
 89  invctwiddle(d[8], d[9], om[2]);
 90  invcitwiddle(d[10], d[11], om[2]);
 91  invctwiddle(d[12], d[13], om[3]);
 92  invcitwiddle(d[14], d[15], om[3]);
 93  // third pass inverse
 94  invctwiddle(d[0], d[2], om[4]);
 95  invctwiddle(d[1], d[3], om[4]);
 96  invcitwiddle(d[4], d[6], om[4]);
 97  invcitwiddle(d[5], d[7], om[4]);
 98  invctwiddle(d[8], d[10], om[5]);
 99  invctwiddle(d[9], d[11], om[5]);
100  invcitwiddle(d[12], d[14], om[5]);
101  invcitwiddle(d[13], d[15], om[5]);
102  // second pass inverse
103  invctwiddle(d[0], d[4], om[6]);
104  invctwiddle(d[1], d[5], om[6]);
105  invctwiddle(d[2], d[6], om[6]);
106  invctwiddle(d[3], d[7], om[6]);
107  invcitwiddle(d[8], d[12], om[6]);
108  invcitwiddle(d[9], d[13], om[6]);
109  invcitwiddle(d[10], d[14], om[6]);
110  invcitwiddle(d[11], d[15], om[6]);
111  // first pass
112  for (uint64_t i = 0; i < 8; ++i) {
113    invctwiddle(d[0 + i], d[8 + i], om[7]);
114  }
115}
116
117void cplx_ifft_ref_bfs_2(CPLX* dat, const CPLX** omg, uint32_t m) {
118  CPLX* const dend = dat + m;
119  for (CPLX* d = dat; d < dend; d += 2) {
120    split_fft_last_ref(d, (*omg)[0]);
121    *omg += 1;
122  }
123#if 0
124  printf("after first: ");
125  for (uint64_t ii=0; ii<nn/2; ++ii) {
126    printf("%.6lf %.6lf ",data[ii][0],data[ii][1]);
127  }
128  printf("\n");
129#endif
130  int32_t Ms2 = m / 2;
131  for (int32_t h = 2; h <= Ms2; h <<= 1) {
132    for (CPLX* d = dat; d < dend; d += 2 * h) {
133      if (memcmp((*omg)[0], (*omg)[1], 8) != 0) abort();
134      cplx_split_fft_ref(h, d, **omg);
135      *omg += 2;
136    }
137#if 0
138    printf("after split %d: ", h);
139    for (uint64_t ii=0; ii<nn/2; ++ii) {
140      printf("%.6lf %.6lf ",data[ii][0],data[ii][1]);
141    }
142    printf("\n");
143#endif
144  }
145}
146
147void cplx_ifft_ref_bfs_16(CPLX* dat, const CPLX** omg, uint32_t m) {
148  const uint64_t log2m = log2(m);
149  CPLX* const dend = dat + m;
150  // h=1,2,4,8 use the 16-dim macroblock
151  for (CPLX* d = dat; d < dend; d += 16) {
152    cplx_ifft16_ref(d, *omg);
153    *omg += 8;
154  }
155#if 0
156  printf("after first: ");
157  for (uint64_t ii=0; ii<nn/2; ++ii) {
158    printf("%.6lf %.6lf ",data[ii][0],data[ii][1]);
159  }
160  printf("\n");
161#endif
162  int32_t h = 16;
163  if (log2m % 2 != 0) {
164    // if parity needs it, uses one regular twiddle
165    for (CPLX* d = dat; d < dend; d += 2 * h) {
166      cplx_split_fft_ref(h, d, **omg);
167      *omg += 1;
168    }
169    h = 32;
170  }
171  // h=16,...,2*floor(Ms2/2) use the bitwiddle
172  for (; h < m; h <<= 2) {
173    for (CPLX* d = dat; d < dend; d += 4 * h) {
174      cplx_bisplit_fft_ref(h, d, *omg);
175      *omg += 2;
176    }
177#if 0
178    printf("after split %d: ", h);
179    for (uint64_t ii=0; ii<nn/2; ++ii) {
180      printf("%.6lf %.6lf ",data[ii][0],data[ii][1]);
181    }
182    printf("\n");
183#endif
184  }
185}
186
187void fill_cplx_ifft_omegas_bfs_16(const double entry_pwr, CPLX** omg, uint32_t m) {
188  const uint64_t log2m = log2(m);
189  double pwr = entry_pwr * 16. / m;
190  {
191    // h=8
192    for (uint32_t i = 0; i < m / 16; i++) {
193      cplx_ifft16_precomp(pwr + fracrevbits(i), omg);
194    }
195  }
196  int32_t h = 16;
197  if (log2m % 2 != 0) {
198    // if parity needs it, uses one regular twiddle
199    for (uint32_t i = 0; i < m / (2 * h); i++) {
200      cplx_set_e2pimx(omg[0][0], pwr + fracrevbits(i) / 2.);
201      *omg += 1;
202    }
203    pwr *= 2.;
204    h = 32;
205  }
206  for (; h < m; h <<= 2) {
207    for (uint32_t i = 0; i < m / (2 * h); i += 2) {
208      cplx_set_e2pimx(omg[0][0], pwr + fracrevbits(i) / 2.);
209      cplx_set_e2pimx(omg[0][1], 2. * pwr + fracrevbits(i));
210      *omg += 2;
211    }
212    pwr *= 4.;
213  }
214}
215
216void fill_cplx_ifft_omegas_bfs_2(const double entry_pwr, CPLX** omg, uint32_t m) {
217  double pom = entry_pwr / m;
218  {
219    // h=1
220    for (uint32_t i = 0; i < m / 2; i++) {
221      cplx_set_e2pimx((*omg)[0], pom + fracrevbits(i) / 2.);
222      *omg += 1;  // optim function reads by 1
223    }
224  }
225  for (int32_t h = 2; h <= m / 2; h <<= 1) {
226    pom *= 2;
227    for (uint32_t i = 0; i < m / (2 * h); i++) {
228      cplx_set_e2pimx(omg[0][0], pom + fracrevbits(i) / 2.);
229      cplx_set(omg[0][1], omg[0][0]);
230      *omg += 2;
231    }
232  }
233}
234
235void fill_cplx_ifft_omegas_rec_16(const double entry_pwr, CPLX** omg, uint32_t m) {
236  if (m == 1) return;
237  if (m <= 8) return fill_cplx_ifft_omegas_bfs_2(entry_pwr, omg, m);
238  if (m <= 2048) return fill_cplx_ifft_omegas_bfs_16(entry_pwr, omg, m);
239  double pom = entry_pwr / 2.;
240  uint32_t h = m / 2;
241  fill_cplx_ifft_omegas_rec_16(pom, omg, h);
242  fill_cplx_ifft_omegas_rec_16(pom + 0.5, omg, h);
243  cplx_set_e2pimx((*omg)[0], pom);
244  cplx_set((*omg)[1], (*omg)[0]);
245  *omg += 2;
246}
247
248void cplx_ifft_ref_rec_16(CPLX* dat, const CPLX** omg, uint32_t m) {
249  if (m == 1) return;
250  if (m <= 8) return cplx_ifft_ref_bfs_2(dat, omg, m);
251  if (m <= 2048) return cplx_ifft_ref_bfs_16(dat, omg, m);
252  const uint32_t h = m / 2;
253  cplx_ifft_ref_rec_16(dat, omg, h);
254  cplx_ifft_ref_rec_16(dat + h, omg, h);
255  if (memcmp((*omg)[0], (*omg)[1], 8) != 0) abort();
256  cplx_split_fft_ref(h, dat, **omg);
257  *omg += 2;
258}
259
260EXPORT void cplx_ifft_ref(const CPLX_IFFT_PRECOMP* precomp, void* d) {
261  CPLX* data = (CPLX*)d;
262  const int32_t m = precomp->m;
263  const CPLX* omg = (CPLX*)precomp->powomegas;
264  if (m == 1) return;
265  if (m <= 8) return cplx_ifft_ref_bfs_2(data, &omg, m);
266  if (m <= 2048) return cplx_ifft_ref_bfs_16(data, &omg, m);
267  cplx_ifft_ref_rec_16(data, &omg, m);
268}
269
270EXPORT CPLX_IFFT_PRECOMP* new_cplx_ifft_precomp(uint32_t m, uint32_t num_buffers) {
271  const uint64_t OMG_SPACE = ceilto64b(2 * m * sizeof(CPLX));
272  const uint64_t BUF_SIZE = ceilto64b(m * sizeof(CPLX));
273  void* reps = malloc(sizeof(CPLX_IFFT_PRECOMP) + 63  // padding
274                      + OMG_SPACE                     // tables
275                      + num_buffers * BUF_SIZE        // buffers
276  );
277  uint64_t aligned_addr = ceilto64b((uint64_t)reps + sizeof(CPLX_IFFT_PRECOMP));
278  CPLX_IFFT_PRECOMP* r = (CPLX_IFFT_PRECOMP*)reps;
279  r->m = m;
280  r->buf_size = BUF_SIZE;
281  r->powomegas = (double*)aligned_addr;
282  r->aligned_buffers = (void*)(aligned_addr + OMG_SPACE);
283  // fill in powomegas
284  CPLX* omg = (CPLX*)r->powomegas;
285  fill_cplx_ifft_omegas_rec_16(0.25, &omg, m);
286  if (((uint64_t)omg) - aligned_addr > OMG_SPACE) abort();
287  {
288    if (m <= 4) {
289      // currently, we do not have any acceletated
290      // implementation for m<=4
291      r->function = cplx_ifft_ref;
292    } else if (CPU_SUPPORTS("fma")) {
293      r->function = cplx_ifft_avx2_fma;
294    } else {
295      r->function = cplx_ifft_ref;
296    }
297  }
298  return reps;
299}
300
301EXPORT void* cplx_ifft_precomp_get_buffer(const CPLX_IFFT_PRECOMP* itables, uint32_t buffer_index) {
302  return (uint8_t*)itables->aligned_buffers + buffer_index * itables->buf_size;
303}
304
305EXPORT void cplx_ifft(const CPLX_IFFT_PRECOMP* itables, void* data) { itables->function(itables, data); }
306
307EXPORT void cplx_ifft_simple(uint32_t m, void* data) {
308  static CPLX_IFFT_PRECOMP* p[31] = {0};
309  CPLX_IFFT_PRECOMP** f = p + log2m(m);
310  if (!*f) *f = new_cplx_ifft_precomp(m, 0);
311  (*f)->function(*f, data);
312}