Skip to content

GitLab

Commit e91402f6 authored by Martin Reinecke's avatar Martin Reinecke
Browse files
Options

bwd -> fwd

parent 6f26ecd1
Hide whitespace changes
Inline Side-by-side
Showing with 91 additions and 91 deletions
pocketfft_hdronly.h
View file @ e91402f6
......@@ -199,12 +199,12 @@ template<typename T> struct cmplx {
template<typename T2> auto operator* (const cmplx<T2> &other) const
-> cmplx<decltype(r+other.r)>
{ return {r*other.r-i*other.i, r*other.i + i*other.r}; }
template<bool bwd, typename T2> auto special_mul (const cmplx<T2> &other) const
template<bool fwd, typename T2> auto special_mul (const cmplx<T2> &other) const
-> cmplx<decltype(r+other.r)>
{
using Tres = cmplx<decltype(r+other.r)>;
return bwd ? Tres(r*other.r-i*other.i, r*other.i + i*other.r)
: Tres(r*other.r+i*other.i, i*other.r - r*other.i);
return fwd ? Tres(r*other.r+i*other.i, i*other.r-r*other.i)
: Tres(r*other.r-i*other.i, r*other.i+i*other.r);
}
};
template<typename T> void PMC(cmplx<T> &a, cmplx<T> &b,
......@@ -215,8 +215,8 @@ template<typename T> cmplx<T> conj(const cmplx<T> &a)
template<typename T> void ROT90(cmplx<T> &a)
{ auto tmp_=a.r; a.r=-a.i; a.i=tmp_; }
template<bool bwd, typename T> void ROTX90(cmplx<T> &a)
{ auto tmp_= bwd ? a.r : -a.r; a.r = bwd ? -a.i : a.i; a.i=tmp_; }
template<bool fwd, typename T> void ROTX90(cmplx<T> &a)
{ auto tmp_= fwd ? -a.r : a.r; a.r = fwd ? a.i : -a.i; a.i=tmp_; }
//
// twiddle factor section
......@@ -522,7 +522,7 @@ template<typename T0> class cfftp
void add_factor(size_t factor)
{ fact.push_back({factor, nullptr, nullptr}); }
template<bool bwd, typename T> void pass2 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=2;
......@@ -548,7 +548,7 @@ template<bool bwd, typename T> void pass2 (size_t ido, size_t l1,
for (size_t i=1; i<ido; ++i)
{
CH(i,k,0) = CC(i,0,k)+CC(i,1,k);
CH(i,k,1) = (CC(i,0,k)-CC(i,1,k)).template special_mul<bwd>(WA(0,i));
CH(i,k,1) = (CC(i,0,k)-CC(i,1,k)).template special_mul<fwd>(WA(0,i));
}
}
}
......@@ -570,15 +570,15 @@ template<bool bwd, typename T> void pass2 (size_t ido, size_t l1,
ca=t0+t1*twr; \
cb=t2*twi; ROT90(cb); \
PMC(da,db,ca,cb); \
CH(i,k,u1) = da.template special_mul<bwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<bwd>(WA(u2-1,i)); \
CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
}
template<bool bwd, typename T> void pass3 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=3;
constexpr T0 tw1r=-0.5,
tw1i= (bwd ? 1: -1) * T0(0.8660254037844386467637231707529362L);
tw1i= (fwd ? -1: 1) * T0(0.8660254037844386467637231707529362L);
auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T&
{ return ch[a+ido*(b+l1*c)]; };
......@@ -612,7 +612,7 @@ template<bool bwd, typename T> void pass3 (size_t ido, size_t l1,
#undef POCKETFFT_PARTSTEP3a
#undef POCKETFFT_PREP3
template<bool bwd, typename T> void pass4 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=4;
......@@ -630,7 +630,7 @@ template<bool bwd, typename T> void pass4 (size_t ido, size_t l1,
T t1, t2, t3, t4;
PMC(t2,t1,CC(0,0,k),CC(0,2,k));
PMC(t3,t4,CC(0,1,k),CC(0,3,k));
ROTX90<bwd>(t4);
ROTX90<fwd>(t4);
PMC(CH(0,k,0),CH(0,k,2),t2,t3);
PMC(CH(0,k,1),CH(0,k,3),t1,t4);
}
......@@ -641,7 +641,7 @@ template<bool bwd, typename T> void pass4 (size_t ido, size_t l1,
T t1, t2, t3, t4;
PMC(t2,t1,CC(0,0,k),CC(0,2,k));
PMC(t3,t4,CC(0,1,k),CC(0,3,k));
ROTX90<bwd>(t4);
ROTX90<fwd>(t4);
PMC(CH(0,k,0),CH(0,k,2),t2,t3);
PMC(CH(0,k,1),CH(0,k,3),t1,t4);
}
......@@ -651,12 +651,12 @@ template<bool bwd, typename T> void pass4 (size_t ido, size_t l1,
T cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k);
PMC(t2,t1,cc0,cc2);
PMC(t3,t4,cc1,cc3);
ROTX90<bwd>(t4);
ROTX90<fwd>(t4);
PMC(CH(i,k,0),c3,t2,t3);
PMC(c2,c4,t1,t4);
CH(i,k,1) = c2.template special_mul<bwd>(WA(0,i));
CH(i,k,2) = c3.template special_mul<bwd>(WA(1,i));
CH(i,k,3) = c4.template special_mul<bwd>(WA(2,i));
CH(i,k,1) = c2.template special_mul<fwd>(WA(0,i));
CH(i,k,2) = c3.template special_mul<fwd>(WA(1,i));
CH(i,k,3) = c4.template special_mul<fwd>(WA(2,i));
}
}
}
......@@ -686,17 +686,17 @@ template<bool bwd, typename T> void pass4 (size_t ido, size_t l1,
cb.i=twai*t4.r twbi*t3.r; \
cb.r=-(twai*t4.i twbi*t3.i); \
PMC(da,db,ca,cb); \
CH(i,k,u1) = da.template special_mul<bwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<bwd>(WA(u2-1,i)); \
CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
}
template<bool bwd, typename T> void pass5 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=5;
constexpr T0 tw1r= T0(0.3090169943749474241022934171828191L),
tw1i= (bwd ? 1: -1) * T0(0.9510565162951535721164393333793821L),
tw1i= (fwd ? -1: 1) * T0(0.9510565162951535721164393333793821L),
tw2r= T0(-0.8090169943749474241022934171828191L),
tw2i= (bwd ? 1: -1) * T0(0.5877852522924731291687059546390728L);
tw2i= (fwd ? -1: 1) * T0(0.5877852522924731291687059546390728L);
auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T&
{ return ch[a+ido*(b+l1*c)]; };
......@@ -756,20 +756,20 @@ template<bool bwd, typename T> void pass5 (size_t ido, size_t l1,
{ \
T da,db; \
POCKETFFT_PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,da,db) \
CH(i,k,u1) = da.template special_mul<bwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<bwd>(WA(u2-1,i)); \
CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
}
template<bool bwd, typename T> void pass7(size_t ido, size_t l1,
template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=7;
constexpr T0 tw1r= T0(0.6234898018587335305250048840042398L),
tw1i= (bwd ? 1 : -1) * T0(0.7818314824680298087084445266740578L),
tw1i= (fwd ? -1 : 1) * T0(0.7818314824680298087084445266740578L),
tw2r= T0(-0.2225209339563144042889025644967948L),
tw2i= (bwd ? 1 : -1) * T0(0.9749279121818236070181316829939312L),
tw2i= (fwd ? -1 : 1) * T0(0.9749279121818236070181316829939312L),
tw3r= T0(-0.9009688679024191262361023195074451L),
tw3i= (bwd ? 1 : -1) * T0(0.433883739117558120475768332848359L);
tw3i= (fwd ? -1 : 1) * T0(0.433883739117558120475768332848359L);
auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T&
{ return ch[a+ido*(b+l1*c)]; };
......@@ -810,26 +810,26 @@ template<bool bwd, typename T> void pass7(size_t ido, size_t l1,
#undef POCKETFFT_PARTSTEP7a
#undef POCKETFFT_PREP7
template <bool bwd, typename T> void ROTX45(T &a)
template <bool fwd, typename T> void ROTX45(T &a)
{
constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
if (bwd)
{ auto tmp_=a.r; a.r=hsqt2*(a.r-a.i); a.i=hsqt2*(a.i+tmp_); }
else
if (fwd)
{ auto tmp_=a.r; a.r=hsqt2*(a.r+a.i); a.i=hsqt2*(a.i-tmp_); }
else
{ auto tmp_=a.r; a.r=hsqt2*(a.r-a.i); a.i=hsqt2*(a.i+tmp_); }
}
template <bool bwd, typename T> void ROTX135(T &a)
template <bool fwd, typename T> void ROTX135(T &a)
{
constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
if (bwd)
{ auto tmp_=a.r; a.r=hsqt2*(-a.r-a.i); a.i=hsqt2*(tmp_-a.i); }
else
if (fwd)
{ auto tmp_=a.r; a.r=hsqt2*(a.i-a.r); a.i=hsqt2*(-tmp_-a.i); }
else
{ auto tmp_=a.r; a.r=hsqt2*(-a.r-a.i); a.i=hsqt2*(tmp_-a.i); }
}
template<typename T> inline void PMINPLACE(T &a, T &b)
{ T t = a; a.r+=b.r; a.i+=b.i; b.r=t.r-b.r; b.i=t.i-b.i; }
template<bool bwd, typename T> void pass8 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=8;
......@@ -849,15 +849,15 @@ template<bool bwd, typename T> void pass8 (size_t ido, size_t l1,
PMC(a1,a5,CC(0,1,k),CC(0,5,k));
PMC(a2,a6,CC(0,2,k),CC(0,6,k));
PMC(a3,a7,CC(0,3,k),CC(0,7,k));
ROTX90<bwd>(a6);
ROTX90<bwd>(a7);
ROTX90<fwd>(a6);
ROTX90<fwd>(a7);
PMINPLACE(a0,a2);
PMINPLACE(a1,a3);
PMINPLACE(a4,a6);
PMINPLACE(a5,a7);
ROTX45<bwd>(a5);
ROTX90<bwd>(a3);
ROTX135<bwd>(a7);
ROTX45<fwd>(a5);
ROTX90<fwd>(a3);
ROTX135<fwd>(a7);
PMC(CH(0,k,0),CH(0,k,4),a0,a1);
PMC(CH(0,k,1),CH(0,k,5),a4,a5);
PMC(CH(0,k,2),CH(0,k,6),a2,a3);
......@@ -871,15 +871,15 @@ template<bool bwd, typename T> void pass8 (size_t ido, size_t l1,
PMC(a1,a5,CC(0,1,k),CC(0,5,k));
PMC(a2,a6,CC(0,2,k),CC(0,6,k));
PMC(a3,a7,CC(0,3,k),CC(0,7,k));
ROTX90<bwd>(a6);
ROTX90<bwd>(a7);
ROTX90<fwd>(a6);
ROTX90<fwd>(a7);
PMINPLACE(a0,a2);
PMINPLACE(a1,a3);
PMINPLACE(a4,a6);
PMINPLACE(a5,a7);
ROTX45<bwd>(a5);
ROTX90<bwd>(a3);
ROTX135<bwd>(a7);
ROTX45<fwd>(a5);
ROTX90<fwd>(a3);
ROTX135<fwd>(a7);
PMC(CH(0,k,0),CH(0,k,4),a0,a1);
PMC(CH(0,k,1),CH(0,k,5),a4,a5);
PMC(CH(0,k,2),CH(0,k,6),a2,a3);
......@@ -892,27 +892,27 @@ template<bool bwd, typename T> void pass8 (size_t ido, size_t l1,
PMC(a1,a5,CC(i,1,k),CC(i,5,k));
PMC(a2,a6,CC(i,2,k),CC(i,6,k));
PMC(a3,a7,CC(i,3,k),CC(i,7,k));
ROTX90<bwd>(a6);
ROTX90<bwd>(a7);
ROTX90<fwd>(a6);
ROTX90<fwd>(a7);
PMINPLACE(a0,a2);
PMINPLACE(a1,a3);
PMINPLACE(a4,a6);
PMINPLACE(a5,a7);
ROTX45<bwd>(a5);
ROTX90<bwd>(a3);
ROTX135<bwd>(a7);
ROTX45<fwd>(a5);
ROTX90<fwd>(a3);
ROTX135<fwd>(a7);
PMINPLACE(a0,a1);
PMINPLACE(a2,a3);
PMINPLACE(a4,a5);
PMINPLACE(a6,a7);
CH(i,k,0) = a0;
CH(i,k,1) = a4.template special_mul<bwd>(WA(0,i));
CH(i,k,2) = a2.template special_mul<bwd>(WA(1,i));
CH(i,k,3) = a6.template special_mul<bwd>(WA(2,i));
CH(i,k,4) = a1.template special_mul<bwd>(WA(3,i));
CH(i,k,5) = a5.template special_mul<bwd>(WA(4,i));
CH(i,k,6) = a3.template special_mul<bwd>(WA(5,i));
CH(i,k,7) = a7.template special_mul<bwd>(WA(6,i));
CH(i,k,1) = a4.template special_mul<fwd>(WA(0,i));
CH(i,k,2) = a2.template special_mul<fwd>(WA(1,i));
CH(i,k,3) = a6.template special_mul<fwd>(WA(2,i));
CH(i,k,4) = a1.template special_mul<fwd>(WA(3,i));
CH(i,k,5) = a5.template special_mul<fwd>(WA(4,i));
CH(i,k,6) = a3.template special_mul<fwd>(WA(5,i));
CH(i,k,7) = a7.template special_mul<fwd>(WA(6,i));
}
}
}
......@@ -942,24 +942,24 @@ template<bool bwd, typename T> void pass8 (size_t ido, size_t l1,
{ \
T da,db; \
PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,da,db) \
CH(i,k,u1) = da.template special_mul<bwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<bwd>(WA(u2-1,i)); \
CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
}
template<bool bwd, typename T> void pass11 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
const T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa)
{
constexpr size_t cdim=11;
constexpr T0 tw1r= T0(0.8412535328311811688618116489193677L),
tw1i= (bwd ? 1 : -1) * T0(0.5406408174555975821076359543186917L),
tw1i= (fwd ? -1 : 1) * T0(0.5406408174555975821076359543186917L),
tw2r= T0(0.4154150130018864255292741492296232L),
tw2i= (bwd ? 1 : -1) * T0(0.9096319953545183714117153830790285L),
tw2i= (fwd ? -1 : 1) * T0(0.9096319953545183714117153830790285L),
tw3r= T0(-0.1423148382732851404437926686163697L),
tw3i= (bwd ? 1 : -1) * T0(0.9898214418809327323760920377767188L),
tw3i= (fwd ? -1 : 1) * T0(0.9898214418809327323760920377767188L),
tw4r= T0(-0.6548607339452850640569250724662936L),
tw4i= (bwd ? 1 : -1) * T0(0.7557495743542582837740358439723444L),
tw4i= (fwd ? -1 : 1) * T0(0.7557495743542582837740358439723444L),
tw5r= T0(-0.9594929736144973898903680570663277L),
tw5i= (bwd ? 1 : -1) * T0(0.2817325568414296977114179153466169L);
tw5i= (fwd ? -1 : 1) * T0(0.2817325568414296977114179153466169L);
auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T&
{ return ch[a+ido*(b+l1*c)]; };
......@@ -1006,7 +1006,7 @@ template<bool bwd, typename T> void pass11 (size_t ido, size_t l1,
#undef PARTSTEP11a
#undef POCKETFFT_PREP11
template<bool bwd, typename T> void passg (size_t ido, size_t ip,
template<bool fwd, typename T> void passg (size_t ido, size_t ip,
size_t l1, T * RESTRICT cc, T * RESTRICT ch, const cmplx<T0> * RESTRICT wa,
const cmplx<T0> * RESTRICT csarr)
{
......@@ -1028,7 +1028,7 @@ template<bool bwd, typename T> void passg (size_t ido, size_t ip,
arr<cmplx<T0>> wal(ip);
wal[0] = cmplx<T0>(1., 0.);
for (size_t i=1; i<ip; ++i)
wal[i]=cmplx<T0>(csarr[i].r,bwd ? csarr[i].i : -csarr[i].i);
wal[i]=cmplx<T0>(csarr[i].r,fwd ? -csarr[i].i : csarr[i].i);
for (size_t k=0; k<l1; ++k)
for (size_t i=0; i<ido; ++i)
......@@ -1106,15 +1106,15 @@ template<bool bwd, typename T> void passg (size_t ido, size_t ip,
T x1, x2;
PMC(x1,x2,CX(i,k,j),CX(i,k,jc));
size_t idij=(j-1)*(ido-1)+i-1;
CX(i,k,j) = x1.template special_mul<bwd>(wa[idij]);
CX(i,k,j) = x1.template special_mul<fwd>(wa[idij]);
idij=(jc-1)*(ido-1)+i-1;
CX(i,k,jc) = x2.template special_mul<bwd>(wa[idij]);
CX(i,k,jc) = x2.template special_mul<fwd>(wa[idij]);
}
}
}
}
template<bool bwd, typename T> void pass_all(T c[], T0 fct)
template<bool fwd, typename T> void pass_all(T c[], T0 fct)
{
if (length==1) { c[0]*=fct; return; }
size_t l1=1;
......@@ -1127,22 +1127,22 @@ template<bool bwd, typename T> void pass_all(T c[], T0 fct)
size_t l2=ip*l1;
size_t ido = length/l2;
if (ip==4)
pass4<bwd> (ido, l1, p1, p2, fact[k1].tw);
pass4<fwd> (ido, l1, p1, p2, fact[k1].tw);
else if(ip==8)
pass8<bwd>(ido, l1, p1, p2, fact[k1].tw);
pass8<fwd>(ido, l1, p1, p2, fact[k1].tw);
else if(ip==2)
pass2<bwd>(ido, l1, p1, p2, fact[k1].tw);
pass2<fwd>(ido, l1, p1, p2, fact[k1].tw);
else if(ip==3)
pass3<bwd> (ido, l1, p1, p2, fact[k1].tw);
pass3<fwd> (ido, l1, p1, p2, fact[k1].tw);
else if(ip==5)
pass5<bwd> (ido, l1, p1, p2, fact[k1].tw);
pass5<fwd> (ido, l1, p1, p2, fact[k1].tw);
else if(ip==7)
pass7<bwd> (ido, l1, p1, p2, fact[k1].tw);
pass7<fwd> (ido, l1, p1, p2, fact[k1].tw);
else if(ip==11)
pass11<bwd> (ido, l1, p1, p2, fact[k1].tw);
pass11<fwd> (ido, l1, p1, p2, fact[k1].tw);
else
{
passg<bwd>(ido, ip, l1, p1, p2, fact[k1].tw, fact[k1].tws);
passg<fwd>(ido, ip, l1, p1, p2, fact[k1].tw, fact[k1].tws);
swap(p1,p2);
}
swap(p1,p2);
......@@ -1164,10 +1164,10 @@ template<bool bwd, typename T> void pass_all(T c[], T0 fct)
public:
template<typename T> void forward(T c[], T0 fct)
{ pass_all<false>(c, fct); }
{ pass_all<true>(c, fct); }
template<typename T> void backward(T c[], T0 fct)
{ pass_all<true>(c, fct); }
{ pass_all<false>(c, fct); }
private:
NOINLINE void factorize()
......@@ -2088,13 +2088,13 @@ template<typename T0> class fftblue
arr<cmplx<T0>> mem;
cmplx<T0> *bk, *bkf;
template<bool bwd, typename T> void fft(cmplx<T> c[], T0 fct)
template<bool fwd, typename T> void fft(cmplx<T> c[], T0 fct)
{
arr<cmplx<T>> akf(n2);
/* initialize a_k and FFT it */
for (size_t m=0; m<n; ++m)
akf[m] = c[m].template special_mul<bwd>(bk[m]);
akf[m] = c[m].template special_mul<fwd>(bk[m]);
auto zero = akf[0]*T0(0);
for (size_t m=n; m<n2; ++m)
akf[m]=zero;
......@@ -2103,14 +2103,14 @@ template<typename T0> class fftblue
/* do the convolution */
for (size_t m=0; m<n2; ++m)
akf[m] = akf[m].template special_mul<!bwd>(bkf[m]);
akf[m] = akf[m].template special_mul<!fwd>(bkf[m]);
/* inverse FFT */
plan.backward (akf.data(),1.);
/* multiply by b_k */
for (size_t m=0; m<n; ++m)
c[m] = akf[m].template special_mul<bwd>(bk[m])*fct;
c[m] = akf[m].template special_mul<fwd>(bk[m])*fct;
}
public:
......@@ -2142,10 +2142,10 @@ template<typename T0> class fftblue
}
template<typename T> void backward(cmplx<T> c[], T0 fct)
{ fft<true>(c,fct); }
{ fft<false>(c,fct); }
template<typename T> void forward(cmplx<T> c[], T0 fct)
{ fft<false>(c,fct); }
{ fft<true>(c,fct); }
template<typename T> void backward_r(T c[], T0 fct)
{
......@@ -2156,7 +2156,7 @@ template<typename T0> class fftblue
if ((n&1)==0) tmp[n/2].i=T0(0)*c[0];
for (size_t m=1; 2*m<n; ++m)
tmp[n-m].Set(tmp[m].r, -tmp[m].i);
fft<true>(tmp.data(),fct);
fft<false>(tmp.data(),fct);
for (size_t m=0; m<n; ++m)
c[m] = tmp[m].r;
}
......@@ -2167,7 +2167,7 @@ template<typename T0> class fftblue
auto zero = T0(0)*c[0];
for (size_t m=0; m<n; ++m)
tmp[m].Set(c[m], zero);
fft<false>(tmp.data(),fct);
fft<true>(tmp.data(),fct);
c[0] = tmp[0].r;
memcpy (c+1, tmp.data()+1, (n-1)*sizeof(T));
}
......
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