1#include <string.h>
 2
 3#include "../commons_private.h"
 4#include "cplx_fft_internal.h"
 5#include "cplx_fft_private.h"
 6
 7EXPORT void cplx_fftvec_addmul_ref(const CPLX_FFTVEC_ADDMUL_PRECOMP* precomp, void* r, const void* a, const void* b) {
 8  const uint32_t m = precomp->m;
 9  const CPLX* aa = (CPLX*)a;
10  const CPLX* bb = (CPLX*)b;
11  CPLX* rr = (CPLX*)r;
12  for (uint32_t i = 0; i < m; ++i) {
13    const double re = aa[i][0] * bb[i][0] - aa[i][1] * bb[i][1];
14    const double im = aa[i][0] * bb[i][1] + aa[i][1] * bb[i][0];
15    rr[i][0] += re;
16    rr[i][1] += im;
17  }
18}
19
20EXPORT void cplx_fftvec_mul_ref(const CPLX_FFTVEC_MUL_PRECOMP* precomp, void* r, const void* a, const void* b) {
21  const uint32_t m = precomp->m;
22  const CPLX* aa = (CPLX*)a;
23  const CPLX* bb = (CPLX*)b;
24  CPLX* rr = (CPLX*)r;
25  for (uint32_t i = 0; i < m; ++i) {
26    const double re = aa[i][0] * bb[i][0] - aa[i][1] * bb[i][1];
27    const double im = aa[i][0] * bb[i][1] + aa[i][1] * bb[i][0];
28    rr[i][0] = re;
29    rr[i][1] = im;
30  }
31}
32
33EXPORT void* init_cplx_fftvec_addmul_precomp(CPLX_FFTVEC_ADDMUL_PRECOMP* r, uint32_t m) {
34  if (m & (m - 1)) return spqlios_error("m must be a power of two");
35  r->m = m;
36  if (m <= 4) {
37    r->function = cplx_fftvec_addmul_ref;
38  } else if (CPU_SUPPORTS("fma")) {
39    r->function = cplx_fftvec_addmul_fma;
40  } else {
41    r->function = cplx_fftvec_addmul_ref;
42  }
43  return r;
44}
45
46EXPORT void* init_cplx_fftvec_mul_precomp(CPLX_FFTVEC_MUL_PRECOMP* r, uint32_t m) {
47  if (m & (m - 1)) return spqlios_error("m must be a power of two");
48  r->m = m;
49  if (m <= 4) {
50    r->function = cplx_fftvec_mul_ref;
51  } else if (CPU_SUPPORTS("fma")) {
52    r->function = cplx_fftvec_mul_fma;
53  } else {
54    r->function = cplx_fftvec_mul_ref;
55  }
56  return r;
57}
58
59EXPORT CPLX_FFTVEC_ADDMUL_PRECOMP* new_cplx_fftvec_addmul_precomp(uint32_t m) {
60  CPLX_FFTVEC_ADDMUL_PRECOMP* r = malloc(sizeof(CPLX_FFTVEC_MUL_PRECOMP));
61  return spqlios_keep_or_free(r, init_cplx_fftvec_addmul_precomp(r, m));
62}
63
64EXPORT CPLX_FFTVEC_MUL_PRECOMP* new_cplx_fftvec_mul_precomp(uint32_t m) {
65  CPLX_FFTVEC_MUL_PRECOMP* r = malloc(sizeof(CPLX_FFTVEC_MUL_PRECOMP));
66  return spqlios_keep_or_free(r, init_cplx_fftvec_mul_precomp(r, m));
67}
68
69EXPORT void cplx_fftvec_mul_simple(uint32_t m, void* r, const void* a, const void* b) {
70  static CPLX_FFTVEC_MUL_PRECOMP p[31] = {0};
71  CPLX_FFTVEC_MUL_PRECOMP* f = p + log2m(m);
72  if (!f->function) {
73    if (!init_cplx_fftvec_mul_precomp(f, m)) abort();
74  }
75  f->function(f, r, a, b);
76}
77
78EXPORT void cplx_fftvec_addmul_simple(uint32_t m, void* r, const void* a, const void* b) {
79  static CPLX_FFTVEC_ADDMUL_PRECOMP p[31] = {0};
80  CPLX_FFTVEC_ADDMUL_PRECOMP* f = p + log2m(m);
81  if (!f->function) {
82    if (!init_cplx_fftvec_addmul_precomp(f, m)) abort();
83  }
84  f->function(f, r, a, b);
85}