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
 10EXPORT void reim4_from_cplx_ref(const REIM4_FROM_CPLX_PRECOMP* tables, double* r, const void* a) {
 11  const double* x = (double*)a;
 12  const uint64_t m = tables->m;
 13  for (uint64_t i = 0; i < m / 4; ++i) {
 14    double r0 = x[0];
 15    double i0 = x[1];
 16    double r1 = x[2];
 17    double i1 = x[3];
 18    double r2 = x[4];
 19    double i2 = x[5];
 20    double r3 = x[6];
 21    double i3 = x[7];
 22    r[0] = r0;
 23    r[1] = r2;
 24    r[2] = r1;
 25    r[3] = r3;
 26    r[4] = i0;
 27    r[5] = i2;
 28    r[6] = i1;
 29    r[7] = i3;
 30    x += 8;
 31    r += 8;
 32  }
 33}
 34
 35void* init_reim4_from_cplx_precomp(REIM4_FROM_CPLX_PRECOMP* res, uint32_t nn) {
 36  res->m = nn / 2;
 37  if (CPU_SUPPORTS("fma")) {
 38    if (nn >= 4) {
 39      res->function = reim4_from_cplx_fma;
 40    } else {
 41      res->function = reim4_from_cplx_ref;
 42    }
 43  } else {
 44    res->function = reim4_from_cplx_ref;
 45  }
 46  return res;
 47}
 48
 49EXPORT REIM4_FROM_CPLX_PRECOMP* new_reim4_from_cplx_precomp(uint32_t m) {
 50  REIM4_FROM_CPLX_PRECOMP* res = malloc(sizeof(*res));
 51  if (!res) return spqlios_error(strerror(errno));
 52  return spqlios_keep_or_free(res, init_reim4_from_cplx_precomp(res, m));
 53}
 54
 55EXPORT void reim4_from_cplx_simple(uint32_t m, double* r, const void* a) {
 56  static REIM4_FROM_CPLX_PRECOMP precomp[32];
 57  REIM4_FROM_CPLX_PRECOMP* p = precomp + log2m(m);
 58  if (!p->function) {
 59    if (!init_reim4_from_cplx_precomp(p, m)) abort();
 60  }
 61  p->function(p, r, a);
 62}
 63
 64EXPORT void reim4_to_cplx_ref(const REIM4_TO_CPLX_PRECOMP* tables, void* r, const double* a) {
 65  double* y = (double*)r;
 66  const uint64_t m = tables->m;
 67  for (uint64_t i = 0; i < m / 4; ++i) {
 68    double r0 = a[0];
 69    double r2 = a[1];
 70    double r1 = a[2];
 71    double r3 = a[3];
 72    double i0 = a[4];
 73    double i2 = a[5];
 74    double i1 = a[6];
 75    double i3 = a[7];
 76    y[0] = r0;
 77    y[1] = i0;
 78    y[2] = r1;
 79    y[3] = i1;
 80    y[4] = r2;
 81    y[5] = i2;
 82    y[6] = r3;
 83    y[7] = i3;
 84    a += 8;
 85    y += 8;
 86  }
 87}
 88
 89void* init_reim4_to_cplx_precomp(REIM4_TO_CPLX_PRECOMP* res, uint32_t m) {
 90  res->m = m;
 91  if (CPU_SUPPORTS("fma")) {
 92    if (m >= 2) {
 93      res->function = reim4_to_cplx_fma;
 94    } else {
 95      res->function = reim4_to_cplx_ref;
 96    }
 97  } else {
 98    res->function = reim4_to_cplx_ref;
 99  }
100  return res;
101}
102
103EXPORT REIM4_TO_CPLX_PRECOMP* new_reim4_to_cplx_precomp(uint32_t m) {
104  REIM4_TO_CPLX_PRECOMP* res = malloc(sizeof(*res));
105  if (!res) return spqlios_error(strerror(errno));
106  return spqlios_keep_or_free(res, init_reim4_to_cplx_precomp(res, m));
107}
108
109EXPORT void reim4_to_cplx_simple(uint32_t m, void* r, const double* a) {
110  static REIM4_TO_CPLX_PRECOMP precomp[32];
111  REIM4_TO_CPLX_PRECOMP* p = precomp + log2m(m);
112  if (!p->function) {
113    if (!init_reim4_to_cplx_precomp(p, m)) abort();
114  }
115  p->function(p, r, a);
116}