1#include <gtest/gtest.h>
  2#include <spqlios/reim/reim_fft_internal.h>
  3
  4#include "testlib/test_commons.h"
  5
  6TEST(reim_conversions, reim_to_tnx) {
  7  for (uint32_t m : {1, 2, 64, 128, 512}) {
  8    for (double divisor : {1, 2, int(m)}) {
  9      for (uint32_t log2overhead : {1, 2, 10, 18, 35, 42}) {
 10        double maxdiff = pow(2., log2overhead - 50);
 11        std::vector<double> data(2 * m);
 12        std::vector<double> dout(2 * m);
 13        for (uint64_t i = 0; i < 2 * m; ++i) {
 14          data[i] = (uniform_f64_01() - 0.5) * pow(2., log2overhead + 1) * divisor;
 15        }
 16        REIM_TO_TNX_PRECOMP* p = new_reim_to_tnx_precomp(m, divisor, 18);
 17        reim_to_tnx_ref(p, dout.data(), data.data());
 18        for (uint64_t i = 0; i < 2 * m; ++i) {
 19          ASSERT_LE(fabs(dout[i]), 0.5);
 20          double diff = dout[i] - data[i] / divisor;
 21          double fracdiff = diff - rint(diff);
 22          ASSERT_LE(fabs(fracdiff), maxdiff);
 23        }
 24        delete_reim_to_tnx_precomp(p);
 25      }
 26    }
 27  }
 28}
 29
 30#ifdef __x86_64__
 31TEST(reim_conversions, reim_to_tnx_ref_vs_avx) {
 32  for (uint32_t m : {8, 16, 64, 128, 512}) {
 33    for (double divisor : {1, 2, int(m)}) {
 34      for (uint32_t log2overhead : {1, 2, 10, 18, 35, 42}) {
 35        // double maxdiff = pow(2., log2overhead - 50);
 36        std::vector<double> data(2 * m);
 37        std::vector<double> dout1(2 * m);
 38        std::vector<double> dout2(2 * m);
 39        for (uint64_t i = 0; i < 2 * m; ++i) {
 40          data[i] = (uniform_f64_01() - 0.5) * pow(2., log2overhead + 1) * divisor;
 41        }
 42        REIM_TO_TNX_PRECOMP* p = new_reim_to_tnx_precomp(m, divisor, 18);
 43        reim_to_tnx_ref(p, dout1.data(), data.data());
 44        reim_to_tnx_avx(p, dout2.data(), data.data());
 45        for (uint64_t i = 0; i < 2 * m; ++i) {
 46          ASSERT_LE(fabs(dout1[i] - dout2[i]), 0.);
 47        }
 48        delete_reim_to_tnx_precomp(p);
 49      }
 50    }
 51  }
 52}
 53#endif
 54
 55typedef typeof(reim_from_znx64_ref) reim_from_znx64_f;
 56
 57void test_reim_from_znx64(reim_from_znx64_f reim_from_znx64, uint64_t maxbnd) {
 58  for (uint32_t m : {4, 8, 16, 64, 16384}) {
 59    REIM_FROM_ZNX64_PRECOMP* p = new_reim_from_znx64_precomp(m, maxbnd);
 60    std::vector<int64_t> data(2 * m);
 61    std::vector<double> dout(2 * m);
 62    for (uint64_t i = 0; i < 2 * m; ++i) {
 63      int64_t magnitude = int64_t(uniform_u64() % (maxbnd + 1));
 64      data[i] = uniform_i64() >> (63 - magnitude);
 65      REQUIRE_DRAMATICALLY(abs(data[i]) <= (INT64_C(1) << magnitude), "pb");
 66    }
 67    reim_from_znx64(p, dout.data(), data.data());
 68    for (uint64_t i = 0; i < 2 * m; ++i) {
 69      ASSERT_EQ(dout[i], double(data[i])) << dout[i] << " " << data[i];
 70    }
 71    delete_reim_from_znx64_precomp(p);
 72  }
 73}
 74
 75TEST(reim_conversions, reim_from_znx64) {
 76  for (uint64_t maxbnd : {50}) {
 77    test_reim_from_znx64(reim_from_znx64, maxbnd);
 78  }
 79}
 80TEST(reim_conversions, reim_from_znx64_ref) { test_reim_from_znx64(reim_from_znx64_ref, 50); }
 81#ifdef __x86_64__
 82TEST(reim_conversions, reim_from_znx64_avx2_bnd50_fma) { test_reim_from_znx64(reim_from_znx64_bnd50_fma, 50); }
 83#endif
 84
 85typedef typeof(reim_to_znx64_ref) reim_to_znx64_f;
 86
 87void test_reim_to_znx64(reim_to_znx64_f reim_to_znx64_fcn, int64_t maxbnd) {
 88  for (uint32_t m : {4, 8, 16, 64, 16384}) {
 89    for (double divisor : {1, 2, int(m)}) {
 90      REIM_TO_ZNX64_PRECOMP* p = new_reim_to_znx64_precomp(m, divisor, maxbnd);
 91      std::vector<double> data(2 * m);
 92      std::vector<int64_t> dout(2 * m);
 93      for (uint64_t i = 0; i < 2 * m; ++i) {
 94        int64_t magnitude = int64_t(uniform_u64() % (maxbnd + 11)) - 10;
 95        data[i] = (uniform_f64_01() - 0.5) * pow(2., magnitude + 1) * divisor;
 96      }
 97      reim_to_znx64_fcn(p, dout.data(), data.data());
 98      for (uint64_t i = 0; i < 2 * m; ++i) {
 99        ASSERT_LE(dout[i] - data[i] / divisor, 0.5) << dout[i] << " " << data[i];
100      }
101      delete_reim_to_znx64_precomp(p);
102    }
103  }
104}
105
106TEST(reim_conversions, reim_to_znx64) {
107  for (uint64_t maxbnd : {63, 50}) {
108    test_reim_to_znx64(reim_to_znx64, maxbnd);
109  }
110}
111TEST(reim_conversions, reim_to_znx64_ref) { test_reim_to_znx64(reim_to_znx64_ref, 63); }
112#ifdef __x86_64__
113TEST(reim_conversions, reim_to_znx64_avx2_bnd63_fma) { test_reim_to_znx64(reim_to_znx64_avx2_bnd63_fma, 63); }
114TEST(reim_conversions, reim_to_znx64_avx2_bnd50_fma) { test_reim_to_znx64(reim_to_znx64_avx2_bnd50_fma, 50); }
115#endif