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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_tft.hpp
* DESCRIPTION: FFT based algorithms
* 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__POLYNOMIAL_TFT__HPP
#define __MMX__POLYNOMIAL_TFT__HPP
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/polynomial_ring_dicho.hpp>
#include <algebramix/fft_blocks.hpp>
#include <algebramix/fft_simd.hpp>
#include <algebramix/fft_truncated.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* TFT based variant polynomials
******************************************************************************/

struct tft_threshold {};
template<typename V> struct polynomial_balanced_tft;

template<typename C>
struct threshold_helper<C,tft_threshold> {
  typedef fixed_value<nat,128> impl;
};

template<typename C, typename V>
struct threshold_helper<polynomial<C,V>,tft_threshold> {
  typedef fixed_value<nat,32> impl;
};

template<typename C, typename U, typename V, typename W>
struct threshold_helper<modular<modulus<polynomial<C,U>,V>,W>,
                        tft_threshold> {
typedef fixed_value<nat,32> impl;
};

template<typename V, typename T= tft_threshold>
struct polynomial_tft_inc: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef typename V::Positive Positive;
  typedef polynomial_tft_inc<typename V::No_simd,T> No_simd;
  typedef polynomial_tft_inc<typename V::No_thread,T> No_thread;
  typedef polynomial_tft_inc<typename V::No_scaled,T> No_scaled;
};

template<typename F, typename V, typename W, typename Th>
struct implementation<F,V,polynomial_tft_inc<W,Th> >:
  public implementation<F,V,W> {};

DEFINE_VARIANT_1 (typename V, V, polynomial_tft,
                  polynomial_dicho<
                   polynomial_balanced_tft<
                     polynomial_tft_inc<
                       polynomial_karatsuba<V> > > >)

STMPL
struct polynomial_variant_helper<modular<modulus<nat>,
                                         modular_fixed<nat,0x0c0000001> > > {
  typedef polynomial_tft<polynomial_naive> PV;
};

template<typename C, typename V>
struct polynomial_variant_helper<polynomial<C,polynomial_tft<V> > > {
  typedef polynomial_gcd_ring_dicho<polynomial_tft<V> > PV;
};

/******************************************************************************
* Multiplication
******************************************************************************/

template<typename V, typename W, typename Th>
struct implementation<polynomial_multiply,V,polynomial_tft_inc<W,Th> >:
  public implementation<polynomial_linear,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,V> Pol;
  typedef implementation<polynomial_multiply,W> Fallback;

TMPL static inline void
mul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  typedef fft_blocks_transformer<C, fft_simd_transformer<C>,
    8, 5, 10, 16> FFTer;
  typedef fft_truncated_transformer<C,FFTer> Tfter;
  //typedef fft_truncated_transformer<C> Tfter;
  
  if (n1 < Threshold(C,Th) || n2 < Threshold(C,Th))
    Fallback::mul (dest, s1, s2, n1, n2);
  else {
    nat n= n1 + n2 - 1;
    nat m= next_power_of_two (n);
    nat l= aligned_size<C,V> (m);
    C* temp0= mmx_new<C> (l);
    C* temp1= mmx_new<C> (l);
    C* tempd= mmx_new<C> (l);
    Tfter tfter (n);
    Vec::copy (temp0, s1, n1);
    Vec::clear (temp0+n1, n-n1);
    Vec::copy (temp1, s2, n2);
    Vec::clear (temp1+n2, n-n2);
    tfter.direct_transform (temp0);
    tfter.direct_transform (temp1);
    Vec::mul (tempd, temp0, temp1, n);
    tfter.inverse_transform (tempd);
    Vec::copy (dest, tempd, n);
    mmx_delete<C> (temp0, l);
    mmx_delete<C> (temp1, l);
    mmx_delete<C> (tempd, l);
  }
}

TMPL static inline void
square (C* dest, const C* s1, nat n1) {
  typedef fft_truncated_transformer<C> Tfter;
  if (n1 < Threshold(C,Th))
    Fallback::square (dest, s1, n1);
  else {
    nat n= 2 * n1 - 1;
    nat m= next_power_of_two (n);
    nat l= aligned_size<C,V> (m);
    C* temp0= mmx_new<C> (l);
    C* tempd= mmx_new<C> (l);
    Tfter tfter (n);
    Vec::copy (temp0, s1, n1);
    Vec::clear (temp0+n1, n-n1);
    tfter.direct_transform (temp0);
    Vec::mul (tempd, temp0, temp0, n);
    tfter.inverse_transform (tempd);
    Vec::copy (dest, tempd, n);
    mmx_delete<C> (temp0, l);
    mmx_delete<C> (tempd, l);
  }
}

}; // implementation<polynomial_multiply,V,polynomial_tft_inc<W,Th> >

/******************************************************************************
* Degree balancing
******************************************************************************/

template<typename V>
struct polynomial_balanced_tft: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_balanced_tft<typename V::Positive> Positive;
  typedef polynomial_balanced_tft<typename V::No_simd> No_simd;
  typedef polynomial_balanced_tft<typename V::No_thread> No_thread;
  typedef polynomial_balanced_tft<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,polynomial_balanced_tft<W> >:
  public implementation<F,V,W> {};

template<typename V, typename W>
struct implementation<polynomial_multiply,V,polynomial_balanced_tft<W> >:
  public implementation<polynomial_linear,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,W> Pol;
  typedef implementation<polynomial_multiply,W> Fallback;

TMPL static inline void
mul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  if (n1 == 0 && n2 == 0)
    return;
  if (n1 == 0 || n2 == 0) {
    Pol::clear (dest, n1+n2-1);
    return;
  }
  if (n2 < n1) {
    mul (dest, s2, s1, n2, n1);
    return;
  }
  nat n= n1+n2-1;
  nat m= next_power_of_two (2*n1 - 1);
  if (n1 <= n2 && n >= m) {
    Pol::clear (dest, n);
    nat l= aligned_size<C,V> (m);
    C* temp= mmx_new<C> (l);
    while (n1+n2-1 >= m) {
      Fallback::mul (temp, s1, s2, n1, m-n1+1);
      Vec::add (dest, temp, m);
      n2 -= m-n1+1; s2 += m-n1+1; dest += m-n1+1;
    }
    if (n2 > 0) {
      mul (temp, s1, s2, n1, n2);
      Vec::add (dest, temp, n1+n2-1);
    }
    mmx_delete<C> (temp, l);
    return;
  }
  Fallback::mul (dest, s1, s2, n1, n2);
}
  
TMPL static inline void
square (C* dest, const C* s1, nat n1) {
  Fallback::square (dest, s1, n1);
}

}; // implementation<polynomial_multiply,V,polynomial_balanced_tft<W> >

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