refactored ei_kissfft_impl to maintain a cache of cpx fft plans

2 files changed, 317 insertions, 296 deletions

diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
index ce2c9f16e..3580e6c61 100644
--- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
+++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
@@ -24,21 +24,279 @@
 
 #include <complex>
 #include <vector>
+#include <map>
 
 namespace Eigen {
 
+    template <typename _Scalar>
+    struct ei_kiss_cpx_fft
+    {
+        typedef  _Scalar Scalar;
+        typedef  std::complex<Scalar> Complex;
+        std::vector<Complex> m_twiddles;
+        std::vector<int> m_stageRadix;
+        std::vector<int> m_stageRemainder;
+        bool m_inverse;
+
+        ei_kiss_cpx_fft() { }
+
+        void make_twiddles(int nfft,bool inverse)
+        {
+            m_inverse = inverse;
+            m_twiddles.resize(nfft);
+            Scalar phinc =  (inverse?2:-2)* acos( (Scalar) -1)  / nfft;
+            for (int i=0;i<nfft;++i)
+                m_twiddles[i] = exp( Complex(0,i*phinc) );
+        }
+
+        void invert()
+        {
+            m_inverse = !m_inverse;
+            for ( size_t i=0;i<m_twiddles.size() ;++i)
+                m_twiddles[i] = conj( m_twiddles[i] );
+        }
+
+        void factorize(int nfft)
+        {
+            if (m_stageRadix.size()==0 || m_stageRadix[0] * m_stageRemainder[0] != nfft)
+            {
+                m_stageRadix.resize(0);
+                m_stageRemainder.resize(0);
+                //factorize
+                //start factoring out 4's, then 2's, then 3,5,7,9,...
+                int n= nfft;
+                int p=4;
+                do {
+                    while (n % p) {
+                        switch (p) {
+                            case 4: p = 2; break;
+                            case 2: p = 3; break;
+                            default: p += 2; break;
+                        }
+                        if (p*p>n)
+                            p=n;// impossible to have a factor > sqrt(n)
+                    }
+                    n /= p;
+                    m_stageRadix.push_back(p);
+                    m_stageRemainder.push_back(n);
+                }while(n>1);
+            }
+        }
+
+        template <typename _Src>
+            void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride)
+            {
+                int p = m_stageRadix[stage];
+                int m = m_stageRemainder[stage];
+                Complex * Fout_beg = xout;
+                Complex * Fout_end = xout + p*m;
+
+                if (m>1) {
+                    do{
+                        // recursive call:
+                        // DFT of size m*p performed by doing
+                        // p instances of smaller DFTs of size m, 
+                        // each one takes a decimated version of the input
+                        work(stage+1, xout , xin, fstride*p,in_stride);
+                        xin += fstride*in_stride;
+                    }while( (xout += m) != Fout_end );
+                }else{
+                    do{
+                        *xout = *xin;
+                        xin += fstride*in_stride;
+                    }while(++xout != Fout_end );
+                }
+                xout=Fout_beg;
+
+                // recombine the p smaller DFTs 
+                switch (p) {
+                    case 2: bfly2(xout,fstride,m); break;
+                    case 3: bfly3(xout,fstride,m); break;
+                    case 4: bfly4(xout,fstride,m); break;
+                    case 5: bfly5(xout,fstride,m); break;
+                    default: bfly_generic(xout,fstride,m,p); break;
+                }
+            }
+
+        void bfly2( Complex * Fout, const size_t fstride, int m)
+        {
+            for (int k=0;k<m;++k) {
+                Complex t = Fout[m+k] * m_twiddles[k*fstride];
+                Fout[m+k] = Fout[k] - t;
+                Fout[k] += t;
+            }
+        }
+
+        void bfly4( Complex * Fout, const size_t fstride, const size_t m)
+        {
+            Complex scratch[6];
+            int negative_if_inverse = m_inverse * -2 +1;
+            for (size_t k=0;k<m;++k) {
+                scratch[0] = Fout[k+m] * m_twiddles[k*fstride];
+                scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2];
+                scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3];
+                scratch[5] = Fout[k] - scratch[1];
+
+                Fout[k] += scratch[1];
+                scratch[3] = scratch[0] + scratch[2];
+                scratch[4] = scratch[0] - scratch[2];
+                scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse );
+
+                Fout[k+2*m]  = Fout[k] - scratch[3];
+                Fout[k] += scratch[3];
+                Fout[k+m] = scratch[5] + scratch[4];
+                Fout[k+3*m] = scratch[5] - scratch[4];
+            }
+        }
+
+        void bfly3( Complex * Fout, const size_t fstride, const size_t m)
+        {
+            size_t k=m;
+            const size_t m2 = 2*m;
+            Complex *tw1,*tw2;
+            Complex scratch[5];
+            Complex epi3;
+            epi3 = m_twiddles[fstride*m];
+
+            tw1=tw2=&m_twiddles[0];
+
+            do{
+                scratch[1]=Fout[m] * *tw1;
+                scratch[2]=Fout[m2] * *tw2;
+
+                scratch[3]=scratch[1]+scratch[2];
+                scratch[0]=scratch[1]-scratch[2];
+                tw1 += fstride;
+                tw2 += fstride*2;
+                Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() );
+                scratch[0] *= epi3.imag();
+                *Fout += scratch[3];
+                Fout[m2] = Complex(  Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );
+                Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() );
+                ++Fout;
+            }while(--k);
+        }
+
+        void bfly5( Complex * Fout, const size_t fstride, const size_t m)
+        {
+            Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
+            size_t u;
+            Complex scratch[13];
+            Complex * twiddles = &m_twiddles[0];
+            Complex *tw;
+            Complex ya,yb;
+            ya = twiddles[fstride*m];
+            yb = twiddles[fstride*2*m];
+
+            Fout0=Fout;
+            Fout1=Fout0+m;
+            Fout2=Fout0+2*m;
+            Fout3=Fout0+3*m;
+            Fout4=Fout0+4*m;
+
+            tw=twiddles;
+            for ( u=0; u<m; ++u ) {
+                scratch[0] = *Fout0;
+
+                scratch[1]  = *Fout1 * tw[u*fstride];
+                scratch[2]  = *Fout2 * tw[2*u*fstride];
+                scratch[3]  = *Fout3 * tw[3*u*fstride];
+                scratch[4]  = *Fout4 * tw[4*u*fstride];
+
+                scratch[7] = scratch[1] + scratch[4];
+                scratch[10] = scratch[1] - scratch[4];
+                scratch[8] = scratch[2] + scratch[3];
+                scratch[9] = scratch[2] - scratch[3];
+
+                *Fout0 +=  scratch[7];
+                *Fout0 +=  scratch[8];
+
+                scratch[5] = scratch[0] + Complex(
+                        (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ),
+                        (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real())
+                        );
+
+                scratch[6] = Complex(
+                        (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()),
+                        -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag())
+                        );
+
+                *Fout1 = scratch[5] - scratch[6];
+                *Fout4 = scratch[5] + scratch[6];
+
+                scratch[11] = scratch[0] +
+                    Complex(
+                            (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()),
+                            (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real())
+                           );
+
+                scratch[12] = Complex(
+                        -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()),
+                        (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag())
+                        );
+
+                *Fout2=scratch[11]+scratch[12];
+                *Fout3=scratch[11]-scratch[12];
+
+                ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
+            }
+        }
+
+        /* perform the butterfly for one stage of a mixed radix FFT */
+        void bfly_generic(
+                Complex * Fout,
+                const size_t fstride,
+                int m,
+                int p
+                )
+        {
+            int u,k,q1,q;
+            Complex * twiddles = &m_twiddles[0];
+            Complex t;
+            int Norig = m_twiddles.size();
+            Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
+
+            for ( u=0; u<m; ++u ) {
+                k=u;
+                for ( q1=0 ; q1<p ; ++q1 ) {
+                    scratchbuf[q1] = Fout[ k  ];
+                    k += m;
+                }
+
+                k=u;
+                for ( q1=0 ; q1<p ; ++q1 ) {
+                    int twidx=0;
+                    Fout[ k ] = scratchbuf[0];
+                    for (q=1;q<p;++q ) {
+                        twidx += fstride * k;
+                        if (twidx>=Norig) twidx-=Norig;
+                        t=scratchbuf[q] * twiddles[twidx];
+                        Fout[ k ] += t;
+                    }
+                    k += m;
+                }
+            }
+        }
+    };
+
+
   template <typename _Scalar>
   struct ei_kissfft_impl
   {
     typedef _Scalar Scalar;
     typedef std::complex<Scalar> Complex;
-    ei_kissfft_impl() : m_nfft(0) {} 
+    ei_kissfft_impl() {} 
+
+    void clear() 
+    {
+        m_plans.clear();
+        m_realTwiddles.clear();
+    }
 
     template <typename _Src>
     void fwd( Complex * dst,const _Src *src,int nfft)
     {
-        prepare(nfft,false);
-        work(0, dst, src, 1,1);
+        get_plan(nfft,false).work(0, dst, src, 1,1);
     }
 
     // real-to-complex forward FFT
@@ -47,16 +305,16 @@ namespace Eigen {
     // then fill in the conjugate symmetric half
     void fwd( Complex * dst,const Scalar * src,int nfft) 
     {
-        if ( nfft&1 ) {
+        if ( nfft&3  ) {
             // use generic mode for odd
-            prepare(nfft,false);
-            work(0, dst, src, 1,1);
+            get_plan(nfft,false).work(0, dst, src, 1,1);
         }else{
             int ncfft = nfft>>1;
             int ncfft2 = nfft>>2;
+            Complex * rtw = real_twiddles(ncfft2);
+
             // use optimized mode for even real
-            fwd( dst, reinterpret_cast<const Complex*> (src),ncfft);
-            make_real_twiddles(nfft);
+            fwd( dst, reinterpret_cast<const Complex*> (src), ncfft);
             Complex dc = dst[0].real() +  dst[0].imag();
             Complex nyquist = dst[0].real() -  dst[0].imag();
             int k;
@@ -65,8 +323,7 @@ namespace Eigen {
                 Complex fpnk = conj(dst[ncfft-k]);
                 Complex f1k = fpk + fpnk;
                 Complex f2k = fpk - fpnk;
-                //Complex tw = f2k * exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
-                Complex tw= f2k * m_realTwiddles[k-1];
+                Complex tw= f2k * rtw[k-1];
 
                 dst[k] =  (f1k + tw) * Scalar(.5);
                 dst[ncfft-k] =  conj(f1k -tw)*Scalar(.5);
@@ -94,304 +351,67 @@ namespace Eigen {
 
     void inv(Complex * dst,const Complex  *src,int nfft)
     {
-        prepare(nfft,true);
-        work(0, dst, src, 1,1);
-        scale(dst, Scalar(1)/m_nfft );
-    }
-
-    void prepare(int nfft,bool inverse)
-    {
-        make_twiddles(nfft,inverse);
-        factorize(nfft);
-    }
-
-    void make_real_twiddles(int nfft)
-    {
-        int ncfft2 = nfft>>2;
-        if ( m_realTwiddles.size() != ncfft2) {
-            m_realTwiddles.resize(ncfft2);
-            int ncfft= nfft>>1;
-            for (int k=1;k<=ncfft2;++k) 
-                m_realTwiddles[k-1] = exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
-        }
-    }
-
-    void make_twiddles(int nfft,bool inverse)
-    {
-        if ( m_twiddles.size() == nfft) {
-            // reuse the twiddles, conjugate if necessary
-            if (inverse != m_inverse)
-                for (int i=0;i<nfft;++i)
-                    m_twiddles[i] = conj( m_twiddles[i] );
-        }else{
-            m_twiddles.resize(nfft);
-            Scalar phinc =  (inverse?2:-2)* acos( (Scalar) -1)  / nfft;
-            for (int i=0;i<nfft;++i)
-                m_twiddles[i] = exp( Complex(0,i*phinc) );
-        }
-        m_inverse = inverse;
-    }
-
-    void factorize(int nfft)
-    {
-        if (m_stageRadix.size()==0 || m_stageRadix[0] * m_stageRemainder[0] != nfft)
-        {
-            m_stageRadix.resize(0);
-            m_stageRemainder.resize(0);
-            //factorize
-            //start factoring out 4's, then 2's, then 3,5,7,9,...
-            int n= nfft;
-            int p=4;
-            do {
-                while (n % p) {
-                    switch (p) {
-                        case 4: p = 2; break;
-                        case 2: p = 3; break;
-                        default: p += 2; break;
-                    }
-                    if (p*p>n)
-                        p=n;// impossible to have a factor > sqrt(n)
-                }
-                n /= p;
-                m_stageRadix.push_back(p);
-                m_stageRemainder.push_back(n);
-            }while(n>1);
-        }
-        m_nfft = nfft;
-    }
-
-    void scale(Complex *dst,Scalar s) 
-    {
-        for (int k=0;k<m_nfft;++k)
-            dst[k] *= s;
+        get_plan(nfft,true).work(0, dst, src, 1,1);
+        scale(dst, nfft, Scalar(1)/nfft );
     }
 
     private:
 
-    template <typename _Src>
-    void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride)
-    {
-      int p = m_stageRadix[stage];
-      int m = m_stageRemainder[stage];
-      Complex * Fout_beg = xout;
-      Complex * Fout_end = xout + p*m;
-
-      if (m>1) {
-        do{
-          // recursive call:
-          // DFT of size m*p performed by doing
-          // p instances of smaller DFTs of size m, 
-          // each one takes a decimated version of the input
-          work(stage+1, xout , xin, fstride*p,in_stride);
-          xin += fstride*in_stride;
-        }while( (xout += m) != Fout_end );
-      }else{
-          do{
-              *xout = *xin;
-              xin += fstride*in_stride;
-          }while(++xout != Fout_end );
-      }
-      xout=Fout_beg;
-
-      // recombine the p smaller DFTs 
-      switch (p) {
-        case 2: bfly2(xout,fstride,m); break;
-        case 3: bfly3(xout,fstride,m); break;
-        case 4: bfly4(xout,fstride,m); break;
-        case 5: bfly5(xout,fstride,m); break;
-        default: bfly_generic(xout,fstride,m,p); break;
-      }
-    }
-
-    void bfly2( Complex * Fout, const size_t fstride, int m)
-    {
-      for (int k=0;k<m;++k) {
-        Complex t = Fout[m+k] * m_twiddles[k*fstride];
-        Fout[m+k] = Fout[k] - t;
-        Fout[k] += t;
-      }
-    }
-
-    void bfly4( Complex * Fout, const size_t fstride, const size_t m)
-    {
-      Complex scratch[6];
-      int negative_if_inverse = m_inverse * -2 +1;
-      for (size_t k=0;k<m;++k) {
-        scratch[0] = Fout[k+m] * m_twiddles[k*fstride];
-        scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2];
-        scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3];
-        scratch[5] = Fout[k] - scratch[1];
-
-        Fout[k] += scratch[1];
-        scratch[3] = scratch[0] + scratch[2];
-        scratch[4] = scratch[0] - scratch[2];
-        scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse );
-
-        Fout[k+2*m]  = Fout[k] - scratch[3];
-        Fout[k] += scratch[3];
-        Fout[k+m] = scratch[5] + scratch[4];
-        Fout[k+3*m] = scratch[5] - scratch[4];
-      }
-    }
+    typedef ei_kiss_cpx_fft<Scalar> PlanData;
 
-    void bfly3( Complex * Fout, const size_t fstride, const size_t m)
-    {
-      size_t k=m;
-      const size_t m2 = 2*m;
-      Complex *tw1,*tw2;
-      Complex scratch[5];
-      Complex epi3;
-      epi3 = m_twiddles[fstride*m];
-
-      tw1=tw2=&m_twiddles[0];
-
-      do{
-        scratch[1]=Fout[m] * *tw1;
-        scratch[2]=Fout[m2] * *tw2;
-
-        scratch[3]=scratch[1]+scratch[2];
-        scratch[0]=scratch[1]-scratch[2];
-        tw1 += fstride;
-        tw2 += fstride*2;
-        Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() );
-        scratch[0] *= epi3.imag();
-        *Fout += scratch[3];
-        Fout[m2] = Complex(  Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );
-        Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() );
-        ++Fout;
-      }while(--k);
-    }
+    typedef std::map<int,PlanData> PlanMap;
+    PlanMap m_plans;
+    std::map<int, std::vector<Complex> > m_realTwiddles;
 
-    void bfly5( Complex * Fout, const size_t fstride, const size_t m)
-    {
-      Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
-      size_t u;
-      Complex scratch[13];
-      Complex * twiddles = &m_twiddles[0];
-      Complex *tw;
-      Complex ya,yb;
-      ya = twiddles[fstride*m];
-      yb = twiddles[fstride*2*m];
-
-      Fout0=Fout;
-      Fout1=Fout0+m;
-      Fout2=Fout0+2*m;
-      Fout3=Fout0+3*m;
-      Fout4=Fout0+4*m;
-
-      tw=twiddles;
-      for ( u=0; u<m; ++u ) {
-        scratch[0] = *Fout0;
-
-        scratch[1]  = *Fout1 * tw[u*fstride];
-        scratch[2]  = *Fout2 * tw[2*u*fstride];
-        scratch[3]  = *Fout3 * tw[3*u*fstride];
-        scratch[4]  = *Fout4 * tw[4*u*fstride];
-
-        scratch[7] = scratch[1] + scratch[4];
-        scratch[10] = scratch[1] - scratch[4];
-        scratch[8] = scratch[2] + scratch[3];
-        scratch[9] = scratch[2] - scratch[3];
-
-        *Fout0 +=  scratch[7];
-        *Fout0 +=  scratch[8];
-
-        scratch[5] = scratch[0] + Complex(
-            (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ),
-            (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real())
-            );
-
-        scratch[6] = Complex(
-            (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()),
-            -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag())
-            );
-
-        *Fout1 = scratch[5] - scratch[6];
-        *Fout4 = scratch[5] + scratch[6];
-
-        scratch[11] = scratch[0] +
-          Complex(
-              (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()),
-              (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real())
-              );
-
-        scratch[12] = Complex(
-            -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()),
-            (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag())
-            );
-
-        *Fout2=scratch[11]+scratch[12];
-        *Fout3=scratch[11]-scratch[12];
-
-        ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
-      }
-    }
+    int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; }
 
