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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_schonhage_triadic.hpp
* DESCRIPTION: Schonhage's triadic cyclic wrapped 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, Exercise 8.30.

#ifndef __MMX__POLYNOMIAL_SCHONHAGE_TRIADIC__HPP
#define __MMX__POLYNOMIAL_SCHONHAGE_TRIADIC__HPP
#include <algebramix/fft_triadic_naive.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/polynomial_balanced.hpp>

namespace mmx {
#define TMPL template<typename C>

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

struct schonhage_triadic_threshold {};

template<typename V, typename Th= schonhage_triadic_threshold>
struct polynomial_schonhage_triadic_inc: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef typename V::Positive Positive;
  typedef polynomial_schonhage_triadic_inc<typename V::No_simd,Th> No_simd;
  typedef polynomial_schonhage_triadic_inc<typename V::No_thread,Th> No_thread;
  typedef polynomial_schonhage_triadic_inc<typename V::No_scaled,Th> No_scaled;
};

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

DEFINE_VARIANT_1 (typename V, V,
                  polynomial_schonhage_triadic,
                  polynomial_balanced<
                    polynomial_schonhage_triadic_inc<
                      polynomial_karatsuba<V> > >)

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

template<typename V, typename W, typename Th>
struct implementation<polynomial_multiply,V,
                      polynomial_schonhage_triadic_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 m, nat t) {
    nat l = aligned_size<C,V> (m);
    C** dest= mmx_new<C*> (t);
    for (nat j = 0; j < t; j++)
      dest[j]= mmx_new<C> (l);
    return dest; }

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

  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 2n, and dest has size 2t in x^m and 2m in x.
    // We assume that n <= m * t 
    nat i, j, n2= n << 1, m2= m << 1, t2= t << 1;
    for (i = 0, j = 0; i < n2; i += m, j++)
      if (n2-i >= m) {
        Vec::copy  (dest[j], s + i, m);
        Vec::clear (dest[j] + m, m);
      }
      else {
        Vec::copy  (dest[j], s + i, n2 - i);
        Vec::clear (dest[j] + n2 - i, m2 - (n2 - i));
      }
    for (; j < t2; 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) mod (x^(2n) + x^n + 1)
    // dest has size 2n, s has size 2t.
    // We assume that n = m * t
    nat i, n2= n << 1, m2= m << 1, t2= t << 1;
    Vec::clear (dest, n2);
    for (i = 0; i+1 < t2; i++)
      Pol::add (dest + i * m, s[i], m2);
    Pol::add (dest + i * m, s[i]    , m);
    Pol::sub (dest        , s[i] + m, m);
    Pol::sub (dest + n    , s[i] + m, m); }

  TMPL static inline void
  triadic_shift (C* dest, const C* src, nat i, nat m) {
    // multiply by x^i modulo x^(2*m) + x^m + 1
    nat m2= m << 1, m3= m2 + m;
    i= i % m3;
    if (i > m2) {
      Vec::copy (dest         , src + m3 - i, i - m);
      Vec::neg  (dest + i - m , src         , m3 - i);
      Vec::sub  (dest + i - m2, src         , m3 - i);
    }
    else if (i > m) {
      Vec::neg  (dest    , src + m2 - i, m);
      Vec::neg  (dest + m, src + m2 - i, m);
      Vec::add  (dest    , src + m3 - i, i - m);
      Vec::add  (dest + i, src         , m2 - i);
    }
    else {
      Vec::neg  (dest    , src + m2 - i, i);
      Vec::copy (dest + i, src         , m2 - i);
      Vec::sub  (dest + m, src + m2 - i, i); } }
  
  template<typename C>
  struct unptr_helper {};

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

  template<typename Cp>
  struct triadic_roots_helper {
    // modulo x^(2*m) + x^m + 1
    typedef Cp  C;
    typedef nat U;
    typedef typename unptr_helper<Cp>::type CC;
    typedef CC  S;

    static inline nat&
    dyn_m () {
      static nat m= 0;
      return m; }

    static inline C&
    dyn_temp (nat i) {
      static C temp0= NULL;
      static C temp1= NULL;
      static C temp2= NULL;
      static C temp3= NULL;
      switch (i) {
      case 0: return temp0;
      case 1: return temp1;
      case 2: return temp2;
      case 3: return temp3;
      default: ERROR ("index out of range"); return temp0; } }

    static inline nat
    primitive_root (nat t, nat i) {
      if (t == 0) return 1;
      i = i % t;
      nat m3 = 3 * dyn_m ();
      ASSERT (m3 % t == 0, "primitive root out of range");
      return i * (m3 / t); }

    static U*
    create_roots (nat t) {
      nat m2= dyn_m () << 1;
      nat l= aligned_size<CC,V> (m2);
      for (nat i = 0; i <= 3; i++)
        dyn_temp (i)= mmx_new<CC> (l);
      U* roots= mmx_new<U> (t);
      for (nat i = 0; i < t; i++)
        roots[i]  = primitive_root (t, digit_mirror_triadic (i, t));
      return roots; }

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

    static void
    destroy_roots (U* u, nat t) {
      nat m2= dyn_m () << 1;
      nat l= aligned_size<CC,V> (m2);
      for (nat i = 0; i <= 3; i++) {
        C& temp= dyn_temp (i);
        if (temp != NULL) {
          mmx_delete<CC> (temp, l);
          temp= NULL;
        }
      }
      mmx_delete<U> (u, t); }

    static inline void
    dfft_cross (C* c1, C* c2, C* c3, const U* u1, const U* u2, const U* u3) {
      C temp0= dyn_temp (0);
      C temp1= dyn_temp (1);
      C temp2= dyn_temp (2);
      C temp3= dyn_temp (3);
      nat m= dyn_m (), m2= m << 1;

      triadic_shift (temp0, *c3, *u3, m);
      Vec::add (temp0, *c2, m2);
      triadic_shift (temp3, temp0, *u3, m);

      triadic_shift (temp0, *c3, *u2, m);
      Vec::add (temp0, *c2, m2);
      triadic_shift (temp2, temp0, *u2, m);

      triadic_shift (temp0, *c3, *u1, m);
      Vec::add (temp0, *c2, m2);
      triadic_shift (temp1, temp0, *u1, m);

      Vec::add (*c3, *c1, temp3, m2);
      Vec::add (*c2, *c1, temp2, m2);
      Vec::add (*c1, temp1, m2); }

    static inline void
    ifft_cross (C* c1, C* c2, C* c3, const U* u1, const U* u2, const U* u3) {
      C temp0= dyn_temp (0);
      C temp1= dyn_temp (1);
      C temp2= dyn_temp (2);
      C temp3= dyn_temp (3);
      nat m= dyn_m (), m2= m << 1;

      triadic_shift (temp1, *c1, *u1, m);
      Vec::add (*c1, *c2, m2);
      Vec::add (*c1, *c3, m2);

      triadic_shift (temp2, *c2, *u2, m);
      triadic_shift (temp3, *c3, *u3, m);
      triadic_shift (temp0, *c2, *u3, m);
      Vec::add (*c2, temp1, temp2, m2);
      Vec::add (*c2, temp3, m2);

      triadic_shift (temp3, *c3, *u2, m);
      Vec::add (temp0, temp3, m2);
      Vec::add (temp0, temp1, m2);
      triadic_shift (*c3, temp0, *u1, m); }

    static inline void
    fft_shift (C* dest, S v, nat t) {
      nat m2= dyn_m () << 1;
      for (nat i = 0; i < t; i++)
        Vec::mul (dest[i], v, m2); }
  };

  struct triadic_roots {
    template<typename C>
    struct helper {
      typedef triadic_roots_helper<C> roots_type; };
  };

  TMPL static void
  direct_transform (C** b, nat m, nat t) {
    triadic_roots_helper<C*>::dyn_m ()= m;
    fft_triadic_naive_transformer<C*, triadic_roots> ffter (t);
    ffter.direct_transform_triadic (b); } 

  TMPL static void
  inverse_transform (C** b, nat m, nat t, bool shift=true) {
    triadic_roots_helper<C*>::dyn_m ()= m;
    fft_triadic_naive_transformer<C*, triadic_roots> ffter (t);
    ffter.inverse_transform_triadic (b,shift); } 

  TMPL static void
  variable_dilate (C** b, nat eta, nat m, nat t) {
    // substitute eta * y for y in b of size t
    nat m2= m << 1, m3= m2 + m;
    nat l= aligned_size<C,V> (m2);
    C* temp= mmx_new<C> (l);
    nat w= 0;
    for (nat i = 0; i < t; i++) {
      Vec::copy (temp, b[i], m2);
      triadic_shift (b[i], temp, w, m);
      w= (w + eta) % m3;
    }
    mmx_delete<C> (temp, l); }

  TMPL static void
  bivariate_mod (C** dest, const C** h, nat w, nat m, nat t) {
    // dest = h mod (y^t - w)
    nat m2= m << 1, m3= m2 + m;
    w= w % m3;
    for (nat i = 0; i < t; i++) {
      triadic_shift (dest[i], h[i+t], w, m);
      Vec::add (dest[i], h[i], m2); } }

  TMPL static void
  bivariate_crt (C** h, const C** h1, const C** h2, nat w,
                 nat m, nat t, bool shift=true) {
    // h1 = h mod (y^t - w), h2 = h mod (y^t - w^2). 
    nat m2= m << 1, m3= m2 + m, t2= t << 1;
    w= w % m3;
    nat w2= (w << 1) % m3;
    nat l= aligned_size<C,V> (m2);
    C* temp= mmx_new<C> (l);
    for (nat i = 0; i < t; i++) 
      Vec::sub (h[i+t], h2[i], h1[i], m2);
    for (nat i = 0; i < t; i++) {
      triadic_shift (h[i], h1[i], w2, m);
      triadic_shift (temp, h2[i], w , m);
      Vec::sub (h[i], temp, m2);
    }
    for (nat i = 0; i < t2; i++) {
      triadic_shift (temp, h[i], w, m);
      Vec::add (h[i], temp, m2);
      Vec::add (h[i], temp, m2);
      if (shift)
        Vec::mul (h[i], invert (C(3)), m2);
    }
    mmx_delete<C> (temp, l); }
  
