1#include <errno.h>
 2#include <memory.h>
 3#include <stdlib.h>
 4#include <string.h>
 5
 6#include "../commons_private.h"
 7#include "reim4_fftvec_internal.h"
 8#include "reim4_fftvec_private.h"
 9
10void* init_reim4_fftvec_addmul_precomp(REIM4_FFTVEC_ADDMUL_PRECOMP* res, uint32_t m) {
11  res->m = m;
12  if (CPU_SUPPORTS("fma")) {
13    if (m >= 2) {
14      res->function = reim4_fftvec_addmul_fma;
15    } else {
16      res->function = reim4_fftvec_addmul_ref;
17    }
18  } else {
19    res->function = reim4_fftvec_addmul_ref;
20  }
21  return res;
22}
23
24EXPORT REIM4_FFTVEC_ADDMUL_PRECOMP* new_reim4_fftvec_addmul_precomp(uint32_t m) {
25  REIM4_FFTVEC_ADDMUL_PRECOMP* res = malloc(sizeof(*res));
26  if (!res) return spqlios_error(strerror(errno));
27  return spqlios_keep_or_free(res, init_reim4_fftvec_addmul_precomp(res, m));
28}
29
30EXPORT void reim4_fftvec_addmul_ref(const REIM4_FFTVEC_ADDMUL_PRECOMP* precomp, double* r, const double* a,
31                                    const double* b) {
32  const uint64_t m = precomp->m;
33  for (uint64_t j = 0; j < m / 4; ++j) {
34    for (uint64_t i = 0; i < 4; ++i) {
35      double re = a[i] * b[i] - a[i + 4] * b[i + 4];
36      double im = a[i] * b[i + 4] + a[i + 4] * b[i];
37      r[i] += re;
38      r[i + 4] += im;
39    }
40    a += 8;
41    b += 8;
42    r += 8;
43  }
44}
45
46EXPORT void reim4_fftvec_addmul_simple(uint32_t m, double* r, const double* a, const double* b) {
47  static REIM4_FFTVEC_ADDMUL_PRECOMP precomp[32];
48  REIM4_FFTVEC_ADDMUL_PRECOMP* p = precomp + log2m(m);
49  if (!p->function) {
50    if (!init_reim4_fftvec_addmul_precomp(p, m)) abort();
51  }
52  p->function(p, r, a, b);
53}
54
55void* init_reim4_fftvec_mul_precomp(REIM4_FFTVEC_MUL_PRECOMP* res, uint32_t m) {
56  res->m = m;
57  if (CPU_SUPPORTS("fma")) {
58    if (m >= 4) {
59      res->function = reim4_fftvec_mul_fma;
60    } else {
61      res->function = reim4_fftvec_mul_ref;
62    }
63  } else {
64    res->function = reim4_fftvec_mul_ref;
65  }
66  return res;
67}
68
69EXPORT REIM4_FFTVEC_MUL_PRECOMP* new_reim4_fftvec_mul_precomp(uint32_t m) {
70  REIM4_FFTVEC_MUL_PRECOMP* res = malloc(sizeof(*res));
71  if (!res) return spqlios_error(strerror(errno));
72  return spqlios_keep_or_free(res, init_reim4_fftvec_mul_precomp(res, m));
73}
74
75EXPORT void reim4_fftvec_mul_ref(const REIM4_FFTVEC_MUL_PRECOMP* precomp, double* r, const double* a, const double* b) {
76  const uint64_t m = precomp->m;
77  for (uint64_t j = 0; j < m / 4; ++j) {
78    for (uint64_t i = 0; i < 4; ++i) {
79      double re = a[i] * b[i] - a[i + 4] * b[i + 4];
80      double im = a[i] * b[i + 4] + a[i + 4] * b[i];
81      r[i] = re;
82      r[i + 4] = im;
83    }
84    a += 8;
85    b += 8;
86    r += 8;
87  }
88}
89
90EXPORT void reim4_fftvec_mul_simple(uint32_t m, double* r, const double* a, const double* b) {
91  static REIM4_FFTVEC_MUL_PRECOMP precomp[32];
92  REIM4_FFTVEC_MUL_PRECOMP* p = precomp + log2m(m);
93  if (!p->function) {
94    if (!init_reim4_fftvec_mul_precomp(p, m)) abort();
95  }
96  p->function(p, r, a, b);
97}