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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_schonhage.hpp
* DESCRIPTION: Fast product for any ring with unity
* 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.
******************************************************************************/

// Ref. "Modern Computer Algebra", Chapter 8, Exercise 8.30.

#ifndef __MMX__POLYNOMIAL_SCHONHAGE__HPP
#define __MMX__POLYNOMIAL_SCHONHAGE__HPP
#include <numerix/modulus.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/polynomial_schonhage_strassen.hpp>
#include <algebramix/polynomial_schonhage_triadic.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Variant for Schonhage multiplication
******************************************************************************/

struct schonhage_threshold {};
struct schonhage_internal_threshold {};

template<typename V,
         typename Th1= schonhage_threshold,
         typename Th2= schonhage_internal_threshold>
struct polynomial_schonhage_inc: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef typename V::Positive Positive;
  typedef polynomial_schonhage_inc<typename V::No_simd,Th1,Th2> No_simd;
  typedef polynomial_schonhage_inc<typename V::No_thread,Th1,Th2> No_thread;
  typedef polynomial_schonhage_inc<typename V::No_scaled,Th1,Th2> No_scaled;
};

template<typename F, typename V, typename W, typename Th1, typename Th2>
struct implementation<F,V,polynomial_schonhage_inc<W,Th1,Th2> >:
  public implementation<F,V,W> {};

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_schonhage,
                  polynomial_balanced_tft<
                    polynomial_schonhage_inc<
                      polynomial_karatsuba<V> > >)

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

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

  TMPL static inline void
  mul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
    typedef polynomial_schonhage_strassen_inc<V,Th_rec> V2;
    typedef polynomial_schonhage_triadic_inc<V,Th_rec> V3;
    typedef implementation<polynomial_multiply,V2> Pol2;
    typedef implementation<polynomial_multiply,V3> Pol3;
    typedef modulus_int_naive<8*sizeof(nat)> MV;
    typedef modulus<nat, MV> Modulus;
    if (n1 < Threshold(C,Th) || n2 < Threshold(C,Th)) {
      Fallback::mul (dest, s1, s2, n1, n2);
      return;
    }
    nat n = n1 + n2 - 1;
    nat l = aligned_size<C,V> (n);
    C* temp = mmx_new<C> (l);
    nat k2= Pol2::mul (dest, s1, s2, n1, n2, false);
    if (k2 == 0) return;
    nat k3= Pol3::mul (temp, s1, s2, n1, n2, false);
    if (k3 == 0) {
      Vec::copy (dest, temp, n);
      return;
    }
    nat u2, u3;
    nat c2= ((nat) 1) << k2, c3= binpow ((nat) 3, k3);
    MV::inv_mod (u2, binpow ((nat) 2, k2), Modulus (c3));
    MV::inv_mod (u3, binpow ((nat) 3, k3), Modulus (c2));
    u3= c2 - u3;
    VERIFY (u2*c2 - u3*c3 == 1, "bug");
    Vec::mul (dest, C(u2), n);
    Vec::mul (temp, C(u3), n);
    Vec::sub (dest, temp , n);
    mmx_delete<C> (temp, l); }

  TMPL static inline void
  square (C* dest, const C* s1, nat n1) {
    mul (dest, s1, s1, n1, n1); }

}; // implementation<polynomial_multiply,V,polynomial_schonhage_inc<W,Th,Th_rec> >

/******************************************************************************
* Schonhage Strassen multiplication for generic coefficients
******************************************************************************/

DEFINE_VARIANT (polynomial_generic_schonhage,
                polynomial_ring_dicho<
                  polynomial_schonhage<
                    polynomial_ring_naive<
                      polynomial_naive> > >)

STMPL
struct threshold_helper<generic,
                        polynomial_multiply_threshold<
                          polynomial_generic_schonhage> > {
  typedef fixed_value<nat,16> impl;
};

STMPL
struct threshold_helper<generic, schonhage_threshold> {
  typedef fixed_value<nat,16384> impl;
};

STMPL
struct threshold_helper<generic, schonhage_internal_threshold> {
  typedef fixed_value<nat,128> impl;
};

#define POLYNOMIAL_GENERIC_USES_SCHONHAGE \
STMPL \
struct polynomial_variant_helper<generic> { \
  typedef polynomial_generic_schonhage PV; \
};

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