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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_carry_relaxed.hpp
* DESCRIPTION: relaxed p-adic numbers
* COPYRIGHT  : (C) 2009  Jeremy Berthomieu and 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_SERIES_CARRY_RELAXED_HPP
#define __MMX_SERIES_CARRY_RELAXED_HPP
#include <algebramix/series_relaxed.hpp>
#include <algebramix/series_carry_naive.hpp>

namespace mmx {
#define Series     series<M,V>
#define Series_rep series_rep<M,V>
#define C       typename M::C
#define Modulus typename M::modulus
#define TMPL    template<typename M, typename V>

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

template<typename V>
struct series_carry_relaxed {};

template<typename U,typename V,typename W>
struct implementation<V,W,series_carry_relaxed<U> >:
  public implementation<V,W,U> {};

/******************************************************************************
* Relaxed multiplication
******************************************************************************/

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

public:
TMPL
class mul_series_rep: public Series_rep {
public:
  typedef Lift_type(M) L;

protected:
  Series f, g;

  // NOTE: The following fields keep track of large segments of zeros
  // in the coefficients. This makes the multiplication algorithm linear
  // in the case when we multiply by a polynomial. This is for instance
  // interesting when computing coefficients of holonomic functions.
  int lnz_f;    // last non-zero coefficient in f
  int lnz_g;    // last non-zero coefficient in g
  nat zeros_f;  // bit 2^p set if f[2^p-1 ... 2^{p+1}-1] = 0
  nat zeros_g;  // bit 2^p set if g[2^p-1 ... 2^{p+1}-1] = 0

  vector<L> b_f, d_f, b_g, d_g; // binary representations for the segments
  // b_f[q] is the last segment of f of size 2^q
  // d_f[q] is the "diagonal segment" of size 2^q
  vector<L> powers_of_p;  // p^(2^i)
  vector<L> carry; // carry[q] is the carry in size 2^q

L get_power_of_p (const Modulus& p, nat i) {
  // lazy access to p^(2^i), for i >= 0
  // with its own memory allocation
  static const L zero (0);
  while (i >= N (powers_of_p))
    powers_of_p << vector<L> (zero, max ((nat) 1, N (powers_of_p)));
  if (powers_of_p [i] == zero) {
    if (i == 0) powers_of_p[0]= L(* p);
    else powers_of_p[i]= square (get_power_of_p (p, i - 1));
  }  
  return powers_of_p [i];
}

public:
  mul_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f (f2), g (g2), lnz_f (-1), lnz_g (-1),
    zeros_f (0), zeros_g (0) {}
  ~mul_series_rep () {}

  syntactic expression (const syntactic& z) const {
    return flatten (f, z) * flatten (g, z); }

  void Set_order (nat l2) {
    Series_rep::Set_order (l2);
    nat m= log_2 (next_power_of_two (l2 + 1));
    if (N(b_f) < m) {
      vector<L> tmp (L(0), m - N(b_f));
      b_f << tmp; d_f << tmp;
      b_g << tmp; d_g << tmp;
      carry << tmp; } }

  void Increase_order (nat l) {
    Set_order (l);
    increase_order (f, l);
    increase_order (g, l); }

  M next () {
    //mmout << "--> next " << this->n << "\n";
    Modulus p= M::get_modulus ();
    nat k= this->n;
    if (exact_neq (f[k], M(0))) lnz_f= k;
    if (exact_neq (g[k], M(0))) lnz_g= k;
    if (k == 0) {
      C r= 0;
      if (lnz_f < 0) zeros_f |= 1;
      if (lnz_g < 0) zeros_g |= 1;
      if ((zeros_f&1) == 0 && (zeros_g&1) == 0) {
        C q= 0;
        mul_mod (r, * f[0], * g[0], p, q);
        carry[0]= as<L> (q);
      }
      b_f[0]= d_f[0]= lift (f[0]);
      b_g[0]= d_g[0]= lift (g[0]);
      return M (r, true);
    }
    else {
      k = 2 * (this->n + 2);
      L t= 0, e_f= lift (f[this->n]), e_g= lift (g[this->n]);
      nat q= -1;
      while ((k&1) == 0 && k != 2) {
        q++; k= k >> 1;
        //mmout << "k= " << k << "\n";
        if (q > 0) {
          e_f= b_f[q-1] + get_power_of_p (p, q-1) * e_f;
          e_g= b_g[q-1] + get_power_of_p (p, q-1) * e_g;
        }
        if (k == 2) {
          if (lnz_f < ((int) (1<<q)-1)) zeros_f |= (1<<q);
          if (lnz_g < ((int) (1<<q)-1)) zeros_g |= (1<<q);
          d_f[q]= e_f; d_g[q]= e_g;
        }
        t += carry[q];
        if (q == 0 ||
            ((zeros_f & (1<<q)) == 0 && lnz_g >= ((int) (((k-1)<<q)-1))))
          t += d_f[q] * e_g;
        if (k == 2) break;
        if (q == 0 ||
            ((zeros_g & (1<<q)) == 0 && lnz_f >= ((int) (((k-1)<<q)-1)))) {
          t += e_f * d_g[q];
        }
      }
      b_f[q]= e_f; b_g[q]= e_g;
      while (q != (nat) -1) {
        t = rem (t, get_power_of_p (p, q), carry[q]);
        q--;
      }
      return M (as<C> (t), true); } }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  //mmout << "f= " << f << "\n";
  //mmout << "g= " << g << "\n";
  typedef mul_series_rep<M,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<M,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_carry_relaxed>:
   
#undef Series
#undef Series_rep
#undef C
#undef Modulus
#undef TMPL

} // namespace mmx
#endif // __MMX_SERIES_CARRY_RELAXED_HPP