/*

* Copyright (c) 2003, 2007-14 Matteo Frigo

* Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology

*

* 128-bit AVX support by Erik Lindahl, 2015.

* Erik Lindahl hereby places his modifications in the public domain.

*

* This program is free software; you can redistribute it and/or modify

* it under the terms of the GNU General Public License as published by

* the Free Software Foundation; either version 2 of the License, or

* (at your option) any later version.

*

* This program is distributed in the hope that it will be useful,

* but WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

* GNU General Public License for more details.

*

* You should have received a copy of the GNU General Public License

* along with this program; if not, write to the Free Software

* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

*

*/

#if defined(FFTW_LDOUBLE) || defined(FFTW_QUAD)

#error "AVX only works in single or double precision"

#endif

#ifdef FFTW_SINGLE

# define DS(d,s) s /* single-precision option */

# define SUFF(name) name ## s

#else

# define DS(d,s) d /* double-precision option */

# define SUFF(name) name ## d

#endif

#define SIMD_SUFFIX _avx_128_fma /* for renaming */

#define VL DS(1,2) /* SIMD vector length, in term of complex numbers */

#define SIMD_VSTRIDE_OKA(x) DS(1,((x) == 2))

#define SIMD_STRIDE_OKPAIR SIMD_STRIDE_OK

#ifdef _MSC_VER

#ifndef inline

#define inline __inline

#endif

#endif

#include <immintrin.h>

#ifdef _MSC_VER

# include <intrin.h>

#elif defined (__GNUC__)

# include <x86intrin.h>

#endif

#if !(defined(__AVX__) && defined(__FMA4__)) /* sanity check */

#error "compiling simd-avx-128-fma.h without -mavx or -mfma4"

#endif

typedef DS(__m128d,__m128) V;

#define VADD SUFF(_mm_add_p)

#define VSUB SUFF(_mm_sub_p)

#define VMUL SUFF(_mm_mul_p)

#define VXOR SUFF(_mm_xor_p)

#define SHUF SUFF(_mm_shuffle_p)

#define VPERM1 SUFF(_mm_permute_p)

#define UNPCKL SUFF(_mm_unpacklo_p)

#define UNPCKH SUFF(_mm_unpackhi_p)

#define SHUFVALS(fp0,fp1,fp2,fp3) \

(((fp3) << 6) | ((fp2) << 4) | ((fp1) << 2) | ((fp0)))

#define VDUPL(x) DS(_mm_permute_pd(x,0), _mm_moveldup_ps(x))

#define VDUPH(x) DS(_mm_permute_pd(x,3), _mm_movehdup_ps(x))

#define LOADH(addr, val) _mm_loadh_pi(val, (const __m64 *)(addr))

#define LOADL(addr, val) _mm_loadl_pi(val, (const __m64 *)(addr))

#define STOREH(a, v) DS(_mm_storeh_pd(a, v), _mm_storeh_pi((__m64 *)(a), v))

#define STOREL(a, v) DS(_mm_storel_pd(a, v), _mm_storel_pi((__m64 *)(a), v))

#define VLIT(x0, x1) DS(_mm_set_pd(x0, x1), _mm_set_ps(x0, x1, x0, x1))

#define DVK(var, val) V var = VLIT(val, val)

#define LDK(x) x

static inline V LDA(const R *x, INT ivs, const R *aligned_like)

{

(void)aligned_like; /* UNUSED */

(void)ivs; /* UNUSED */

return *(const V *)x;

}

static inline void STA(R *x, V v, INT ovs, const R *aligned_like)

{

(void)aligned_like; /* UNUSED */

(void)ovs; /* UNUSED */

*(V *)x = v;

}

#ifdef FFTW_SINGLE

static inline V LD(const R *x, INT ivs, const R *aligned_like)

{

V var;

#if defined(__ICC) || (__GNUC__ > 4) || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)

var = LOADL(x, SUFF(_mm_undefined_p)());

var = LOADH(x + ivs, var);

#else

var = LOADL(x, var);

var = LOADH(x + ivs, var);

#endif

return var;

}

