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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_schonhage_strassen.hpp
* DESCRIPTION: Schonhage-Strassen fast product
* 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, Algorithm 8.20.

#ifndef __MMX__POLYNOMIAL_SCHONHAGE_STRASSEN__HPP
#define __MMX__POLYNOMIAL_SCHONHAGE_STRASSEN__HPP
#include <algebramix/polynomial.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/polynomial_tft.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Variant for Schonhage-Strassen multiplication
******************************************************************************/

struct schonhage_strassen_threshold {};

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

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_schonhage_strassen,
                  polynomial_balanced_tft<
                    polynomial_schonhage_strassen_inc<
                      polynomial_karatsuba<V> > >)

/******************************************************************************
* Schonhage-Strassen multiplication
******************************************************************************/

template<typename V, typename W, typename Th>
struct implementation<polynomial_multiply,V,
                      polynomial_schonhage_strassen_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> Inner;

private:  
  TMPL static inline C**
  bivariate_new (nat m2, nat t) {
    nat l = aligned_size<C ,V> (m2);
    nat ly= aligned_size<C*,V> (t);
    C** dest= mmx_new<C*> (ly);
    for (nat j = 0; j < t; j++)
      dest[j]= mmx_new<C> (l);
    return dest; }

  TMPL static inline void
  bivariate_delete (C** dest, nat m2, nat t) {
    nat l = aligned_size<C ,V> (m2);
    nat ly= aligned_size<C*,V> (t);
    for (nat j = 0; j < t; j++)
      mmx_delete<C> (dest[j], l);
    mmx_delete<C*> (dest, ly); }

  TMPL static inline void
  bivariate_encode (C** dest, const C*s, nat n, nat m, nat t) {
    // dest (x, x^m) = s (x)
    // s has size n, and dest has size t in x^m and 2m in x.
    // We assume that n <= m * t 
    nat i, j, m2= m << 1;
    for (i = 0, j = 0; i < n; i += m, j++)
      if (n-i >= m) {
        Vec::copy  (dest[j], s + i, m);
        Vec::clear (dest[j] + m, m);
      }
      else {
        Vec::copy  (dest[j], s + i, n - i);
        Vec::clear (dest[j] + n - i, m2 - (n - i));
      }
    for (; j < t; j++)
      Vec::clear (dest[j], m2); }
  
  TMPL static inline void
  bivariate_decode (C* dest, const C** s, nat n, nat m, nat t) {
    // dest (x) = s (x, x^m), dest has size n
    // We assume that n = m * t
    nat i, m2= m << 1;
    Vec::clear (dest, n);
    for (i = 0; i+1 < t; i++)
      Pol::add (dest + i * m, s[i], m2);
    Pol::add (dest + i * m, s[i]    , m);
    Pol::sub (dest        , s[i] + m, m); }

  TMPL static inline void
  negative_cyclic_shift (C* dest, const C* src, nat m2, nat i) {
    // multiply by x^i modulo x^m2 + 1
    bool negate= ((i / m2) & 1);
    i= i % m2;
    if (negate) {
      Vec::neg   (dest + i, src         , m2 - i);
      Vec::copy  (dest    , src + m2 - i, i);
    }
    else {
      Vec::copy  (dest + i, src         , m2 - i);
      Vec::neg   (dest    , src + m2 - i, i); } }
  
  TMPL static inline void
  negative_cyclic_shift (C* dest, nat m2, nat i) {
    nat l= aligned_size<C,V> (m2);
    C* temp= mmx_new<C> (l);
    Vec::copy (temp, dest, m2);
    negative_cyclic_shift (dest, temp, m2, i);
    mmx_delete<C> (temp, l); }

  template<typename C>
  struct unptr_helper {};

  template<typename C>
  struct unptr_helper<C*> {
    typedef C type; };

  template<typename Cp>
  struct negative_cyclic_roots_helper {
    // modulo x^m2 + 1
    typedef Cp  C;
    typedef nat U;
    typedef typename unptr_helper<Cp>::type CC;
    typedef CC  S;

    static inline nat&
    dyn_modulus () {
      static nat m2= 0;
      return m2; }

    static inline C&
    get_temp () {
      static C temp= NULL;
      return temp; }

    static inline nat
    primitive_root (nat t, nat i) {
      ASSERT (t != 0, "unexpected zero root order");
      i = i % t;
      nat m2 = dyn_modulus ();
      return i * (m2 << 1) / t; }

    static U*
    create_roots (nat t) {
      nat m2= dyn_modulus ();
      nat l= aligned_size<CC,V> (m2);
      get_temp ()= mmx_new<CC> (l);
      U* roots= mmx_new<U> (t);
      for (nat i = 0; i < t; i += 2) {
        roots[i]  = primitive_root (t, bit_mirror (i, t));
        roots[i+1]= primitive_root (t, i == 0 ? 0 : t - bit_mirror (i, t));
      }
      return roots; }

    static void
    destroy_roots (U* u, nat t) {
      C& temp= get_temp ();
      if (temp != NULL) {
        nat m2= dyn_modulus ();
        nat l= aligned_size<CC,V> (m2);
        mmx_delete<CC> (temp, l);
        temp= NULL;
      }
      mmx_delete<U> (u, t); }

    static inline void
    fft_cross (C* c1, C* c2) {
      C temp= get_temp ();
      nat m2= dyn_modulus ();
      Vec::copy (temp, *c2, m2);
      Vec::sub  (*c2, *c1, temp, m2);
      Vec::add  (*c1, temp, m2); }
    
    static inline void
    dfft_cross (C* c1, C* c2, const U* u) {
      C temp= get_temp ();
      nat m2= dyn_modulus ();
      negative_cyclic_shift (temp, *c2, m2, *u);
      Vec::sub  (*c2, *c1, temp, m2);
      Vec::add  (*c1, temp, m2); }

    static inline void
    ifft_cross (C* c1, C* c2, const U* u) {
      C temp= get_temp ();
      nat m2= dyn_modulus ();
      Vec::copy (temp, *c1, m2);
      Vec::add  (*c1, *c2, m2);
      Vec::sub  (*c2, temp, m2);
      negative_cyclic_shift (*c2, m2, (*u) + m2); }

    static inline void
    dtft_cross (C* c1, C* c2) {
      static CC h= invert (CC(2));
      nat m2= dyn_modulus ();
      fft_cross (c1, c2);
      Vec::mul (*c1, h, m2);
      Vec::mul (*c2, h, m2); }
    
    static inline void
    dtft_cross (C* c1, C* c2, const U* u) {
      static CC h= invert (CC(2));
      nat m2= dyn_modulus ();
      dfft_cross (c1, c2, u);
      Vec::mul (*c1, h, m2);
      Vec::mul (*c2, h, m2); }
    
    static inline void
    itft_flip (C* c1, C* c2, const U* u) {
      static CC h= invert (CC(2));
      C temp= get_temp ();
      nat m2= dyn_modulus ();
      negative_cyclic_shift (temp, *c2, m2, *u);
      Vec::add (*c1, *c1, m2);
      Vec::sub (*c1, temp, m2);
      Vec::sub (*c2, *c1, temp, m2);
      Vec::mul (*c2, h, m2); }
    
    static inline void
    itft_flip (C* c1, C* c2) {
      static CC h= invert (CC(2));
      nat m2= dyn_modulus ();
      Vec::add (*c1, *c1, m2);
      Vec::sub (*c1, *c2, m2);
      Vec::sub (*c2, *c1, *c2, m2);
      Vec::mul (*c2, h, m2); }
    
  struct fft_mul_sc_op : mul_op {
    static inline void
    set_op (C& x, const S& y) {
      Vec::mul (x, y, dyn_modulus ()); }
  };
  };

  TMPL
  struct negative_cyclic_roots {
    typedef negative_cyclic_roots_helper<C> roots_type;
  };

  TMPL static void
  direct_transform (C** b, nat m2, nat t) {
    negative_cyclic_roots_helper<C*>::dyn_modulus ()= m2;
    fft_naive_transformer<C*, negative_cyclic_roots<C*> > ffter (t);
    ffter.direct_transform (b); } 

  TMPL static void
  inverse_transform (C** b, nat m2, nat t, bool divide=true) {
    negative_cyclic_roots_helper<C*>::dyn_modulus ()= m2;
    fft_naive_transformer<C*, negative_cyclic_roots<C*> > ffter (t);
    ffter.inverse_transform (b, divide); } 

  TMPL static void
  direct_transform_truncated (C** b, nat m2, nat len) {
    negative_cyclic_roots_helper<C*>::dyn_modulus ()= m2;
    fft_truncated_transformer<C*,
      fft_naive_transformer<C*, negative_cyclic_roots<C*> > > ffter (len);
    ffter.dtft (b, 1, len, 0); }

  TMPL static void
  inverse_transform_truncated (C** b, nat m2, nat len) {
    nat i;
    for (i = len; i < next_power_of_two (len); i++)
      Vec::clear (b[i], m2);
    negative_cyclic_roots_helper<C*>::dyn_modulus ()= m2;
    fft_truncated_transformer<C*,
      fft_naive_transformer<C*, negative_cyclic_roots<C*> > > ffter (len);
    ffter.itft (b, 1, len, 0);
    C h= invert (C(next_power_of_two (len)));
    for (i = 0; i < len; i++)
      Vec::mul (b[i], h, m2);
    for (; i < next_power_of_two (len); i++)
      Vec::clear (b[i], m2); }

  TMPL static void
  variable_dilate (C** b, nat eta, nat m2, nat t) {
    // substitute eta * y for y in b of size t
    nat m4= m2 << 1;
    nat w= 0;
    for (nat i = 0; i < t; i++) {
      negative_cyclic_shift (b[i], m2, w);
      w= (w + eta) % m4; } }