public:
  TMPL static inline nat
  mul_triadic (C* dest, const C* src, nat n, bool shift=true) {
    // dest *= src mod (x^(2*n) + x^n + 1)
    nat n2= n << 1;
    nat u = next_power_of_three (n);
    ASSERT (u != 0, "maximum size exceeded");
    ASSERT (n == u, "power of three expected");
    if (n <= max ((nat) 9, (nat) Threshold(C,Th))) {
      nat l= aligned_size<C,V> (2*n2-1);
      C* temp= mmx_new<C> (l);
      Inner::mul (temp, dest, src, n2, n2);
      Vec::copy (dest    , temp         , n2);
      Vec::add  (dest    , temp + n2 + n, n - 1);
      Vec::sub  (dest    , temp + n2    , n);
      Vec::sub  (dest + n, temp + n2    , n);
      mmx_delete<C> (temp, l);
      return 0;
    }
    else {
      nat r = 0;
      nat k = log_3 (n);
      nat m = binpow (3, (k+1) >> 1), m2= m << 1, m3= m2 + m;
      nat t = n / m, t2= t << 1;
      nat eta= (t == m) ? 1 : 3;
      // mmout << "n = " << n << ", m = " << m << ", t = " << t << "\n";
      C** hh= bivariate_new<C> (m2, t2);
      C** gg= bivariate_new<C> (m2, t2);
      C** h1= bivariate_new<C> (m2, t);
      C** h2= bivariate_new<C> (m2, t);
      C** g = bivariate_new<C> (m2, t);
      bivariate_encode (hh, src , n, m, t);
      bivariate_encode (gg, dest, n, m, t);
      for (nat j = 1; j <= 2; j++) {
        C** h= j == 1 ? h1 : h2;
        bivariate_mod (h, (const C**) hh, (j*eta*t) % m3, m, t);
        bivariate_mod (g, (const C**) gg, (j*eta*t) % m3, m, t);
        variable_dilate (h, (j*eta) % m3, m, t);
        variable_dilate (g, (j*eta) % m3, m, t);
        direct_transform (h, m, t);
        direct_transform (g, m, t);
        for (nat i = 0; i < t; i++)
          r = mul_triadic (h[i], g[i], m, shift);
        inverse_transform (h, m, t, shift);
        variable_dilate (h, m3 - ((j*eta) % m3), m, t);
      }
      bivariate_delete<C> (g , m2 , t);
      bivariate_delete<C> (gg, m2, t2);
      bivariate_crt (hh, (const C**) h1, (const C**) h2, (eta*t) % m3, m, t, shift);
      bivariate_delete<C> (h1, m2, t);
      bivariate_delete<C> (h2, m2, t); 
      bivariate_decode (dest, (const C**) hh, n, m, t);
      bivariate_delete<C> (hh, m2, t2);
      return (shift) ? 0 : r + 1 + (k >> 1); }}

  TMPL static inline nat
  mul (C* dest, const C* s1, const C* s2, nat n1, nat n2, bool shift=true) {
    nat len = n1 + n2 - 1;
    nat n = next_power_of_three ((len + 1) >> 1);
    nat nn= n << 1;
    nat l = aligned_size<C,V> (nn);
    C* t1 = mmx_new<C> (l);
    C* t2 = mmx_new<C> (l);
    Vec::copy (t1, s1, n1); Vec::clear (t1 + n1, nn - n1);
    Vec::copy (t2, s2, n2); Vec::clear (t2 + n2, nn - n2);
    nat k = mul_triadic (t1, t2, n, shift);
    Vec::copy (dest, t1, len);
    mmx_delete<C> (t1, l);
    mmx_delete<C> (t2, l);
    return k; }

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

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