rfft.cpp

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

      t5.re=x5[0];
      s0.re=t0.re+t3.re;
      s1.re=t0.re-t3.re;
      s2.re=t1.re+t2.re;
      s3.re=t1.re-t2.re;
      s4.re=t4.re+t5.re;
      s5.re=t4.re-t5.re;
      y0[ 0]=s0.re+s2.re+s4.re;     
      y1[-1]=s1.re+A*(s3.re-s5.re); 
      y1[ 0]=      B*(s2.re-s4.re); 
      y2[-1]=s0.re-A*(s2.re+s4.re); 
      y2[ 0]=      B*(s3.re+s5.re); 
      y3[-1]=s1.re-s3.re+s5.re;     
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  t0.re=x0[2*k-1];
	  t0.im=x0[2*k  ];
	  t1.re=x1[2*k-1];
	  t1.im=x1[2*k  ];
	  t2.re=x2[2*k-1];
	  t2.im=x2[2*k  ];
	  t3.re=x3[2*k-1];
	  t3.im=x3[2*k  ];
	  t4.re=x4[2*k-1];
	  t4.im=x4[2*k  ];
	  t5.re=x5[2*k-1];
	  t5.im=x5[2*k  ];
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w2.re*w.re-w2.im*w.im;
	  w3.im=w2.re*w.im+w2.im*w.re; 
	  w4.re=w3.re*w.re-w3.im*w.im;
	  w4.im=w3.re*w.im+w3.im*w.re; 
	  w5.re=w4.re*w.re-w4.im*w.im;
	  w5.im=w4.re*w.im+w4.im*w.re; 
	  s0.re=t0.re;
	  s0.im=t0.im;
	  s1.re=t1.re*w.re-t1.im*w.im;
	  s1.im=t1.re*w.im+t1.im*w.re;
	  s2.re=t2.re*w2.re-t2.im*w2.im;
	  s2.im=t2.re*w2.im+t2.im*w2.re;
	  s3.re=t3.re*w3.re-t3.im*w3.im;
	  s3.im=t3.re*w3.im+t3.im*w3.re;
	  s4.re=t4.re*w4.re-t4.im*w4.im;
	  s4.im=t4.re*w4.im+t4.im*w4.re;
	  s5.re=t5.re*w5.re-t5.im*w5.im;
	  s5.im=t5.re*w5.im+t5.im*w5.re;
	  t0.re=s0.re+s3.re;
	  t0.im=s0.im+s3.im;
	  t1.re=s1.re+s2.re+s4.re+s5.re;
	  t1.im=s1.im+s2.im+s4.im+s5.im;
	  t2.re=s1.re-s2.re-s4.re+s5.re;
	  t2.im=s1.im-s2.im-s4.im+s5.im;
	  t3.re=s1.re+s2.re-s4.re-s5.re;
	  t3.im=s1.im+s2.im-s4.im-s5.im;
	  t4.re=s1.re-s2.re+s4.re-s5.re;
	  t4.im=s1.im-s2.im+s4.im-s5.im;
	  t5.re=s0.re-s3.re;
	  t5.im=s0.im-s3.im;
	  s0.re=t0.re+t1.re;
	  s0.im=t0.im+t1.im;
	  s5.re=t5.re-t2.re;
	  s5.im=t5.im-t2.im;
	  s1.re=t5.re+A*t2.re;
	  s1.im=t5.im+A*t2.im;
	  s3.re=B*t3.re;
	  s3.im=B*t3.im;
	  s2.re=t0.re-A*t1.re;
	  s2.im=t0.im-A*t1.im;
	  s4.re=B*t4.re;
	  s4.im=B*t4.im;
	  y0[ 2*k-1] = s0.re;
	  y0[ 2*k  ] = s0.im;
	  y1[ 2*k-1] = s1.re-s3.im;
	  y1[ 2*k  ] = s1.im+s3.re;
	  y1[-2*k-1] = s1.re+s3.im;
	  y1[-2*k  ] =-s1.im+s3.re;
	  y2[ 2*k-1] = s2.re-s4.im;
	  y2[ 2*k  ] = s2.im+s4.re;
	  y2[-2*k-1] = s2.re+s4.im;
	  y2[-2*k  ] =-s2.im+s4.re;
	  y3[-2*k-1] = s5.re;
	  y3[-2*k  ] =-s5.im;
	  z.re=omega.re*w.re-omega.im*w.im+w.re;
	  z.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=z.re;
	  w.im=z.im;
        }
      if(L%2==0)
        {
	  t0.re=x0[L-1];
	  t1.re=x1[L-1];
	  t2.re=x2[L-1];
	  t3.re=x3[L-1];
	  t4.re=x4[L-1];
	  t5.re=x5[L-1];
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w2.re*w.re-w2.im*w.im;
	  w3.im=w2.re*w.im+w2.im*w.re; 
	  w4.re=w3.re*w.re-w3.im*w.im;
	  w4.im=w3.re*w.im+w3.im*w.re; 
	  w5.re=w4.re*w.re-w4.im*w.im;
	  w5.im=w4.re*w.im+w4.im*w.re; 
	  s0.re=t0.re;
	  s0.im=0.0;
	  s1.re=t1.re*w.re;
	  s1.im=t1.re*w.im;
	  s2.re=t2.re*w2.re;
	  s2.im=t2.re*w2.im;
	  s3.re=t3.re*w3.re;
	  s3.im=t3.re*w3.im;
	  s4.re=t4.re*w4.re;
	  s4.im=t4.re*w4.im;
	  s5.re=t5.re*w5.re;
	  s5.im=t5.re*w5.im;
	  t0.re=s0.re+s3.re;
	  t0.im=s0.im+s3.im;
	  t1.re=s0.re-s3.re;
	  t1.im=s0.im-s3.im;
	  t2.re=s1.re+s2.re;
	  t2.im=s1.im+s2.im;
	  t3.re=s1.re-s2.re;
	  t3.im=s1.im-s2.im;
	  t4.re=s4.re+s5.re;
	  t4.im=s4.im+s5.im;
	  t5.re=s4.re-s5.re;
	  t5.im=s4.im-s5.im;
	  s0.re=t0.re+t2.re+t4.re;
	  s0.im=t0.im+t2.im+t4.im;
	  s1.re=t1.re+A*(t3.re-t5.re)-B*(t2.im-t4.im);
	  s1.im=t1.im+A*(t3.im-t5.im)+B*(t2.re-t4.re);
	  s2.re=t0.re-A*(t2.re+t4.re)-B*(t3.im+t5.im);
	  s2.im=t0.im-A*(t2.im+t4.im)+B*(t3.re+t5.re);
	  y0[L-1] = s0.re;
	  y0[L  ] = s0.im;
	  y1[L-1] = s1.re;
	  y1[L  ] = s1.im;
	  y2[L-1] = s2.re;
	  y2[L  ] = s2.im;
        }
    }
}

/// Internal radix-6 butterfly backward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis_radix6(real *y,real *x,int L,int N)
{
  complex z,w1,w2,w3,w4,w5;
  complex t0,t1,t2,t3,t4,t5;
  complex s0,s1,s2,s3,s4,s5;
  const real A= 0.5L;
  const real B=-0.5L*root3;
  real *x0,*x1,*x2,*x3,*x4,*x5;
  real *y0,*y1,*y2,*y3;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      x4=x+j*L+4*N*L;
      x5=x+j*L+5*N*L;
      y0=y+j*6*L;
      y1=y+j*6*L+2*L;
      y2=y+j*6*L+4*L;
      y3=y+j*6*L+6*L;
      t0.re=y0[ 0];
      t1.re=y1[-1];
      t1.im=y1[ 0];
      t2.re=y2[-1];
      t2.im=y2[ 0];
      t3.re=y3[-1];
      s0.re=t0.re+t3.re;
      s1.re=2.0*(t1.re+t2.re);
      s1.im=2.0*(t1.im+t2.im);
      s2.re=2.0*(t1.re-t2.re);
      s2.im=2.0*(t1.im-t2.im);
      s3.re=t0.re-t3.re;
      s4.re=s3.re+A*s2.re;
      s4.im=B*s1.im;
      s5.re=s0.re-A*s1.re;
      s5.im=B*s2.im;
      x0[0]=s0.re+s1.re;
      x1[0]=s4.re+s4.im;
      x2[0]=s5.re+s5.im;
      x3[0]=s3.re-s2.re;
      x4[0]=s5.re-s5.im;
      x5[0]=s4.re-s4.im;
      w1.re=omega.re+1.0;                   
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  t0.re= y0[ 2*k-1];
	  t0.im= y0[ 2*k  ];
	  t1.re= y1[ 2*k-1];
	  t1.im= y1[ 2*k  ];
	  t3.re= y1[-2*k-1];
	  t3.im=-y1[-2*k  ];
	  t2.re= y2[ 2*k-1];
	  t2.im= y2[ 2*k  ];
	  t4.re= y2[-2*k-1];
	  t4.im=-y2[-2*k  ];
	  t5.re= y3[-2*k-1];
	  t5.im=-y3[-2*k  ];
	  s0.re=t0.re+t5.re;
	  s0.im=t0.im+t5.im;
	  s1.re=t0.re-t5.re;
	  s1.im=t0.im-t5.im;
	  s2.re=t1.re-t2.re+t3.re-t4.re;
	  s2.im=t1.im-t2.im+t3.im-t4.im;
	  s3.re=t1.re+t2.re+t3.re+t4.re;
	  s3.im=t1.im+t2.im+t3.im+t4.im;
	  s4.re=t1.re+t2.re-t3.re-t4.re;
	  s4.im=t1.im+t2.im-t3.im-t4.im;
	  s5.re=t1.re-t2.re-t3.re+t4.re;
	  s5.im=t1.im-t2.im-t3.im+t4.im;
	  t0.re=s0.re+s3.re;
	  t0.im=s0.im+s3.im;
	  t1.re=s1.re+A*s2.re;
	  t1.im=s1.im+A*s2.im;
	  t2.re=s0.re-A*s3.re;
	  t2.im=s0.im-A*s3.im;
	  t3.re=s1.re-s2.re;
	  t3.im=s1.im-s2.im;
	  t4.re=B*s5.re;
	  t4.im=B*s5.im;
	  t5.re=B*s4.re;
	  t5.im=B*s4.im;
	  s0.re=t0.re;
	  s0.im=t0.im;
	  s3.re=t3.re;
	  s3.im=t3.im;
	  s1.re=t1.re+t5.im;
	  s1.im=t1.im-t5.re;
	  s2.re=t2.re+t4.im;
	  s2.im=t2.im-t4.re;
	  s5.re=t1.re-t5.im;
	  s5.im=t1.im+t5.re;
	  s4.re=t2.re-t4.im;
	  s4.im=t2.im+t4.re;
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  w3.re=w2.re*w1.re-w2.im*w1.im;
	  w3.im=w2.re*w1.im+w2.im*w1.re; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re;
	  x0[2*k-1]=s0.re;
	  x0[2*k  ]=s0.im;
	  x1[2*k-1]=s1.re*w1.re+s1.im*w1.im;
	  x1[2*k  ]=s1.im*w1.re-s1.re*w1.im;
	  x2[2*k-1]=s2.re*w2.re+s2.im*w2.im;
	  x2[2*k  ]=s2.im*w2.re-s2.re*w2.im;
	  x3[2*k-1]=s3.re*w3.re+s3.im*w3.im;
	  x3[2*k  ]=s3.im*w3.re-s3.re*w3.im;
	  x4[2*k-1]=s4.re*w4.re+s4.im*w4.im;
	  x4[2*k  ]=s4.im*w4.re-s4.re*w4.im;
	  x5[2*k-1]=s5.re*w5.re+s5.im*w5.im;
	  x5[2*k  ]=s5.im*w5.re-s5.re*w5.im;
	  z.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  z.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=z.re;
	  w1.im=z.im;
        }
      if(L%2==0)
        {
	  t0.re=y0[L-1];
	  t0.im=y0[L  ];
	  t1.re=y1[L-1];
	  t1.im=y1[L  ];
	  t2.re=y2[L-1];
	  t2.im=y2[L  ];
	  t5.re= t0.re;
	  t5.im=-t0.im;
	  t4.re= t1.re;
	  t4.im=-t1.im;
	  t3.re= t2.re;
	  t3.im=-t2.im;
	  s0.re=t0.re+t3.re;
	  s0.im=t0.im+t3.im;
	  s1.re=t1.re+t2.re+t4.re+t5.re;
	  s1.im=t1.im+t2.im+t4.im+t5.im;
	  s2.re=t0.re-t3.re;
	  s2.im=t0.im-t3.im;
	  s3.re=t1.re-t2.re-t4.re+t5.re;
	  s3.im=t1.im-t2.im-t4.im+t5.im;
	  s4.re=t1.re+t2.re-t4.re-t5.re;
	  s4.im=t1.im+t2.im-t4.im-t5.im;
	  s5.re=t1.re-t2.re+t4.re-t5.re;
	  s5.im=t1.im-t2.im+t4.im-t5.im;
	  t0.re=s0.re+s1.re;
	  t0.im=s0.im+s1.im;
	  t1.re=s2.re+A*s3.re+B*s4.im;
	  t1.im=s2.im+A*s3.im-B*s4.re;
	  t2.re=s0.re-A*s1.re+B*s5.im;
	  t2.im=s0.im-A*s1.im-B*s5.re;
	  t3.re=s2.re-s3.re;
	  t3.im=s2.im-s3.im;
	  t4.re=s0.re-A*s1.re-B*s5.im;
	  t4.im=s0.im-A*s1.im+B*s5.re;
	  t5.re=s2.re+A*s3.re-B*s4.im;
	  t5.im=s2.im+A*s3.im+B*s4.re;
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  w3.re=w2.re*w1.re-w2.im*w1.im;
	  w3.im=w2.re*w1.im+w2.im*w1.re; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re;
	  x0[L-1]=t0.re;
	  x1[L-1]=t1.re*w1.re+t1.im*w1.im;
	  x2[L-1]=t2.re*w2.re+t2.im*w2.im;
	  x3[L-1]=t3.re*w3.re+t3.im*w3.im;
	  x4[L-1]=t4.re*w4.re+t4.im*w4.im;
	  x5[L-1]=t5.re*w5.re+t5.im*w5.im;
            
        }
    }
}

/// Internal radix-8 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix8(real *x,real *y,int L,int N)
{
  complex w1,w2,w3,w4,w5,w6,w7;
  complex s0,s1,s2,s3,s4,s5,s6,s7;
  complex t0,t1,t2,t3,t4,t5,t6,t7;
  real *x0,*x1,*x2,*x3,*x4,*x5,*x6,*x7;
  real *y0,*y1,*y2,*y3,*y4;
  const real c=0.5L*root2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x0+N*L;
      x2=x1+N*L;
      x3=x2+N*L;
      x4=x3+N*L;  
      x5=x4+N*L;
      x6=x5+N*L;
      x7=x6+N*L;
      y0=y    +j*8*L;
      y1=y+2*L+j*8*L;
      y2=y+4*L+j*8*L;
      y3=y+6*L+j*8*L;
      y4=y+8*L+j*8*L;
      t0.re=x0[0];
      t1.re=x1[0];
      t2.re=x2[0];
      t3.re=x3[0];
      t4.re=x4[0];
      t5.re=x5[0];
      t6.re=x6[0];
      t7.re=x7[0];
      s0.re=t0.re+t4.re;
      s1.re=t2.re+t6.re;
      s2.re=t1.re+t5.re;
      s3.re=t3.re+t7.re;
      s4.re=t0.re-t4.re;
      s5.re=t2.re-t6.re;
      s6.re=t1.re-t5.re;
      s7.re=t3.re-t7.re;    
      t0.re=  (s0.re+s1.re);
      t1.re=  (s2.re+s3.re);
      t2.re=  s4.re;
      t2.im= -s5.re;
      t3.re=c*(s6.re-s7.re);
      t3.im=c*(-s6.re-s7.re);
      t4.re=  (s0.re-s1.re); 
      t5.im= -(s2.re-s3.re);
      t6.re=  s4.re;
      t6.im=  s5.re;
      t7.re=-t3.re;
      t7.im=t3.im;
      y0[ 0]=t0.re+t1.re;
      y1[-1]=t2.re+t3.re;
      y1[ 0]=t2.im+t3.im;
      y2[-1]=t4.re;
      y2[ 0]=t5.im;
      y3[-1]=t6.re+t7.re;
      y3[ 0]=t6.im+t7.im;
      y4[-1]=t0.re-t1.re;
      w1.re=omega.re+1.0;                   
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  w3.re=w2.re*w1.re-w2.im*w1.im;
	  w3.im=w2.re*w1.im+w2.im*w1.re; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re; 
	  w6.re=w5.re*w1.re-w5.im*w1.im;
	  w6.im=w5.re*w1.im+w5.im*w1.re; 
	  w7.re=w6.re*w1.re-w6.im*w1.im;
	  w7.im=w6.re*w1.im+w6.im*w1.re; 
	  t0.re=x0[2*k-1];
	  t0.im=x0[2*k  ];
	  s1.re=x1[2*k-1];
	  s1.im=x1[2*k  ];
	  s2.re=x2[2*k-1];
	  s2.im=x2[2*k  ];
	  s3.re=x3[2*k-1];
	  s3.im=x3[2*k  ];
	  s4.re=x4[2*k-1];
	  s4.im=x4[2*k  ];
	  s5.re=x5[2*k-1];
	  s5.im=x5[2*k  ];
	  s6.re=x6[2*k-1];
	  s6.im=x6[2*k  ];
	  s7.re=x7[2*k-1];
	  s7.im=x7[2*k  ];
	  t1.re=s1.re*w1.re-s1.im*w1.im;
	  t1.im=s1.im*w1.re+s1.re*w1.im;
	  t2.re=s2.re*w2.re-s2.im*w2.im;
	  t2.im=s2.im*w2.re+s2.re*w2.im;
	  t3.re=s3.re*w3.re-s3.im*w3.im;
	  t3.im=s3.im*w3.re+s3.re*w3.im;
	  t4.re=s4.re*w4.re-s4.im*w4.im;
	  t4.im=s4.im*w4.re+s4.re*w4.im;
	  t5.re=s5.re*w5.re-s5.im*w5.im;
	  t5.im=s5.im*w5.re+s5.re*w5.im;      
	  t6.re=s6.re*w6.re-s6.im*w6.im;
	  t6.im=s6.im*w6.re+s6.re*w6.im;
	  t7.re=s7.re*w7.re-s7.im*w7.im;
	  t7.im=s7.im*w7.re+s7.re*w7.im;
	  s0.re=t0.re+t4.re;
	  s0.im=t0.im+t4.im;
	  s1.re=t2.re+t6.re;
	  s1.im=t2.im+t6.im;
	  s2.re=t1.re+t5.re;
	  s2.im=t1.im+t5.im;
	  s3.re=t3.re+t7.re;
	  s3.im=t3.im+t7.im;
	  s4.re=t0.re-t4.re;
	  s4.im=t0.im-t4.im;
	  s5.re=t2.re-t6.re;
	  s5.im=t2.im-t6.im;
	  s6.re=t1.re-t5.re;
	  s6.im=t1.im-t5.im;
	  s7.re=t3.re-t7.re;
	  s7.im=t3.im-t7.im;
	  t0.re=  (s0.re+s1.re);
	  t0.im=  (s0.im+s1.im);
	  t1.re=  (s2.re+s3.re);
	  t1.im=  (s2.im+s3.im);
	  t2.re=  (s4.re+s5.im);
	  t2.im=  (s4.im-s5.re);
	  t3.re=c*(s6.re+s6.im-s7.re+s7.im);
	  t3.im=c*(s6.im-s6.re-s7.im-s7.re);
	  t4.re=  (s0.re-s1.re); 
	  t4.im=  (s0.im-s1.im); 
	  t5.re=  (s2.im-s3.im); 
	  t5.im= -(s2.re-s3.re); 
	  t6.re=  (s4.re-s5.im); 
	  t6.im=  (s4.im+s5.re); 
	  t7.re=c*(-s6.re+s6.im+s7.re+s7.im);
	  t7.im=c*(-s6.im-s6.re+s7.im-s7.re);
	  y0[ 2*k-1]= (t0.re+t1.re);
	  y0[ 2*k  ]= (t0.im+t1.im);
	  y1[ 2*k-1]= (t2.re+t3.re);
	  y1[ 2*k  ]= (t2.im+t3.im);
	  y2[ 2*k-1]= (t4.re+t5.re);
	  y2[ 2*k  ]= (t4.im+t5.im);
	  y3[ 2*k-1]= (t6.re+t7.re);
	  y3[ 2*k  ]= (t6.im+t7.im);
	  y1[-2*k-1]= (t6.re-t7.re);
	  y1[-2*k  ]=-(t6.im-t7.im);
	  y2[-2*k-1]= (t4.re-t5.re);
	  y2[-2*k  ]=-(t4.im-t5.im);
	  y3[-2*k-1]= (t2.re-t3.re);
	  y3[-2*k  ]=-(t2.im-t3.im);
	  y4[-2*k-1]= (t0.re-t1.re);
	  y4[-2*k  ]=-(t0.im-t1.im);
	  t0.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  t0.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=t0.re;
	  w1.im=t0.im;
        }
      if(L%2==0)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  w3.re=w2.re*w1.re-w2.im*w1.im;
	  w3.im=w2.re*w1.im+w2.im*w1.re; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re; 
	  w6.re=w5.re*w1.re-w5.im*w1.im;
	  w6.im=w5.re*w1.im+w5.im*w1.re; 
	  w7.re=w6.re*w1.re-w6.im*w1.im;
	  w7.im=w6.re*w1.im+w6.im*w1.re; 
	  t0.re=x0[L-1];
	  t0.im=0.0;
	  t1.re=x1[L-1]*w1.re;
	  t1.im=x1[L-1]*w1.im;
	  t2.re=x2[L-1]*w2.re;
	  t2.im=x2[L-1]*w2.im;
	  t3.re=x3[L-1]*w3.re;
	  t3.im=x3[L-1]*w3.im;
	  t4.re=x4[L-1]*w4.re;
	  t4.im=x4[L-1]*w4.im;
	  t5.re=x5[L-1]*w5.re;
	  t5.im=x5[L-1]*w5.im;
	  t6.re=x6[L-1]*w6.re;
	  t6.im=x6[L-1]*w6.im;
	  t7.re=x7[L-1]*w7.re;
	  t7.im=x7[L-1]*w7.im;
	  s0.re=t0.re+t4.re;
	  s0.im=t0.im+t4.im;
	  s1.re=t2.re+t6.re;
	  s1.im=t2.im+t6.im;
	  s2.re=t1.re+t5.re;
	  s2.im=t1.im+t5.im;
	  s3.re=t3.re+t7.re;
	  s3.im=t3.im+t7.im;
	  s4.re=t0.re-t4.re;
	  s4.im=t0.im-t4.im;
	  s5.re=t2.re-t6.re;
	  s5.im=t2.im-t6.im;
	  s6.re=t1.re-t5.re;
	  s6.im=t1.im-t5.im;
	  s7.re=t3.re-t7.re;
	  s7.im=t3.im-t7.im;
	  t0.re=  ( s0.re+s1.re);
	  t0.im=  ( s0.im+s1.im);
	  t1.re=  ( s2.re+s3.re);
	  t1.im=  ( s2.im+s3.im);
	  t2.re=  ( s4.re+s5.im);
	  t2.im=  ( s4.im-s5.re);
	  t3.re=c*( s6.re+s6.im-s7.re+s7.im);
	  t3.im=c*( s6.im-s6.re-s7.im-s7.re);
	  t4.re=  ( s0.re-s1.re);
	  t4.im=  ( s0.im-s1.im);
	  t5.re=  ( s2.im-s3.im);
	  t5.im=  (-s2.re+s3.re);
	  t6.re=  ( s4.re-s5.im);
	  t6.im=  ( s4.im+s5.re);
	  t7.re=c*(-s6.re+s6.im+s7.re+s7.im);
	  t7.im=c*(-s6.im-s6.re+s7.im-s7.re);
	  y0[L-1] = (t0.re+t1.re);
	  y0[L  ] = (t0.im+t1.im);
	  y1[L-1] = (t2.re+t3.re);
	  y1[L  ] = (t2.im+t3.im);
	  y2[L-1] = (t4.re+t5.re);
	  y2[L  ] = (t4.im+t5.im);