public:
  TMPL static inline nat
  mul_negative_cyclic (C* dest, const C* src, nat n, bool shift=true) {
    // dest *= src mod (x^n + 1)
    nat u= next_power_of_two (n);
    ASSERT (u != 0, "maximum size exceeded");
    ASSERT (n == u, "power of two expected");
    if (n <= max ((nat) 4, (nat) Threshold(C,Th))) {
      nat l= aligned_size<C,V> (2*n-1);
      C* temp= mmx_new<C> (l);
      Inner::mul (temp, dest, src, n, n);
      Vec::copy (dest, temp, n);
      Vec::sub  (dest, temp + n, n - 1);      
      mmx_delete<C> (temp, l);
      return 0;
    }
    else {
      nat r = 0;
      nat k = log_2 (n);
      nat m = (nat) 1 << (k >> 1);
      nat t = u / m;
      nat m2= m << 1;
      C** b1= bivariate_new<C> (m2, t);
      C** b2= bivariate_new<C> (m2, t);
      bivariate_encode (b1, src , n, m, t);
      bivariate_encode (b2, dest, n, m, t);
      nat eta= (t == m2) ? 1 : 2;
      variable_dilate (b1, eta, m2, t);
      variable_dilate (b2, eta, m2, t);
      direct_transform (b1, m2, t);
      direct_transform (b2, m2, t);
      for (nat i = 0; i < t; i++)
        r= mul_negative_cyclic (b1[i], b2[i], m2, shift);
      bivariate_delete<C> (b2, m2, t);
      inverse_transform (b1, m2, t, shift);
      variable_dilate (b1, (m2 << 1) - eta, m2, t);
      bivariate_decode (dest, (const C**) b1, n, m, t);
      bivariate_delete<C> (b1, m2, t);
      return r + ((k+1) >> 1); } }

  TMPL static inline void
  mul_negative_cyclic_truncated (C* dest, const C* src, nat len) {
    // dest *= src mod (x^n + 1)
    // TFT is used in first round
    nat n= next_power_of_two (len);
    ASSERT (n != 0, "maximum size exceeded");
    if (n <= max ((nat) 4, (nat) Threshold(C,Th))) {
      nat l= aligned_size<C,V> (2*n-1);
      C* temp= mmx_new<C> (l);
      Inner::mul (temp, dest, src, n, n);
      Vec::copy (dest, temp, n);
      Vec::sub  (dest, temp + n, n - 1);      
      mmx_delete<C> (temp, l);
      return;
    }
    else {
      nat k = log_2 (n);
      nat m = (nat) 1 << (k >> 1);
      nat t = n / m;
      nat r = (len + m - 1) / m;
      nat m2= m << 1;
      C** b1= bivariate_new<C> (m2, t);
      C** b2= bivariate_new<C> (m2, t);
      bivariate_encode (b1, src , n, m, t);
      bivariate_encode (b2, dest, n, m, t);
      nat eta= (t == m2) ? 1 : 2;
      variable_dilate (b1, eta, m2, r);
      variable_dilate (b2, eta, m2, r);
      direct_transform_truncated (b1, m2, r);
      direct_transform_truncated (b2, m2, r);
      for (nat i = 0; i < r; i++)
        mul_negative_cyclic (b1[i], b2[i], m2);
      bivariate_delete<C> (b2, m2, t);
      inverse_transform_truncated (b1, m2, r);
      variable_dilate (b1, (m2 << 1) - eta, m2, r);
      bivariate_decode (dest, (const C**) b1, n, m, t);
      bivariate_delete<C> (b1, m2, t); } }

  TMPL static inline nat
  mul (C* dest, const C* s1, const C* s2, nat n1, nat n2, bool shift=true) {
    nat ret = 0;
    nat len = n1 + n2 - 1;
    nat n = next_power_of_two (len);
    nat l = aligned_size<C,V> (n);
    C* t1 = mmx_new<C> (l);
    C* t2 = mmx_new<C> (l);
    Vec::copy (t1, s1, n1); Vec::clear (t1 + n1, n - n1);
    Vec::copy (t2, s2, n2); Vec::clear (t2 + n2, n - n2);
    if (shift)
      mul_negative_cyclic_truncated (t1, t2, len);
    else
      ret= mul_negative_cyclic (t1, t2, n, shift);
    Vec::copy (dest, t1, len);
    mmx_delete<C> (t1, l);
    mmx_delete<C> (t2, l);
    return ret; }

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

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