1#include "fft64_dft.h"
  2
  3#include <cstring>
  4
  5#include "../../spqlios/reim/reim_fft.h"
  6#include "../../spqlios/reim/reim_fft_internal.h"
  7
  8reim_fft64vec::reim_fft64vec(uint64_t n) : v(n, 0) {}
  9reim4_elem reim_fft64vec::get_blk(uint64_t blk) const {
 10  return reim_view(v.size() / 2, (double*)v.data()).get_blk(blk);
 11}
 12double* reim_fft64vec::data() { return v.data(); }
 13const double* reim_fft64vec::data() const { return v.data(); }
 14uint64_t reim_fft64vec::nn() const { return v.size(); }
 15reim_fft64vec::reim_fft64vec(uint64_t n, const double* data) : v(data, data + n) {}
 16void reim_fft64vec::save_as(double* dest) const { memcpy(dest, v.data(), nn() * sizeof(double)); }
 17reim_fft64vec reim_fft64vec::zero(uint64_t n) { return reim_fft64vec(n); }
 18void reim_fft64vec::set_blk(uint64_t blk, const reim4_elem& value) {
 19  reim_view(v.size() / 2, (double*)v.data()).set_blk(blk, value);
 20}
 21reim_fft64vec reim_fft64vec::dft_random(uint64_t n, uint64_t log2bound) {
 22  return simple_fft64(znx_i64::random_log2bound(n, log2bound));
 23}
 24reim_fft64vec reim_fft64vec::random(uint64_t n, double log2bound) {
 25  double bound = pow(2., log2bound);
 26  reim_fft64vec res(n);
 27  for (uint64_t i = 0; i < n; ++i) {
 28    res.v[i] = uniform_f64_bounds(-bound, bound);
 29  }
 30  return res;
 31}
 32
 33reim_fft64vec operator+(const reim_fft64vec& a, const reim_fft64vec& b) {
 34  uint64_t nn = a.nn();
 35  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "ring dimension mismatch");
 36  reim_fft64vec res(nn);
 37  double* rv = res.data();
 38  const double* av = a.data();
 39  const double* bv = b.data();
 40  for (uint64_t i = 0; i < nn; ++i) {
 41    rv[i] = av[i] + bv[i];
 42  }
 43  return res;
 44}
 45reim_fft64vec operator-(const reim_fft64vec& a, const reim_fft64vec& b) {
 46  uint64_t nn = a.nn();
 47  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "ring dimension mismatch");
 48  reim_fft64vec res(nn);
 49  double* rv = res.data();
 50  const double* av = a.data();
 51  const double* bv = b.data();
 52  for (uint64_t i = 0; i < nn; ++i) {
 53    rv[i] = av[i] - bv[i];
 54  }
 55  return res;
 56}
 57reim_fft64vec operator*(const reim_fft64vec& a, const reim_fft64vec& b) {
 58  uint64_t nn = a.nn();
 59  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "ring dimension mismatch");
 60  REQUIRE_DRAMATICALLY(nn >= 2, "test not defined for nn=1");
 61  uint64_t m = nn / 2;
 62  reim_fft64vec res(nn);
 63  double* rv = res.data();
 64  const double* av = a.data();
 65  const double* bv = b.data();
 66  for (uint64_t i = 0; i < m; ++i) {
 67    rv[i] = av[i] * bv[i] - av[m + i] * bv[m + i];
 68    rv[m + i] = av[i] * bv[m + i] + av[m + i] * bv[i];
 69  }
 70  return res;
 71}
 72reim_fft64vec& operator+=(reim_fft64vec& a, const reim_fft64vec& b) {
 73  uint64_t nn = a.nn();
 74  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "ring dimension mismatch");
 75  double* av = a.data();
 76  const double* bv = b.data();
 77  for (uint64_t i = 0; i < nn; ++i) {
 78    av[i] = av[i] + bv[i];
 79  }
 80  return a;
 81}
 82reim_fft64vec& operator-=(reim_fft64vec& a, const reim_fft64vec& b) {
 83  uint64_t nn = a.nn();
 84  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "ring dimension mismatch");
 85  double* av = a.data();
 86  const double* bv = b.data();
 87  for (uint64_t i = 0; i < nn; ++i) {
 88    av[i] = av[i] - bv[i];
 89  }
 90  return a;
 91}
 92
 93reim_fft64vec simple_fft64(const znx_i64& polynomial) {
 94  const uint64_t nn = polynomial.nn();
 95  const uint64_t m = nn / 2;
 96  reim_fft64vec res(nn);
 97  double* dat = res.data();
 98  for (uint64_t i = 0; i < nn; ++i) dat[i] = polynomial.get_coeff(i);
 99  reim_fft_simple(m, dat);
100  return res;
101}
102
103znx_i64 simple_rint_ifft64(const reim_fft64vec& fftvec) {
104  const uint64_t nn = fftvec.nn();
105  const uint64_t m = nn / 2;
106  std::vector<double> vv(fftvec.data(), fftvec.data() + nn);
107  double* v = vv.data();
108  reim_ifft_simple(m, v);
109  znx_i64 res(nn);
110  for (uint64_t i = 0; i < nn; ++i) {
111    res.set_coeff(i, rint(v[i] / m));
112  }
113  return res;
114}
115
116rnx_f64 naive_ifft64(const reim_fft64vec& fftvec) {
117  const uint64_t nn = fftvec.nn();
118  const uint64_t m = nn / 2;
119  std::vector<double> vv(fftvec.data(), fftvec.data() + nn);
120  double* v = vv.data();
121  reim_ifft_simple(m, v);
122  rnx_f64 res(nn);
123  for (uint64_t i = 0; i < nn; ++i) {
124    res.set_coeff(i, v[i] / m);
125  }
126  return res;
127}
128double infty_dist(const reim_fft64vec& a, const reim_fft64vec& b) {
129  const uint64_t n = a.nn();
130  REQUIRE_DRAMATICALLY(b.nn() == a.nn(), "dimensions mismatch");
131  const double* da = a.data();
132  const double* db = b.data();
133  double d = 0;
134  for (uint64_t i = 0; i < n; ++i) {
135    double di = abs(da[i] - db[i]);
136    if (di > d) d = di;
137  }
138  return d;
139}
140
141reim_fft64vec simple_fft64(const rnx_f64& polynomial) {
142  const uint64_t nn = polynomial.nn();
143  const uint64_t m = nn / 2;
144  reim_fft64vec res(nn);
145  double* dat = res.data();
146  for (uint64_t i = 0; i < nn; ++i) dat[i] = polynomial.get_coeff(i);
147  reim_fft_simple(m, dat);
148  return res;
149}
150
151reim_fft64vec operator*(double coeff, const reim_fft64vec& v) {
152  const uint64_t nn = v.nn();
153  reim_fft64vec res(nn);
154  double* rr = res.data();
155  const double* vv = v.data();
156  for (uint64_t i = 0; i < nn; ++i) rr[i] = coeff * vv[i];
157  return res;
158}
159
160rnx_f64 simple_ifft64(const reim_fft64vec& v) {
161  const uint64_t nn = v.nn();
162  const uint64_t m = nn / 2;
163  rnx_f64 res(nn);
164  double* dat = res.data();
165  memcpy(dat, v.data(), nn * sizeof(double));
166  reim_ifft_simple(m, dat);
167  return res;
168}