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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_fast.hpp
* DESCRIPTION: Fast univariate power series
* COPYRIGHT  : (C) 2006  Joris van der Hoeven
*******************************************************************************
* 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__SERIES_FAST__HPP
#define __MMX__SERIES_FAST__HPP
#include <algebramix/series.hpp>
#include <algebramix/series_vector.hpp>
#include <algebramix/fkt_transform.hpp>
#include <algebramix/fft_naive.hpp>

namespace mmx {
#define TMPL template<typename C,typename V>
#define Series series<C,V>
#define Series_rep series_rep<C,V>
#define Vector vector<C>
#define Series_vector series<Vector,V>
#define Series_vector_rep series_rep<Vector,V>

/******************************************************************************
* Variant
******************************************************************************/

struct series_fast {};

template<typename F, typename V>
struct implementation<F,V,series_fast>:
    public implementation<F,V,series_naive> {};

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

template<typename U>
struct implementation<series_multiply,U,series_fast>:
  public implementation<series_abstractions,U>
{

// Setting up the different levels for fast multiplication
#define TRANSFORM_NAIVE      0
#define TRANSFORM_KARATSUBA  1
#define TRANSFORM_FFT        2

template <typename C>
struct level_info {
MMX_ALLOCATORS
  typedef C* C_ptr;

  nat n;          // size of blocks
  nat k;          // multiplier
  nat tsz;        // size of transformed blocks
  nat type;       // type of the transformation
  C_ptr head[2];  // buffer for transformed heads
  C_ptr tail[2];  // buffer for transformed tails

  level_info () {
    head[0]= head[1]= NULL;
    tail[0]= tail[1]= NULL;
  }
  ~level_info () {
    if (head[0] != NULL) mmx_classical_delete<C> (head[0]);
    if (head[1] != NULL) mmx_classical_delete<C> (head[1]);
    if (tail[0] != NULL) mmx_classical_delete<C> (tail[0]);
    if (tail[1] != NULL) mmx_classical_delete<C> (tail[1]);
  }
};

static nat
get_inside_multiplier (nat n) {
  if (n <= (1 << 7)) return 16;
  if (n <= (1 << 9)) return 4;
  if (n <= (1 << 10)) return 8;
  if (n <= (1 << 15)) return 16;
  if (n <= (1 << 23)) return 32;
  return 64;
}

static nat
get_border_multiplier (nat n) {
  if (n <= (1 << 4)) return 2;
  if (n <= (1 << 8)) return 4;
  if (n <= (1 << 13)) return 8;
  if (n <= (1 << 21)) return 16;
  return 32;
}

  /* Better for semi-relaxed multiplication
static nat
get_border_multiplier (nat n) {
  if (n <= (1 << 4)) return 2;
  if (n <= (1 << 6)) return 4;
  if (n <= (1 << 9)) return 8;
  if (n <= (1 << 17)) return 16;
  return 32;
}
  */

  /* Better for (semi-)relaxed mult. over F_p and k-bit compl. floats
static nat
get_border_multiplier (nat n) {
  if (n <= (1 << 3)) return 2;
  if (n <= (1 << 4)) return 4;
  if (n <= (1 << 8)) return 8;
  if (n <= (1 << 16)) return 16;
  return 32;
}
  */

static vector<nat>
determine_sizes (nat n) {
  if (n <= 16) return vector<nat> (1, 1);
  vector<nat> v;
  nat k= get_inside_multiplier (n);
  // nat k= get_inside_multiplier (n*2); for F_p and k-bit compl. floats
  n /= k;
  while (n != 1) {
    k= get_border_multiplier (n);
    v << k;
    n /= k;
  }
  v << 1;
  for (nat i=0; i<N(v)/2; i++)
    swap (v[i], v[N(v) - 1 - i]);
  return v;
}


// Fast multiplication
//static nat mul_count= 0;

TMPL
class nrelax_mul_series_rep: public Series_rep {
public:
  typedef typename series_polynomial_helper<C,V >::PV PV;
  typedef implementation<polynomial_linear,PV> Pol;
  typedef fkt_package<polynomial_naive> Fkt;
  
protected:
  Series f[2];           // the series being multiplied
  nat sh[2];             // sh[i]=1 if f[1-i] is fast and 0 otherwise
  nat capacity[2];       // capacity for inner transforms
  nat nr_levels;         // number of levels
  level_info<C>* info;   // information about each level
  nat xnr_levels;        // length of xinfo
  level_info<C>* xinfo;  // storage of info

public:
  nrelax_mul_series_rep (const Series& f2, const Series& g2, nat n):
    Series_rep (CF(f2))
  {
    f[0]= f2; f[1]= g2;
    sh[0]= 1; sh[1]= 1;
    //sh[0]= 1; sh[1]= 0;
    const vector<nat> v= determine_sizes (n);
    nr_levels= N(v);
    info= mmx_new<level_info<C> > (nr_levels);
    nat sz= 1;
    for (nat level=0; level< nr_levels; level++) {
      bool last= (level == nr_levels - 1);
      nat n= sz * v[level];
      nat k= (last? 1: v[level+1]);
      info[level].n = n;
      info[level].k = k;
      nat tsz, type;
      if (v[level] == 1) {
        tsz = 1;
        type= TRANSFORM_NAIVE;
      }
      else if (v[level] == 2) {
        tsz = 3 * info[level-1].tsz;
        type= TRANSFORM_KARATSUBA;
      }
      else {
        tsz= 2 * n;
        type= TRANSFORM_FFT;
      }
      info[level].tsz = tsz;
      info[level].type= type;
      //mmout << "Level " << level << "\n";
      //mmout << "  " << "n   = " << n << "\n";
      //mmout << "  " << "k   = " << k << "\n";
      //mmout << "  " << "tsz = " << tsz << "\n";
      //mmout << "  " << "type= " << type << "\n";
      if (last) {
        info[level].tail[0]= mmx_classical_new<C> (2 * tsz);
        info[level].tail[1]= mmx_classical_new<C> (2 * tsz);
        capacity[0]= 2;
        capacity[1]= 2;
      }
      else {
        if (sh[0] != 0) {
          info[level].head[0]= mmx_classical_new<C> (k * tsz);
          info[level].tail[1]= mmx_classical_new<C> (k * tsz);
        }
        if (sh[1] != 0) {
          info[level].head[1]= mmx_classical_new<C> (k * tsz);
          info[level].tail[0]= mmx_classical_new<C> (k * tsz);
        }
      }
      sz *= v[level];
      //mmout << "  " << "head= " << info[level].head[0] << ", " << info[level].head[1] << "\n";
      //mmout << "  " << "tail= " << info[level].tail[0] << ", " << info[level].tail[1] << "\n";
    }
    xinfo= info;
    xnr_levels= nr_levels;
    if (sh[0] == 0 && sh[1] == 0) {
      info += nr_levels-1;
      nr_levels= 1;
    }
  }

  ~nrelax_mul_series_rep () {
    this->l= allocated (this->l); // not very nice, but necessary correction
    mmx_delete<level_info<C> > (xinfo, xnr_levels); }

  syntactic expression (const syntactic& z) const {
    return flatten (f[0], z) * flatten (f[1], z); }

  nat allocated (nat l) {
    if (l == 0) return 0;
    else return l + 2 * info[nr_levels-1].n; }
  void Set_order (nat l2) { // not very nice, but OK...
    if (l2 <= this->l) return;
    nat old_allocated= allocated (this->l);
    nat new_allocated= allocated (l2);
    C* b= mmx_new<C> (new_allocated);
    Pol::copy (b, this->a, old_allocated);
    Pol::clear (b + old_allocated, new_allocated - old_allocated);
    mmx_delete<C> (this->a, old_allocated);
    this->a= b;
    this->l= l2;
  }

  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f[0], l);
    increase_order (f[1], l); }

