Commit 598e9765 authored by Martin Reinecke's avatar Martin Reinecke

more cleanups

parent 30c887d1
......@@ -198,6 +198,8 @@ template<typename T> struct cmplx {
cmplx(T r_, T i_) : r(r_), i(i_) {}
void Set(T r_, T i_) { r=r_; i=i_; }
void Set(T r_) { r=r_; i=T(0); }
void Split(T &r_, T &i_) const { r_=r; i_=i; }
void SplitConj(T &r_, T &i_) const { r_=r; i_=-i; }
cmplx &operator+= (const cmplx &other)
{ r+=other.r; i+=other.i; return *this; }
template<typename T2>cmplx &operator*= (T2 other)
......@@ -2348,10 +2350,7 @@ template<typename T0> class T_dst1
arr<T> tmp(N);
tmp[0] = tmp[n+1] = c[0]*0;
for (size_t i=0; i<n; ++i)
{
tmp[i+1] = c[i];
tmp[N-1-i] = -c[i];
}
{ tmp[i+1]=c[i]; tmp[N-1-i]=-c[i]; }
fftplan.forward(tmp.data(), fct);
for (size_t i=0; i<n; ++i)
c[i] = -tmp[2*i+2];
......@@ -2366,63 +2365,6 @@ template<typename T0> class T_dcst23
pocketfft_r<T0> fftplan;
vector<T0> twiddle;
template<typename T> POCKETFFT_NOINLINE void exec2(T c[], T0 fct,
bool ortho, bool cosine) const
{
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length();
size_t NS2 = (N+1)/2;
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
for (size_t k=1; k<N-1; k+=2)
{
T tmp = T0(0.5)*(c[k]+c[k+1]);
c[k+1] = T0(0.5)*(c[k+1]-c[k]);
c[k] = tmp;
}
fftplan.backward(c, fct);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1 = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
T t2 = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = t1+t2; c[kc]=t1-t2;
}
if ((N&1)==0)
c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
c[0] *= 2;
if (!cosine)
for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
swap(c[k], c[kc]);
if (ortho) c[0]/=sqrt2;
}
template<typename T> POCKETFFT_NOINLINE void exec3(T c[], T0 fct,
bool ortho, bool cosine) const
{
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length();
size_t NS2 = (N+1)/2;
if (ortho) c[0]*=sqrt2;
if (!cosine)
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
swap(c[k], c[kc]);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1=c[k]+c[kc], t2=c[k]-c[kc];
c[k] = twiddle[k-1]*t2+twiddle[kc-1]*t1;
c[kc]= twiddle[k-1]*t1-twiddle[kc-1]*t2;
}
if ((N&1)==0)
c[NS2] = 2*twiddle[NS2-1]*c[NS2];
fftplan.forward(c, fct);
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k], c[k+1]);
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
}
public:
POCKETFFT_NOINLINE T_dcst23(size_t length)
: fftplan(length), twiddle(length)
......@@ -2434,7 +2376,55 @@ template<typename T0> class T_dcst23
template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho,
int type, bool cosine) const
{ type==2 ? exec2(c, fct, ortho, cosine) : exec3(c, fct, ortho, cosine); }
{
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length();
size_t NS2 = (N+1)/2;
if (type==2)
{
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
c[0] *= 2;
if ((N&1)==0) c[N-1]*=2;
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k+1], c[k]);
fftplan.backward(c, fct);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1 = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
T t2 = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = T0(0.5)*(t1+t2); c[kc]=T0(0.5)*(t1-t2);
}
if ((N&1)==0)
c[NS2] *= twiddle[NS2-1];
if (!cosine)
for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
swap(c[k], c[kc]);
if (ortho) c[0]/=sqrt2;
}
else
{
if (ortho) c[0]*=sqrt2;
if (!cosine)
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
swap(c[k], c[kc]);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1=c[k]+c[kc], t2=c[k]-c[kc];
c[k] = twiddle[k-1]*t2+twiddle[kc-1]*t1;
c[kc]= twiddle[k-1]*t1-twiddle[kc-1]*t2;
}
if ((N&1)==0)
c[NS2] *= 2*twiddle[NS2-1];
fftplan.forward(c, fct);
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k], c[k+1]);
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
}
}
size_t length() const { return fftplan.length(); }
};
......@@ -3120,17 +3110,11 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
if (forward)
for (; i<len-1; i+=2, ++ii)
for (size_t j=0; j<vlen; ++j)
{
tdatav[i ][j] = in[it.iofs(j,ii)].r;
tdatav[i+1][j] = -in[it.iofs(j,ii)].i;
}
in[it.iofs(j,ii)].SplitConj(tdatav[i][j], tdatav[i+1][j]);
else
for (; i<len-1; i+=2, ++ii)
for (size_t j=0; j<vlen; ++j)
{
tdatav[i ][j] = in[it.iofs(j,ii)].r;
tdatav[i+1][j] = in[it.iofs(j,ii)].i;
}
in[it.iofs(j,ii)].Split(tdatav[i][j], tdatav[i+1][j]);
if (i<len)
for (size_t j=0; j<vlen; ++j)
tdatav[i][j] = in[it.iofs(j,ii)].r;
......@@ -3150,16 +3134,10 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
size_t i=1, ii=1;
if (forward)
for (; i<len-1; i+=2, ++ii)
{
tdata[i ] = in[it.iofs(ii)].r;
tdata[i+1] = -in[it.iofs(ii)].i;
}
in[it.iofs(ii)].SplitConj(tdata[i], tdata[i+1]);
else
for (; i<len-1; i+=2, ++ii)
{
tdata[i ] = in[it.iofs(ii)].r;
tdata[i+1] = in[it.iofs(ii)].i;
}
in[it.iofs(ii)].Split(tdata[i], tdata[i+1]);
if (i<len)
tdata[i] = in[it.iofs(ii)].r;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment