/Users/mourrain/Devel/mmx/algebramix/include/algebramix/fft_naive.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_naive.hpp
* DESCRIPTION: generic low-level FFT multiplication
* COPYRIGHT  : (C) 2005  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX__FFT_NAIVE__HPP
#define __MMX__FFT_NAIVE__HPP
#include <algebramix/fft_roots.hpp>

namespace mmx {

/******************************************************************************
* Naive root arithmetic
******************************************************************************/

template<typename CC, typename UU, typename SS>
struct roots_helper {
  typedef CC C;
  typedef UU U;
  typedef SS S;

  static U*
  create_roots (nat n) {
    nat k= primitive_root_max_order<C> (2); (void) k;
    VERIFY (k == 0 || n <= k, "maximum order exceeded");
    VERIFY (n >= 2, "size must be at least two");
    U* roots= mmx_new<U> (n);
    for (nat i=0; i<n; i+=2) {
      U temp = primitive_root<C> (n, bit_mirror (i, n));
      roots[i] = temp;
      roots[i+1]= primitive_root<U> (n, i==0? 0: n - bit_mirror (i, n));
    }
    return roots; }

  static void
  destroy_roots (U* u, nat n) {
    mmx_delete<U> (u, n); }

  static inline void
  fft_cross (C* c1, C* c2) {
    C temp= (*c2);
    *c2   = (*c1) - temp;
    *c1   = (*c1) + temp; }

  static inline void
  dfft_cross (C* c1, C* c2, const U* u) {
    C temp= (*u) * (*c2);
    *c2   = (*c1) - temp;
    *c1   = (*c1) + temp; }

  static inline void
  ifft_cross (C* c1, C* c2, const U* u) {
    C temp= *c2;
    *c2   = (*u ) * ((*c1) - temp);
    *c1   = (*c1) + temp; }

  static inline void
  dtft_cross (C* c1, C* c2) {
    static S h= invert (S (2));
    fft_cross (c1, c2);
    *c1 *= h;
    *c2 *= h; }
  
  static inline void
  dtft_cross (C* c1, C* c2, const U* u) {
    static S h= invert (S (2));
    dfft_cross (c1, c2, u);
    *c1 *= h;
    *c2 *= h; }
  
  static inline void
  itft_flip (C* c1, C* c2, const U* u) {
    static S h= invert (S(2));
    C temp= (*u) * (*c2);
    *c1  += (*c1) - temp;
    *c2   = h * ((*c1) - temp); }
  
  static inline void
  itft_flip (C* c1, C* c2) {
    static S h= invert (S(2));
    *c1  += (*c1) - (*c2);
    *c2   = h * ((*c1) - (*c2)); }

  struct fft_mul_sc_op : mul_op {};
};

/******************************************************************************
* The FFT transformer class
******************************************************************************/

template<typename C, typename V= std_roots<C> >
class fft_naive_transformer {
public:
  typedef implementation<vector_linear,vector_naive> NVec;
  typedef typename V::roots_type R;
  typedef typename R::U U;
  typedef typename R::S S;

  nat depth;
  nat len;
  U*  roots;

public:
  inline fft_naive_transformer (nat n):
    depth (log_2 (n)), len (n), roots (R::create_roots (n)) {
    VERIFY (n == ((nat) 1 << depth), "power of two expected"); }
    
  inline ~fft_naive_transformer () {
    R::destroy_roots (roots, len); }

  inline void
  dfft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    // In place direct fft of c[0], c[stride], ..., c[(2^steps-1) stride]
    // Only perform steps from step1 until step2-1
    // If shift != 0, then roots start at roots + (shift<<1)
    for (nat step= step1; step < step2; step++) {
      //mmout << "step " << step << ": " << flush_now;
      if (step == 0 && shift == 0) {
        nat todo= steps - 1;
        C* cc= c;
        for (nat k= 0; k < ((nat) 1<<todo); k++) {
          R::fft_cross (cc, cc + (stride<<todo));
          cc += stride;
        }
      }
      else {
        nat todo= steps - 1 - step;
        C* cc= c;
        U * uu= roots + ((shift >> todo) << 1);
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          for (nat k= 0; k < ((nat) 1<<todo); k++) {
            R::dfft_cross (cc, cc + (stride<<todo), uu);
            cc += stride;
          }
          cc += (stride<<todo);
          uu += 2;
        }
      }
    }
  }

  inline void
  ifft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    // In place inverse fft of c[0], c[stride], ..., c[(2^steps-1) stride]
    // Only perform steps from step2-1 until step1
    // If shift != 0, then roots start at roots + (shift<<1)
    for (int step= step2-1; (int) step >= ((int) step1); step--) {
      //mmout << "step " << step << ": " << flush_now;
      if (step == 0 && shift == 0) {
        nat todo= steps - 1;
        C* cc= c;
        for (nat k= 0; k < ((nat) 1<<todo); k++) {
          R::fft_cross (cc, cc + (stride<<todo));
          cc += stride;
        }
      }
      else {
        nat todo= steps - 1 - step;
        C* cc= c;
        U * uu= roots + 1 + ((shift >> todo) << 1);
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          for (nat k= 0; k < ((nat) 1<<todo); k++) {
            R::ifft_cross (cc, cc + (stride<<todo), uu);
            cc += stride;
          }
          cc += (stride<<todo);
          uu += 2;
        }
      }
    }
  }

  inline void
  dfft (C* c, nat stride, nat shift, nat steps) {
    dfft (c, stride, shift, steps, 0, steps); }

  inline void
  ifft (C* c, nat stride, nat shift, nat steps) {
    ifft (c, stride, shift, steps, 0, steps); }

  inline void
  direct_transform (C* c) {
    dfft (c, 1, 0, depth); }

  inline void
  inverse_transform (C* c, bool divide=true) {
    ifft (c, 1, 0, depth);
    if (divide) {
      S x= binpow (S (2), depth);
      x= invert (x);
      NVec::template vec_unary_scalar<typename R::fft_mul_sc_op> (c, x, len);
    }
  }
};

/******************************************************************************
* The FFT transformer class
******************************************************************************/

template<typename C> inline void
direct_fft (C* dest, nat n) {
  fft_naive_transformer<C> ffter (n);
  ffter.direct_transform (dest);
}

template<typename C> inline void
inverse_fft (C* dest, nat n) {
  fft_naive_transformer<C> ffter (n);
  ffter.inverse_transform (dest);
}

} // namespace mmx
#endif //__MMX__FFT_NAIVE__HPP