Commit 189bffbf authored by Martin Reinecke's avatar Martin Reinecke

don't store second half of unit circle

parent 61a3e544
......@@ -438,6 +438,7 @@ template<typename T> class sincos_2pibyn
{
private:
using Thigh = typename conditional<(sizeof(T)>sizeof(double)), T, double>::type;
size_t sz;
arr<cmplx<T>> data;
POCKETFFT_NOINLINE void calc_first_octant(size_t den,
......@@ -474,8 +475,8 @@ template<typename T> class sincos_2pibyn
void calc_first_quadrant(size_t n, cmplx<T> * POCKETFFT_RESTRICT res)
{
cmplx<T> * POCKETFFT_RESTRICT p = res+n/2; // n is always even here
calc_first_octant(n<<1, p);
arr<cmplx<T>> p(n/2); // n is always even here
calc_first_octant(n<<1, p.data());
size_t ndone=(n+2)>>2;
size_t i=0, idx1=0, idx2=ndone-1;
for (; i+1<ndone; i+=2, ++idx1, --idx2)
......@@ -490,8 +491,8 @@ template<typename T> class sincos_2pibyn
void calc_first_half(size_t n, cmplx<T> * POCKETFFT_RESTRICT res)
{
int ndone=int(n+1)>>1;
cmplx<T> * p = res+(n-1)/2; // n is always odd here
calc_first_octant(n<<2, p);
arr<cmplx<T>> p((n+1)/2); // n is always odd here
calc_first_octant(n<<2, p.data());
int i4=0, in=int(n), i=0;
for (; i4<=in-i4; ++i, i4+=4) // octant 0
res[i] = p[i4];
......@@ -522,16 +523,7 @@ template<typename T> class sincos_2pibyn
else
for (size_t i=1, j=half-1; 2*i<half; ++i, --j)
res[j].Set(-res[i].r, res[i].i);
}
void fill_second_half(size_t n, cmplx<T> * POCKETFFT_RESTRICT res)
{
if ((n&1)==0)
for (size_t i=0; i<n/2; ++i)
res[i+n/2].Set(-res[i].r, -res[i].i);
else
for (size_t i=1, j=n-1; 2*i<n; ++i, --j)
res[j].Set(res[i].r, -res[i].i);
res[half].Set(T(-1),T(0));
}
POCKETFFT_NOINLINE void sincos_2pibyn_half(size_t n, cmplx<T> * POCKETFFT_RESTRICT res)
......@@ -552,14 +544,18 @@ template<typename T> class sincos_2pibyn
}
public:
POCKETFFT_NOINLINE sincos_2pibyn(size_t n, bool half)
: data(n)
POCKETFFT_NOINLINE sincos_2pibyn(size_t n)
: sz(n), data(n/2+1)
{ sincos_2pibyn_half(n, data.data()); }
cmplx<T> operator[](size_t idx) const
{
sincos_2pibyn_half(n, data.data());
if (!half) fill_second_half(n, data.data());
if (idx<data.size())
return data[idx];
auto c = data[sz-idx];
c.i = -c.i;
return c;
}
cmplx<T> operator[](size_t idx) const { return data[idx]; }
};
struct util // hack to avoid duplicate symbols
......@@ -1613,7 +1609,7 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
void comp_twiddle()
{
sincos_2pibyn<T0> twiddle(length, false);
sincos_2pibyn<T0> twiddle(length);
size_t l1=1;
size_t memofs=0;
for (size_t k=0; k<fact.size(); ++k)
......@@ -2421,7 +2417,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
void comp_twiddle()
{
sincos_2pibyn<T0> twid(length, true);
sincos_2pibyn<T0> twid(length);
size_t l1=1;
T0 *ptr=mem.data();
for (size_t k=0; k<fact.size(); ++k)
......@@ -2515,7 +2511,7 @@ template<typename T0> class fftblue
bk(mem.data()), bkf(mem.data()+n)
{
/* initialize b_k */
sincos_2pibyn<T0> tmp(2*n, false);
sincos_2pibyn<T0> tmp(2*n);
bk[0].Set(1, 0);
size_t coeff=0;
......
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