# ifdef _MSC_VER

# pragma warning(default : 4700)

# pragma runtime_checks("u", restore)

# endif

static inline void ST(R *x, V v, INT ovs, const R *aligned_like)

{

(void)aligned_like; /* UNUSED */

/* WARNING: the extra_iter hack depends upon STOREL occurring

after STOREH */

STOREH(x + ovs, v);

STOREL(x, v);

}

#else /* ! FFTW_SINGLE */

# define LD LDA

# define ST STA

#endif

#define STM2 DS(STA,ST)

#define STN2(x, v0, v1, ovs) /* nop */

#ifdef FFTW_SINGLE

# define STM4(x, v, ovs, aligned_like) /* no-op */

/* STN4 is a macro, not a function, thanks to Visual C++ developers

deciding "it would be infrequent that people would want to pass more

than 3 [__m128 parameters] by value." 3 parameters ought to be enough

for anybody. */

# define STN4(x, v0, v1, v2, v3, ovs) \

{ \

V xxx0, xxx1, xxx2, xxx3; \

xxx0 = UNPCKL(v0, v2); \

xxx1 = UNPCKH(v0, v2); \

xxx2 = UNPCKL(v1, v3); \

xxx3 = UNPCKH(v1, v3); \

STA(x, UNPCKL(xxx0, xxx2), 0, 0); \

STA(x + ovs, UNPCKH(xxx0, xxx2), 0, 0); \

STA(x + 2 * ovs, UNPCKL(xxx1, xxx3), 0, 0); \

STA(x + 3 * ovs, UNPCKH(xxx1, xxx3), 0, 0); \

}

#else /* !FFTW_SINGLE */

static inline void STM4(R *x, V v, INT ovs, const R *aligned_like)

{

(void)aligned_like; /* UNUSED */

STOREL(x, v);

STOREH(x + ovs, v);

}

# define STN4(x, v0, v1, v2, v3, ovs) /* nothing */

#endif

static inline V FLIP_RI(V x)

{

return VPERM1(x, DS(1, SHUFVALS(1, 0, 3, 2)));

}

static inline V VCONJ(V x)

{

/* Produce a SIMD vector[VL] of (0 + -0i).

We really want to write this:

V pmpm = VLIT(-0.0, 0.0);

but historically some compilers have ignored the distiction

between +0 and -0. It looks like 'gcc-8 -fast-math' treats -0

as 0 too.

*/

union uvec {

unsigned u[4];

V v;

};

static const union uvec pmpm = {

};

return VXOR(pmpm.v, x);

}

static inline V VBYI(V x)

{

x = VCONJ(x);

x = FLIP_RI(x);

return x;

}

/* FMA support */

#define VFMA(a, b, c) SUFF(_mm_macc_p)(a,b,c)

#define VFNMS(a, b, c) SUFF(_mm_nmacc_p)(a,b,c)

#define VFMS(a, b, c) SUFF(_mm_msub_p)(a,b,c)

#define VFMAI(b, c) SUFF(_mm_addsub_p)(c,FLIP_RI(b))

#define VFNMSI(b, c) VSUB(c, VBYI(b))

#define VFMACONJ(b,c) VADD(VCONJ(b),c)

#define VFMSCONJ(b,c) VSUB(VCONJ(b),c)

#define VFNMSCONJ(b,c) SUFF(_mm_addsub_p)(c,b)

static inline V VZMUL(V tx, V sr)

{

V tr = VDUPL(tx);

V ti = VDUPH(tx);

tr = VMUL(tr, sr);

ti = VMUL(ti, FLIP_RI(sr));

return SUFF(_mm_addsub_p)(tr,ti);

}

static inline V VZMULJ(V tx, V sr)

{

V tr = VDUPL(tx);

V ti = VDUPH(tx);

tr = VMUL(tr, sr);

sr = VBYI(sr);

return VFNMS(ti, sr, tr);

}

static inline V VZMULI(V tx, V sr)

{

V tr = VDUPL(tx);

V ti = VDUPH(tx);

ti = VMUL(ti, sr);

sr = VBYI(sr);

return VFMS(tr, sr, ti);

}

static inline V VZMULIJ(V tx, V sr)

{

V tr = VDUPL(tx);

V ti = VDUPH(tx);

ti = VMUL(ti, sr);

tr = VMUL(tr, FLIP_RI(sr));

return SUFF(_mm_addsub_p)(ti,tr);

}

/* twiddle storage #1: compact, slower */

#ifdef FFTW_SINGLE

# define VTW1(v,x) \

{TW_COS, v, x}, {TW_COS, v+1, x}, {TW_SIN, v, x}, {TW_SIN, v+1, x}

static inline V BYTW1(const R *t, V sr)

{

const V *twp = (const V *)t;

V tx = twp[0];

V tr = UNPCKL(tx, tx);

V ti = UNPCKH(tx, tx);

tr = VMUL(tr, sr);

ti = VMUL(ti, FLIP_RI(sr));

return SUFF(_mm_addsub_p)(tr,ti);

}

static inline V BYTWJ1(const R *t, V sr)

{

const V *twp = (const V *)t;

V tx = twp[0];

V tr = UNPCKL(tx, tx);

V ti = UNPCKH(tx, tx);

tr = VMUL(tr, sr);

sr = VBYI(sr);

return VFNMS(ti, sr, tr);

}

#else /* !FFTW_SINGLE */

# define VTW1(v,x) {TW_CEXP, v, x}

static inline V BYTW1(const R *t, V sr)

{

V tx = LD(t, 1, t);

return VZMUL(tx, sr);

}

static inline V BYTWJ1(const R *t, V sr)

{

V tx = LD(t, 1, t);

return VZMULJ(tx, sr);

}

#endif

#define TWVL1 (VL)

/* twiddle storage #2: twice the space, faster (when in cache) */

#ifdef FFTW_SINGLE

# define VTW2(v,x) \

{TW_COS, v, x}, {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+1, x}, \

{TW_SIN, v, -x}, {TW_SIN, v, x}, {TW_SIN, v+1, -x}, {TW_SIN, v+1, x}

#else /* !FFTW_SINGLE */

# define VTW2(v,x) \

{TW_COS, v, x}, {TW_COS, v, x}, {TW_SIN, v, -x}, {TW_SIN, v, x}

#endif

#define TWVL2 (2 * VL)

static inline V BYTW2(const R *t, V sr)

{

const V *twp = (const V *)t;

V si = FLIP_RI(sr);

V tr = twp[0], ti = twp[1];

return VFMA(tr, sr, VMUL(ti, si));

}

static inline V BYTWJ2(const R *t, V sr)

{

const V *twp = (const V *)t;

V si = FLIP_RI(sr);

V tr = twp[0], ti = twp[1];

return VFNMS(ti, si, VMUL(tr, sr));

}

/* twiddle storage #3 */

#ifdef FFTW_SINGLE

# define VTW3(v,x) {TW_CEXP, v, x}, {TW_CEXP, v+1, x}

# define TWVL3 (VL)

#else

# define VTW3(v,x) VTW1(v,x)

# define TWVL3 TWVL1

#endif

/* twiddle storage for split arrays */

#ifdef FFTW_SINGLE

# define VTWS(v,x) \

{TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+2, x}, {TW_COS, v+3, x}, \

{TW_SIN, v, x}, {TW_SIN, v+1, x}, {TW_SIN, v+2, x}, {TW_SIN, v+3, x}

#else

# define VTWS(v,x) \

{TW_COS, v, x}, {TW_COS, v+1, x}, {TW_SIN, v, x}, {TW_SIN, v+1, x}

#endif

#define TWVLS (2 * VL)

#define VLEAVE() /* nothing */

#include "simd-common.h"