1#ifndef SPQLIOS_REIM_FFT_CORE_TEMPLATE_H
  2#define SPQLIOS_REIM_FFT_CORE_TEMPLATE_H
  3
  4// this file contains the main template for the fft strategy.
  5// it is meant to be included once for each specialization (ref, avx, neon)
  6// all the leaf functions it uses shall be defined in the following macros,
  7// before including this header.
  8
  9#if !defined(reim_fft16_f) || !defined(reim_fft16_pom_offset) || !defined(fill_reim_fft16_omegas_f)
 10#error "missing reim16 definitions"
 11#endif
 12#if !defined(reim_fft8_f) || !defined(reim_fft8_pom_offset) || !defined(fill_reim_fft8_omegas_f)
 13#error "missing reim8 definitions"
 14#endif
 15#if !defined(reim_fft4_f) || !defined(reim_fft4_pom_offset) || !defined(fill_reim_fft4_omegas_f)
 16#error "missing reim4 definitions"
 17#endif
 18#if !defined(reim_fft2_f) || !defined(reim_fft2_pom_offset) || !defined(fill_reim_fft2_omegas_f)
 19#error "missing reim2 definitions"
 20#endif
 21#if !defined(reim_twiddle_fft_f) || !defined(reim_twiddle_fft_pom_offset) || !defined(fill_reim_twiddle_fft_omegas_f)
 22#error "missing twiddle definitions"
 23#endif
 24#if !defined(reim_bitwiddle_fft_f) || !defined(reim_bitwiddle_fft_pom_offset) || \
 25    !defined(fill_reim_bitwiddle_fft_omegas_f)
 26#error "missing bitwiddle definitions"
 27#endif
 28#if !defined(reim_fft_bfs_16_f) || !defined(fill_reim_fft_bfs_16_omegas_f)
 29#error "missing bfs_16 definitions"
 30#endif
 31#if !defined(reim_fft_rec_16_f) || !defined(fill_reim_fft_rec_16_omegas_f)
 32#error "missing rec_16 definitions"
 33#endif
 34#if !defined(reim_fft_f) || !defined(fill_reim_fft_omegas_f)
 35#error "missing main definitions"
 36#endif
 37
 38void reim_fft_bfs_16_f(uint64_t m, double* re, double* im, double** omg) {
 39  uint64_t log2m = log2(m);
 40  uint64_t mm = m;
 41  if (log2m & 1) {
 42    uint64_t h = mm >> 1;
 43    // do the first twiddle iteration normally
 44    reim_twiddle_fft_f(h, re, im, *omg);
 45    *omg += reim_twiddle_fft_pom_offset;
 46    mm = h;
 47  }
 48  while (mm > 16) {
 49    uint64_t h = mm >> 2;
 50    for (uint64_t off = 0; off < m; off += mm) {
 51      reim_bitwiddle_fft_f(h, re + off, im + off, *omg);
 52      *omg += reim_bitwiddle_fft_pom_offset;
 53    }
 54    mm = h;
 55  }
 56  if (mm != 16) abort();  // bug!
 57  for (uint64_t off = 0; off < m; off += 16) {
 58    reim_fft16_f(re + off, im + off, *omg);
 59    *omg += reim_fft16_pom_offset;
 60  }
 61}
 62
 63void fill_reim_fft_bfs_16_omegas_f(uint64_t m, double entry_pwr, double** omg) {
 64  uint64_t log2m = log2(m);
 65  uint64_t mm = m;
 66  double ss = entry_pwr;
 67  if (log2m % 2 != 0) {
 68    uint64_t h = mm >> 1;
 69    double s = ss / 2.;
 70    // do the first twiddle iteration normally
 71    fill_reim_twiddle_fft_omegas_f(s, omg);
 72    mm = h;
 73    ss = s;
 74  }
 75  while (mm > 16) {
 76    uint64_t h = mm >> 2;
 77    double s = ss / 4.;
 78    for (uint64_t off = 0; off < m; off += mm) {
 79      double rs0 = s + fracrevbits(off / mm) / 4.;
 80      fill_reim_bitwiddle_fft_omegas_f(rs0, omg);
 81    }
 82    mm = h;
 83    ss = s;
 84  }
 85  if (mm != 16) abort();  // bug!
 86  for (uint64_t off = 0; off < m; off += 16) {
 87    double s = ss + fracrevbits(off / 16);
 88    fill_reim_fft16_omegas_f(s, omg);
 89  }
 90}
 91
 92void reim_fft_rec_16_f(uint64_t m, double* re, double* im, double** omg) {
 93  if (m <= 2048) return reim_fft_bfs_16_f(m, re, im, omg);
 94  const uint32_t h = m >> 1;
 95  reim_twiddle_fft_f(h, re, im, *omg);
 96  *omg += reim_twiddle_fft_pom_offset;
 97  reim_fft_rec_16_f(h, re, im, omg);
 98  reim_fft_rec_16_f(h, re + h, im + h, omg);
 99}
100
101void fill_reim_fft_rec_16_omegas_f(uint64_t m, double entry_pwr, double** omg) {
102  if (m <= 2048) return fill_reim_fft_bfs_16_omegas_f(m, entry_pwr, omg);
103  const uint64_t h = m / 2;
104  const double s = entry_pwr / 2;
105  fill_reim_twiddle_fft_omegas_f(s, omg);
106  fill_reim_fft_rec_16_omegas_f(h, s, omg);
107  fill_reim_fft_rec_16_omegas_f(h, s + 0.5, omg);
108}
109
110void reim_fft_f(const REIM_FFT_PRECOMP* precomp, double* dat) {
111  const int32_t m = precomp->m;
112  double* omg = precomp->powomegas;
113  double* re = dat;
114  double* im = dat + m;
115  if (m <= 16) {
116    switch (m) {
117      case 1:
118        return;
119      case 2:
120        return reim_fft2_f(re, im, omg);
121      case 4:
122        return reim_fft4_f(re, im, omg);
123      case 8:
124        return reim_fft8_f(re, im, omg);
125      case 16:
126        return reim_fft16_f(re, im, omg);
127      default:
128        abort();  // m is not a power of 2
129    }
130  }
131  if (m <= 2048) return reim_fft_bfs_16_f(m, re, im, &omg);
132  return reim_fft_rec_16_f(m, re, im, &omg);
133}
134
135void fill_reim_fft_omegas_f(uint64_t m, double entry_pwr, double** omg) {
136  if (m <= 16) {
137    switch (m) {
138      case 1:
139        break;
140      case 2:
141        fill_reim_fft2_omegas_f(entry_pwr, omg);
142        break;
143      case 4:
144        fill_reim_fft4_omegas_f(entry_pwr, omg);
145        break;
146      case 8:
147        fill_reim_fft8_omegas_f(entry_pwr, omg);
148        break;
149      case 16:
150        fill_reim_fft16_omegas_f(entry_pwr, omg);
151        break;
152      default:
153        abort();  // m is not a power of 2
154    }
155  } else if (m <= 2048) {
156    fill_reim_fft_bfs_16_omegas_f(m, entry_pwr, omg);
157  } else {
158    fill_reim_fft_rec_16_omegas_f(m, entry_pwr, omg);
159  }
160}
161
162#endif  // SPQLIOS_REIM_FFT_CORE_TEMPLATE_H