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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_triadic_naive.hpp
* DESCRIPTION: generic low-level triadic FFT
* COPYRIGHT  : (C) 2008  Gregoire Lecerf
*******************************************************************************
* 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_TRIADIC_NAIVE__HPP
#define __MMX__FFT_TRIADIC_NAIVE__HPP
#include <basix/vector_naive.hpp>
#include <algebramix/fft_roots.hpp>

namespace mmx {

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

template<typename CC, typename UU=CC, typename SS=CC>
struct roots_triadic_helper {
  typedef implementation<vector_linear,vector_naive> NVec;
  typedef CC C; // type of entries
  typedef UU U; // type of roots
  typedef SS S; // type of scalars

  static U*
  create_roots (nat n) {
    U* roots= mmx_new<U> (n);
    for (nat i=0; i<n; i++)
      roots[i] = primitive_root<C> (n, digit_mirror_triadic (i, n));
    return roots; }

  static U*
  create_stoor (nat n) {
    U* stoor= mmx_new<U> (n);
    for (nat i=0; i<n; i++)
      stoor[i]= primitive_root<C> (n, i==0 ? 0 :
                                   n - digit_mirror_triadic (i, n));
    return stoor; }

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

  static inline void
  dfft_cross (C* c1, C* c2, C* c3, const U* u1, const U* u2, const U* u3) {
    C temp1= *c1, temp2= *c2;
    *c1 = temp1 + (*u1) * (temp2 + (*u1) * (*c3));
    *c2 = temp1 + (*u2) * (temp2 + (*u2) * (*c3));
    *c3 = temp1 + (*u3) * (temp2 + (*u3) * (*c3));
  }

  static inline void
  ifft_cross (C* c1, C* c2, C* c3, const U* u1, const U* u2, const U* u3) {
    C temp1= *c1, temp2= *c2;
    *c1 = temp1 + temp2 + (*c3);
    *c2 = (*u1) * temp1 + (*u2) * temp2 + (*u3) * (*c3);
    *c3 = (*u1) * ((*u1) * temp1 + (*u3) * temp2 + (*u2) * (*c3));
  }

  static inline void
  fft_shift (C* dest, S v, nat t) {
    NVec::mul (dest, C(v), t); }
};

struct std_roots_triadic {
  template<typename C>
  struct helper {
    typedef cached_roots_helper<roots_triadic_helper<C> > roots_type;
  };
};

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

template<typename C, typename VV= std_roots_triadic>
class fft_triadic_naive_transformer {
public:
  typedef VV V;
  typedef typename V::template helper<C>::roots_type R;
  typedef typename R::U U;
  typedef typename R::S S;

  nat depth;
  nat len;
  U*  roots; // roots in triadic digit reverse order
  U*  stoor; // inverses of the roots

public:
  inline fft_triadic_naive_transformer (nat n):
    depth (log_3 (n)), len (n), roots (R::create_roots (n)),
                                stoor (R::create_stoor (n)) {
    VERIFY (n == binpow ((nat) 3, depth), "power of three expected"); }
  inline ~fft_triadic_naive_transformer () {
    R::destroy_roots (roots, len);
    R::destroy_roots (stoor, len); }

  template<typename CC> inline void
  dfft_triadic (CC* c, nat stride, nat shift,
                nat steps, nat step1, nat step2) {
    // In place direct triadic fft of
    // c[0], c[stride], ..., c[(3^steps-1) stride]
    // Only performs steps from step1 until step2-1
    // Roots start at roots + shift
    typedef typename V::template helper<CC>::roots_type RR;
    for (nat step= step1; step < step2; step++) {
      //mmout << "step " << step << ": " << flush_now;
      nat todo= binpow (3, steps - 1 - step);
      CC* cc  = c;
      U * uu  = roots + (shift / todo);
      for (nat j= 0; j < (nat) binpow (3, step); j++) {
        for (nat k= 0; k < todo; k++) {
          RR::dfft_cross (cc, cc + stride * todo,
                          cc + ((stride * todo) << 1), uu, uu+1, uu+2);
          cc += stride;
        }
        cc += (stride * todo) << 1;
        uu += 3;
      }
    }
  }

  template<typename CC> inline void
  ifft_triadic (CC* c, nat stride, nat shift,
                nat steps, nat step1, nat step2) {
    // In place inverse triadic fft of
    // c[0], c[stride], ..., c[(3^steps-1) stride]
    // Only performs steps from step2-1 until step1
    // Roots start at roots + shift
    typedef typename V::template helper<CC>::roots_type RR;
    for (int step= step2-1; step >= ((int) step1); step--) {
      //mmout << "step " << step << ": " << flush_now;
      nat todo= binpow (3, steps - 1 - step);
      CC* cc  = c;
      U * uu  = stoor + (shift / todo);
      for (nat j= 0; j < (nat) binpow (3, step); j++) {
        for (nat k= 0; k < todo; k++) {
          RR::ifft_cross (cc, cc + stride * todo,
                          cc + ((stride * todo) << 1), uu, uu+1, uu+2);
          cc += stride;
        }
        cc += (stride * todo) << 1;
        uu += 3;
      }
    }
  }

  template<typename CC> inline void
  dfft_triadic (CC* c, nat stride, nat shift, nat steps) {
    dfft_triadic (c, stride, shift, steps, 0, steps);
  }

  template<typename CC> inline void
  ifft_triadic (CC* c, nat stride, nat shift, nat steps) {
    ifft_triadic (c, stride, shift, steps, 0, steps);
  }

  inline void
  direct_transform_triadic (C* c) {
    dfft_triadic (c, 1, 0, depth);
  }

  inline void
  inverse_transform_triadic (C* c, bool shift=true) {
    ifft_triadic (c, 1, 0, depth);
    if (shift)
      R::fft_shift (c, invert (S (len)), len);
  }
};

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

template<typename C> inline void
direct_fft_triadic (C* dest, nat n) {
  fft_triadic_naive_transformer<C> ffter (n);
  ffter.direct_transform_triadic (dest);
}

template<typename C> inline void
inverse_fft_triadic (C* dest, nat n) {
  fft_triadic_naive_transformer<C> ffter (n);
  ffter.inverse_transform_triadic (dest);
}
} // namespace mmx
#endif //__MMX__FFT_TRIADIC_NAIVE__HPP