  void direct_transform (C* dest, nat ld, const C* src, nat ls, nat type) {
    //mmout << "    ld= " << ld << ", ls= " << ls << ", type= " << type << "\n";
    if (type == TRANSFORM_NAIVE) dest[0]= src[0];
    else if (type == TRANSFORM_KARATSUBA) {
      Pol::copy (dest, src, ls);
      Fkt::direct_fkt (dest, ls, ld);
    }
    else {
      Pol::copy (dest, src, ls);
      Pol::clear (dest + ls, ls);
      fft_naive_transformer<C> ffter (ld);
      ffter.direct_transform (dest);
    }
  }

  void inverse_transform (C* dest, nat n, nat tsz, nat type) {
    if (type == TRANSFORM_NAIVE);
    else if (type == TRANSFORM_KARATSUBA) Fkt::inverse_fkt (dest, tsz, n);
    else {
      fft_naive_transformer<C> ffter (tsz);
      ffter.inverse_transform (dest);
    }
  }

  void direct_transform (nat which) {
    for (nat level= 0; level < nr_levels; level++) {
      typedef C* C_ptr;
      nat    n   = info[level].n;
      nat    tsh = sh[0] + sh[1];
      nat    cur = this->n + tsh - n * sh[1-which];
      if (cur % n != 0) break;
      if (sh[1-which] != 0 && this->n + 1 < n) break;
      //mmout << "  Transform " << cur/n << " (" << which << ") at level " << level << "\n";
      bool   last= (level == nr_levels - 1);
      nat    k   = info[level].k;
      nat    tsz = info[level].tsz;
      nat    type= info[level].type;
      C_ptr& head= info[level].head[which];
      C_ptr& tail= info[level].tail[which];
      nat Cur= cur / n;
      nat Mod= Cur % k;
      //mmout << "    Mod= " << Mod << "\n";
      if (f[which]->n < this->n + 1 && which == 1) {
        mmerr << "\n>>> n= " << this->n
             << " and only " << f[which]->n << " terms\n";
        ASSERT (f[which]->n >= this->n + 1, "insufficient number of terms");
      }
      const C* seg= (sh[1-which] == 0?
                     f[which] (this->n, this->n + n):
                     f[which] (this->n + 1 - n, this->n + 1));
      //mmout << "    seg= "; Pol::print (mmout, seg, n); mmout << "\n";
      if (last) {
        nat Cap= capacity[which];
        //mmout << "    Cur= " << cur << ", Cap= " << Cap << "\n";
        if (Cur >= Cap) {
          nat old_cap= Cap * tsz;
          nat new_cap= old_cap << 1;
          C* a= mmx_classical_new<C> (new_cap);
          Pol::copy (a, tail, old_cap);
          Pol::clear (a + old_cap, new_cap - old_cap);
          mmx_classical_delete<C> (tail);
          tail= a;
          capacity[which]= Cap << 1;
          //mmout << "    New capacity= " << capacity[which] << "\n";
        }
        //mmout << "    -> tail " << Cur << "\n";
        direct_transform (tail + Cur * tsz, tsz, seg, n, type);
        //mmout << "    tr= "; Pol::print (mmout, tail + Cur * tsz, tsz); mmout << "\n";
      }
      else if (Cur < k * sh[which]) {
        //mmout << "    -> head " << Mod << "\n";
        direct_transform (head + Mod * tsz, tsz, seg, n, type);
        //mmout << "    tr= "; Pol::print (mmout, head + Mod * tsz, tsz); mmout << "\n";
      }
      else if (tail != NULL) {
        //mmout << "    -> tail " << Mod << "\n";
        direct_transform (tail + Mod * tsz, tsz, seg, n, type);
        //mmout << "    tr= "; Pol::print (mmout, tail + Mod * tsz, tsz); mmout << "\n";
      }
    }
  }

  inline void accumulate (C* dest, const C* s1, const C* s2, nat len) {
    //mul_count += len;
    //mmout << "      s1= "; Pol::print (mmout, s1, len); mmout << "\n";
    //mmout << "      s2= "; Pol::print (mmout, s2, len); mmout << "\n";
    Pol::mul_add (dest, s1, s2, len);
    //mmout << "      d = "; Pol::print (mmout, dest, len); mmout << "\n";
  }

  C next () {
    //mmout << "[" << this->n << "]" << flush_now;
    //mmout << "Coefficient: " << this->n << "\n";
    (void) f[0][this->n];
    (void) f[1][this->n];
    direct_transform (0);
    direct_transform (1);
    for (nat level= 0; level < nr_levels; level++) {
      nat   tsh = sh[0] + sh[1];
      nat   cur = this->n + tsh;
      nat   n   = info[level].n;
      if (cur % n != 0) break;
      nat   Cur = cur/n;
      if (Cur < tsh) break;
      //mmout << "  Product at level " << level << "\n";
      bool  last= (level == nr_levels - 1);
      nat   msh = min (sh[0], sh[1]);
      nat   k   = info[level].k;
      nat   tsz = info[level].tsz;
      nat   type= info[level].type;
      const C* h0 = info[level].head[0];
      const C* h1 = info[level].head[1];
      const C* t0 = info[level].tail[0];
      const C* t1 = info[level].tail[1];
      nat   Mod = Cur % k;

      //mmout << "    tsz= " << tsz << ", l= " << 2*n-1 << "\n";
      C* acc= mmx_new<C> (tsz);
      Pol::clear (acc, tsz);
      if (last)
        for (nat i=sh[0]; i<=Cur-sh[1]; i++)
          accumulate (acc, t0 + i*tsz, t1 + (Cur-i)*tsz, tsz);
      else {
        if (Cur < 2 * msh * k) {
          nat start= 1      ; if (Cur > k) start= Cur - k + 1;
          nat end  = Cur - 1; if (Cur > k) end  = k - 1;
          for (nat i=start; i<=end; i++)
            accumulate (acc, h0 + i*tsz, h1 + (Cur-i)*tsz, tsz);
        }
        if (Cur >= msh * k && Mod != 0) {
          //mmout << "    Normal\n";
          if (sh[0] != 0)
            for (nat i=1; i<=Mod; i++)
              accumulate (acc, h0 + i*tsz, t1 + (Mod-i)*tsz, tsz);
          if (sh[1] != 0)
            for (nat i=0; i<=Mod-1; i++)
              accumulate (acc, t0 + i*tsz, h1 + (Mod-i)*tsz, tsz);
        }
      }
      inverse_transform (acc, n, tsz, type);
      //mmout << "    acc= "; Pol::print (mmout, acc, 2*n-1); mmout << "\n";
      Pol::add (this->a + this->n, acc, 2*n - 1);

      if (Cur >= msh * k && Mod == 0 && !last) {
        //mmout << "    Extra " << k << "\n";
        for (nat Mod= 0; Mod<k-1; Mod++) {
          Pol::clear (acc, tsz);
          for (nat i=Mod+1; i<k; i++) {
            if (sh[0] != 0)
              accumulate (acc, h0 + i*tsz, t1 + (k+Mod-i)*tsz, tsz);
            if (sh[1] != 0)
              accumulate (acc, t0 + i*tsz, h1 + (k+Mod-i)*tsz, tsz);
          }
          inverse_transform (acc, n, tsz, type);
          Pol::add (this->a + this->n + Mod*n, acc, 2*n - 1);
        }
      }
      mmx_delete<C> (acc, tsz);
    }
    return this->a [this->n];
  }
};

