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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_truncated.hpp
* DESCRIPTION: generic low-level univariate TFT multiplication
* COPYRIGHT  : (C) 2005  Joris van der Hoeven and 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_TRUNCATED__HPP
#define __MMX__FFT_TRUNCATED__HPP
#include <basix/vector_naive.hpp>
#include <algebramix/fft_naive.hpp>

namespace mmx {

#define TMPL template<class C>

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

template<typename C, typename Ffter= fft_naive_transformer<C> >
class fft_truncated_transformer {
public:
  typedef implementation<vector_linear,vector_naive> NVec;
  typedef typename Ffter::R R;
  typedef typename Ffter::U U;
  typedef typename Ffter::S S;
  nat len;
  Ffter* ffter;

public:
  inline fft_truncated_transformer (nat s): len (s) {
    nat n= next_power_of_two (s);
    VERIFY(s <= n, "maximum size exceeded");
    ffter= new Ffter (n); }

  inline ~fft_truncated_transformer () { delete ffter; }

  // Direct transformation
  inline void
  fft_cross_range (C* lp, C* rp, nat stride, nat shift, nat nr) {
    U* a= ffter->roots + (shift << 1);
    if (shift == 0)
      for (nat i=nr; i!=0; i--) {
        R::fft_cross (lp, rp);
        lp += stride; rp +=stride;
      }
    else
      for (nat i=nr; i!=0; i--) {
        R::dfft_cross (lp, rp, a);
        lp += stride; rp += stride;
      }
  }

  inline void
  dtft (C* c, nat stride, nat s, nat shift, nat steps) {
    // In place direct tft of c[0], c[stride], ..., c[(s-1) stride]
    // Roots start at roots + (shift<<1)
    if (s == 0 || steps == 0) return;
    nat todo= steps - 1;
    nat w   = (nat) 1 << todo;
    fft_cross_range (c, c + (stride<<todo), stride, shift>>todo, w);
    steps--;
    if (s >= w) {
      ffter->dfft (c, stride, shift, steps);
      s -= w;  shift += w >> 1;  c += stride<<todo;
    }
    dtft (c, stride, s, shift, steps);
  }

  // Inverse transformation
  inline void
  itft_flip_range (C* lp, C* rp, nat stride, nat shift, nat nr) {
    U* a= ffter->roots + (shift << 1);
    if (shift == 0)
      for (nat i=nr; i!=0; i--) {
        R::itft_flip (lp, rp);
        lp += stride; rp += stride;
      }
    else
      for (nat i=nr; i!=0; i--) {
        R::itft_flip (lp, rp, a);
        lp += stride; rp += stride;
      }
  }
  
  inline void
  tft_cross_range (C* lp, C* rp, nat stride, nat shift, nat nr) {
    U* a= ffter->roots + (shift << 1);
    if (shift == 0)
      for (nat i=nr; i!=0; i--) {
        R::dtft_cross (lp, rp);
        lp += stride; rp += stride;
      }
    else
      for (nat i=nr; i!=0; i--) {
        R::dtft_cross (lp, rp, a);
        lp += stride; rp += stride;
      }
  }
  
  inline void
  ifft_cross_range (C* lp, C* rp, nat stride, nat shift, nat nr) {
    U* a= ffter->roots + (shift << 1);
    if (shift == 0)
      for (nat i=nr; i!=0; i--) {
        R::fft_cross (lp, rp);
        lp += stride; rp += stride;
      }
    else
      for (nat i=nr; i!=0; i--) {
        R::ifft_cross (lp, rp, a+1);
        lp += stride; rp += stride;
      }
  }
  
  inline void
  itft (C* c, nat stride, nat s, nat shift, nat steps) {
    // In place inverse tft of c[0], c[stride], ..., c[(s-1) stride]
    // Roots start at roots + (shift<<1)
    if (s == 0 || steps == 0) return;
    nat todo= steps - 1;
    nat l   = (nat) 1 << steps;
    nat w   = (nat) 1 << todo;
    if (s < w) {
      tft_cross_range (c + stride*s, c + stride*(s+w),
                       stride, shift>>todo, w-s);
      itft (c, stride, s, shift, steps-1);
      itft_flip_range (c, c + (stride<<todo), stride, shift>>todo, s);
    }
    else {
      ffter->ifft (c, stride, shift, steps-1);
      itft_flip_range (c + stride*(s-w), c + stride*s,
                       stride, shift>>todo, l-s);
      if (s == w) return;
      itft (c + (stride<<todo), stride, s-w, shift + (w>>1), steps-1);
      ifft_cross_range (c, c + (stride<<todo), stride, shift>>todo, s-w);
    }
  }

  inline void
  dtft (C* c, nat stride, nat s, nat shift) {
    dtft (c, stride, s, shift, ffter->depth); }

  inline void
  itft (C* c, nat stride, nat s, nat shift) {
    itft (c, stride, s, shift, ffter->depth); }

  inline void
  direct_transform (C* c) {
    nat n= ((nat) 1) << ffter->depth;
    if (len == n)
      ffter->direct_transform (c);
    else {
      NVec::clear (c + len, n - len);
      dtft (c, 1, len, 0);
      NVec::clear (c + len, n - len);
    }
  }

  inline void
  inverse_transform (C* c) {
    nat n= ((nat) 1) << ffter->depth;
    if (len == n)
      ffter->inverse_transform (c);
    else {
      NVec::clear (c + len, n - len);
      itft (c, 1, len, 0);
      S x= invert (binpow (S (2), ffter->depth));
      NVec::template vec_unary_scalar<typename R::fft_mul_sc_op> (c, x, len);
      NVec::clear (c + len, n - len);
    }
  }
};

#undef TMPL
} // namespace mmx
#endif //__MMX__FFT_TRUNCATED__HPP