1#include <gtest/gtest.h>
 2
 3#include <cmath>
 4
 5#include "spqlios/cplx/cplx_fft_internal.h"
 6#include "spqlios/cplx/cplx_fft_private.h"
 7
 8#ifdef __x86_64__
 9TEST(fft, cplx_from_znx32_ref_vs_fma) {
10  const uint32_t m = 128;
11  int32_t* src = (int32_t*)spqlios_alloc_custom_align(32, 10 * m * sizeof(int32_t));
12  CPLX* dst1 = (CPLX*)(src + 2 * m);
13  CPLX* dst2 = (CPLX*)(src + 6 * m);
14  for (uint64_t i = 0; i < 2 * m; ++i) {
15    src[i] = rand() - RAND_MAX / 2;
16  }
17  CPLX_FROM_ZNX32_PRECOMP precomp;
18  precomp.m = m;
19  cplx_from_znx32_ref(&precomp, dst1, src);
20  // cplx_from_znx32_simple(m, 32, dst1, src);
21  cplx_from_znx32_avx2_fma(&precomp, dst2, src);
22  for (uint64_t i = 0; i < m; ++i) {
23    ASSERT_EQ(dst1[i][0], dst2[i][0]);
24    ASSERT_EQ(dst1[i][1], dst2[i][1]);
25  }
26  spqlios_free(src);
27}
28#endif
29
30#ifdef __x86_64__
31TEST(fft, cplx_from_tnx32_ref_vs_fma) {
32  const uint32_t m = 128;
33  int32_t* src = (int32_t*)spqlios_alloc_custom_align(32, 10 * m * sizeof(int32_t));
34  CPLX* dst1 = (CPLX*)(src + 2 * m);
35  CPLX* dst2 = (CPLX*)(src + 6 * m);
36  for (uint64_t i = 0; i < 2 * m; ++i) {
37    src[i] = rand() + (rand() << 20);
38  }
39  CPLX_FROM_TNX32_PRECOMP precomp;
40  precomp.m = m;
41  cplx_from_tnx32_ref(&precomp, dst1, src);
42  // cplx_from_tnx32_simple(m, dst1, src);
43  cplx_from_tnx32_avx2_fma(&precomp, dst2, src);
44  for (uint64_t i = 0; i < m; ++i) {
45    ASSERT_EQ(dst1[i][0], dst2[i][0]);
46    ASSERT_EQ(dst1[i][1], dst2[i][1]);
47  }
48  spqlios_free(src);
49}
50#endif
51
52#ifdef __x86_64__
53TEST(fft, cplx_to_tnx32_ref_vs_fma) {
54  for (const uint32_t m : {8, 128, 1024, 65536}) {
55    for (const double divisor : {double(1), double(m), double(0.5)}) {
56      CPLX* src = (CPLX*)spqlios_alloc_custom_align(32, 10 * m * sizeof(int32_t));
57      int32_t* dst1 = (int32_t*)(src + m);
58      int32_t* dst2 = (int32_t*)(src + 2 * m);
59      for (uint64_t i = 0; i < 2 * m; ++i) {
60        src[i][0] = (rand() / double(RAND_MAX) - 0.5) * pow(2., 19 - (rand() % 60)) * divisor;
61        src[i][1] = (rand() / double(RAND_MAX) - 0.5) * pow(2., 19 - (rand() % 60)) * divisor;
62      }
63      CPLX_TO_TNX32_PRECOMP precomp;
64      precomp.m = m;
65      precomp.divisor = divisor;
66      cplx_to_tnx32_ref(&precomp, dst1, src);
67      cplx_to_tnx32_avx2_fma(&precomp, dst2, src);
68      // cplx_to_tnx32_simple(m, divisor, 18, dst2, src);
69      for (uint64_t i = 0; i < 2 * m; ++i) {
70        double truevalue =
71            (src[i % m][i / m] / divisor - floor(src[i % m][i / m] / divisor + 0.5)) * (INT64_C(1) << 32);
72        if (fabs(truevalue - floor(truevalue)) == 0.5) {
73          // ties can differ by 0, 1 or -1
74          ASSERT_LE(abs(dst1[i] - dst2[i]), 0)
75              << i << " " << dst1[i] << " " << dst2[i] << " " << truevalue << std::endl;
76        } else {
77          // otherwise, we should have equality
78          ASSERT_LE(abs(dst1[i] - dst2[i]), 0)
79              << i << " " << dst1[i] << " " << dst2[i] << " " << truevalue << std::endl;
80        }
81      }
82      spqlios_free(src);
83    }
84  }
85}
86#endif