TMPL static Series
nrelax_mul (const Series& f, const Series& g, nat n) {
  return (Series_rep*) new nrelax_mul_series_rep<C,V> (f, g, n);
}

// Top level interface for fast multiplication

TMPL
class mul_series_rep: public 
  implementation<series_abstractions,V>
    ::template binary_series_rep<mul_op,C,V> {
protected:
  Series prod;
  nat    N;
public:
  inline mul_series_rep (const Series& f, const Series& g):
    implementation<series_abstractions,V>
      ::template binary_series_rep<mul_op,C,V > (f, g),
    prod (nrelax_mul (f, g, 1)), N (1) {}
  C next () { return prod [this->n]; }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (prod, l);
    if (l < 2*N) return;
    while (l >= 2*N) N= 2*N;
    //while (l >= 2*N) N= 4*N;
    prod= nrelax_mul (this->f, this->g, N);
  }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  typedef mul_series_rep<C,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g); }

TMPL static inline Series
ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  typedef mul_series_rep<C,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
    return Series (CF(f));
  return (Series_rep*)
    new Mul_rep (piecewise (f, Series (CF(f)), nf),
                 piecewise (g, Series (CF(g)), ng)); }

}; // implementation<series_multiply,U,series_relaxed>

#undef TMPL
#undef Series
#undef Series_rep
#undef Vector
#undef Series_vector
#undef Series_vector_rep
} // namespace mmx
#endif // __MMX__SERIES_FAST__HPP