-    /* perform the butterfly for one stage of a mixed radix FFT */
-    void bfly_generic(
-        Complex * Fout,
-        const size_t fstride,
-        int m,
-        int p
-        )
+    PlanData & get_plan(int nfft,bool inverse)
     {
-      int u,k,q1,q;
-      Complex * twiddles = &m_twiddles[0];
-      Complex t;
-      int Norig = m_nfft;
-      Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
-
-      for ( u=0; u<m; ++u ) {
-        k=u;
-        for ( q1=0 ; q1<p ; ++q1 ) {
-          scratchbuf[q1] = Fout[ k  ];
-          k += m;
+        /* 
+         * for some reason this does not work
+         *
+        typedef typename std::map<int,PlanData>::iterator MapIt;
+        MapIt it;
+        it = m_plans.find( PlanKey(nfft,inverse) );
+        if (it == m_plans.end() ) {
+            // create new entry
+            it = m_plans.insert( make_pair( PlanKey(nfft,inverse) , PlanData() ) );
+            MapIt it2 = m_plans.find( PlanKey(nfft,!inverse) );
+            if (it2 != m_plans.end() ) {
+                it->second = it2.second;
+                it->second.invert();
+            }else{
+                it->second.make_twiddles(nfft,inverse);
+                it->second.factorize(nfft);
+            }
         }
-
-        k=u;
-        for ( q1=0 ; q1<p ; ++q1 ) {
-          int twidx=0;
-          Fout[ k ] = scratchbuf[0];
-          for (q=1;q<p;++q ) {
-            twidx += fstride * k;
-            if (twidx>=Norig) twidx-=Norig;
-            t=scratchbuf[q] * twiddles[twidx];
-            Fout[ k ] += t;
-          }
-          k += m;
+        return it->second;
+        */
+        PlanData & pd = m_plans[ PlanKey(nfft,inverse) ];
+        if ( pd.m_twiddles.size() == 0 ) {
+            pd.make_twiddles(nfft,inverse);
+            pd.factorize(nfft);
         }
-      }
+        return pd;
     }
 
-    int m_nfft;
-    bool m_inverse;
-    std::vector<Complex> m_twiddles;
-    std::vector<Complex> m_realTwiddles;
-    std::vector<int> m_stageRadix;
-    std::vector<int> m_stageRemainder;
-/*
-    enum {FORWARD,INVERSE,REAL,COMPLEX};
-
-    struct PlanKey
+    Complex * real_twiddles(int ncfft2)
     {
-        PlanKey(int nfft,bool isinverse,bool iscomplex)
-        {
-            _key = (nfft<<2) | (isinverse<<1) | iscomplex;
-        }
-
-        bool operator<(const PlanKey & other) const
-        {
-            return this->_key < other._key;
+        std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there
+        if ( (int)twidref.size() != ncfft2 ) {
+            twidref.resize(ncfft2);
+            int ncfft= ncfft2<<1;
+            Scalar pi =  acos( Scalar(-1) );
+            for (int k=1;k<=ncfft2;++k) 
+                twidref[k-1] = exp( Complex(0,-pi * ((double) (k) / ncfft + .5) ) );
         }
-        int _key;
-    };
+        return &twidref[0];
+    }
 
-    struct PlanData
+    void scale(Complex *dst,int n,Scalar s) 
     {
-        std::vector<Complex> m_twiddles;
-    };
-
-    std::map<PlanKey,
-*/
+        for (int k=0;k<n;++k)
+            dst[k] *= s;
+    }
   };
 }
diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp
index daf397790..32d1393d0 100644
--- a/unsupported/test/FFT.cpp
+++ b/unsupported/test/FFT.cpp
@@ -44,7 +44,7 @@ complex<long double>  promote(long double x) { return complex<long double>( x);
     {
         long double totalpower=0;
         long double difpower=0;
-        cerr <<"idx\ttruth\t\tvalue\n";
+        cerr <<"idx\ttruth\t\tvalue\t|dif|=\n";
         for (size_t k0=0;k0<fftbuf.size();++k0) {
             complex<long double> acc = 0;
             long double phinc = -2.*k0* M_PIl / timebuf.size();
@@ -55,7 +55,7 @@ complex<long double>  promote(long double x) { return complex<long double>( x);
             complex<long double> x = promote(fftbuf[k0]); 
             complex<long double> dif = acc - x;
             difpower += norm(dif);
-            cerr << k0 << "\t" << acc << "\t" <<  x << endl;
+            cerr << k0 << "\t" << acc << "\t" <<  x << "\t" << sqrt(norm(dif)) << endl;
         }
         cerr << "rmse:" << sqrt(difpower/totalpower) << endl;
         return sqrt(difpower/totalpower);
@@ -127,8 +127,9 @@ void test_FFT()
 #endif
 
 #if 1
-  CALL_SUBTEST( test_scalar<float>(45) ); CALL_SUBTEST( test_scalar<double>(45) ); CALL_SUBTEST( test_scalar<long double>(45) );
   CALL_SUBTEST( test_scalar<float>(32) ); CALL_SUBTEST( test_scalar<double>(32) ); CALL_SUBTEST( test_scalar<long double>(32) );
+  CALL_SUBTEST( test_scalar<float>(45) ); CALL_SUBTEST( test_scalar<double>(45) ); CALL_SUBTEST( test_scalar<long double>(45) );
+  CALL_SUBTEST( test_scalar<float>(50) ); CALL_SUBTEST( test_scalar<double>(50) ); CALL_SUBTEST( test_scalar<long double>(50) );
   CALL_SUBTEST( test_scalar<float>(256) ); CALL_SUBTEST( test_scalar<double>(256) ); CALL_SUBTEST( test_scalar<long double>(256) );
   CALL_SUBTEST( test_scalar<float>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<double>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<long double>(2*3*4*5*7) );
 #endif