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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_carry_blocks.hpp
* DESCRIPTION: blocks 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_BLOCKS_HPP
#define __MMX_SERIES_CARRY_BLOCKS_HPP
#include <algebramix/series_carry.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>
#define Monoblock_transformer series_carry_monoblock_transformer<M,V,s,BV>

/******************************************************************************
* Monoblock variant
******************************************************************************/

// V fallback variant
// s is the size of the blocks
// BV block variant
// t is the threshold
template<typename V,
         nat s,
         typename BV,
         nat t=s>
struct series_carry_monoblock {};

template<typename V,typename BV,typename U,typename W,nat s,nat t>
struct implementation<U,W,series_carry_monoblock<V,s,BV,t> >:
  public implementation<U,W,V> {};

/******************************************************************************
* Relaxed blocks variant
******************************************************************************/

template<typename V,
         nat s,
         typename BV,
         nat t=s>
struct series_carry_blocks {};

template<typename V,typename BV,typename U,typename W,nat s,nat t>
struct implementation<U,W,series_carry_blocks<V,s,BV,t> >:
  public implementation<U,W,V> {};

/******************************************************************************
* Blocks transformer
******************************************************************************/

#define Monoblock_modular modular<modulus<Lift_type(M)>,                \
                                  modular_global_series_carry_monoblock \
                                  <M,s,BV> >

#define Monoblock_series_rep series_rep<Monoblock_modular,BV>
#define Monoblock_series     series    <Monoblock_modular,BV>

template<typename M,nat s,typename BV>
struct modular_global_series_carry_monoblock {
  template<typename P>
  class modulus_storage {
    static inline P& dyn_modulus () {
      static P modulus = P ();
      return modulus; }
  public:
    static inline void set_modulus (const P& p) { dyn_modulus () = p; }
    static inline P get_modulus () { return dyn_modulus (); }
  };
};

template<typename M, typename V, nat s, typename BV>
class series_carry_monoblock_transformer {

class to_monoblock_series_rep: public Monoblock_series_rep {
public:
  typedef Lift_type(M) L;
  typedef typename Base_transformer_unsigned(L,C) Baser;
   
protected:
  Series f;
  C* buf;
  Baser* baser;

public:
  to_monoblock_series_rep (const Series& f2):
    Monoblock_series_rep (format<Monoblock_modular > (/*FIXME*/)), f(f2) {
    Monoblock_modular::set_modulus (binpow (L(* M::get_modulus ()), s));
    buf= mmx_new<C> (s);
    baser= mmx_new<Baser> (1, * M::get_modulus ()); }

  ~to_monoblock_series_rep () {
    mmx_delete<Baser> (baser, 1);
    mmx_delete<C> (buf, s); }

  syntactic expression (const syntactic& z) const {
    return apply ("to_monoblock", flatten (f, z), flatten (s)); }

  void Increase_order (nat l) {
    Monoblock_series_rep::Increase_order (l);
    increase_order (f, l * s); }

  Monoblock_modular next () {
    nat start = s * this->n;
    for (nat i= 0; i < s; i++) buf[i]= * f[start + i];
    L e; inverse_base (e, buf, s, * baser);
    return Monoblock_modular (e, true); }
};

class from_monoblock_series_rep: public Series_rep {
public:
  typedef Lift_type(M) L;
  typedef typename Base_transformer_unsigned(L,C) Baser;
  
protected:
  Monoblock_series F;
  C* buf; // p-adic decomposition of current term of F
  Baser* baser;
public:

  from_monoblock_series_rep (const Monoblock_series& F2):
    Series_rep (format<M> (/*FIXME*/)), F(F2) {
    buf= mmx_new<C> (s);
    baser= mmx_new<Baser> (1, * M::get_modulus ()); }

  ~from_monoblock_series_rep () {
    mmx_delete<C> (buf, s);
    mmx_delete<Baser> (baser, 1); }

  syntactic expression (const syntactic& z) const {
    return apply ("from_monoblock", flatten (F, z), flatten (s)); }

  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (F, (l + s - 1) / s); }
  
  M next () {
    nat q= this->n / s, r= this->n % s;
    if (r == 0) {
      nat m= direct_base (buf, s, * F[q], * baser);
      for (nat i= m; i < s; i++) buf[i]= 0;
    }
    return buf[r];
  }
};

public:

inline series_carry_monoblock_transformer () {}

inline Monoblock_series
to_monoblock (const Series& f) const {
  if (is_exact_zero (f))
    return Monoblock_series (Monoblock_modular (0) /*FIXME*/);
  return (Monoblock_series_rep*) new to_monoblock_series_rep (f);
}

inline Series
from_monoblock (const Monoblock_series& f) const {
  if (is_exact_zero (f))
    return Series (C(0) /*FIXME*/);
  return (Series_rep*) new from_monoblock_series_rep (f);
}

};

// Top level conversions

template<typename M, typename V, nat s, typename BV> Monoblock_series
to_monoblock (const Series& f,
              const Monoblock_transformer& blocker) {
  return blocker.to_monoblock (f);
}

template<typename M, typename V, nat s, typename BV> Series
from_monoblock (const Monoblock_series& f,
                const Monoblock_transformer& blocker) {
  return blocker.from_monoblock (f);
}

/******************************************************************************
* Monoblock abstractions
******************************************************************************/

template<typename Op, typename M, typename V, nat s, typename BV, nat t>
class binary_monoblock_series_rep: public Series_rep {
  Monoblock_transformer blocker;
  Series f, g;
  bool h1_init;
  Series h0, h1;  

public:
  binary_monoblock_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f(f2), g(g2), h1_init (false) {
    h0= Op::op (f, g); }
  syntactic expression (const syntactic& z) const {
    return Op::op (flatten (f, z), flatten (g, z)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    if (l <= t)
      increase_order (h0, l);
    else {
      if (! h1_init) {
        h1= from_monoblock (Op::op (to_monoblock (f, blocker),
                                    to_monoblock (g, blocker)),
                             blocker);
        h1_init= true;
      }
      increase_order (h1, l); } }
  M next () {
    return this->n < t ? h0[this->n] : h1[this->n]; }
};

template<typename Op, typename M, typename V, nat s, typename BV, nat t>
inline Series
binary_monoblock_series (const Series& f, const Series& g) {
  typedef binary_monoblock_series_rep<Op,M,V,s,BV,t> Op_rep;
  return (Series_rep*) new Op_rep (f, g); }

/******************************************************************************
* Truncate and multiply
******************************************************************************/

template<typename M, typename V, nat s, typename BV, nat t>
class truncate_mul_monoblock_series_rep: public Series_rep {
  Monoblock_transformer blocker;
  Series f, g;
  nat nf, ng;
  bool h1_init;
  Series h0, h1;  

public:
  truncate_mul_monoblock_series_rep (const Series& f2, const Series& g2,
                                     nat nf2, nat ng2):
    Series_rep (CF(f2)), f(f2), g(g2), nf (nf2), ng (ng2), h1_init (false) {
    h0= truncate_mul (f, g, nf, ng); }
  syntactic expression (const syntactic& z) const {
    return apply ("truncate_mul", flatten (f, z), flatten (g, z),
                  flatten (nf), flatten (ng)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    if (l <= t)
      increase_order (h0, l);
    else {
      if (! h1_init) {
        Series f0= piecewise (f, Series (CF(f)), nf);
        Series g0= piecewise (g, Series (CF(g)), ng);
        h1= from_monoblock (truncate_mul (to_monoblock (f0, blocker),
                                          to_monoblock (g0, blocker),
                                          (nf + s - 1) / s,
                                          (ng + s - 1) / s),
                            blocker);
        h1_init= true;
      }
      increase_order (h1, l); } }
  M next () {
    return this->n < t ? h0[this->n] : h1[this->n]; }
};

template<typename M, typename V, nat s, typename BV, nat t>
inline Series
truncate_mul_monoblock_series (const Series& f, const Series& g,
                               nat nf, nat ng) {
  typedef truncate_mul_monoblock_series_rep<M,V,s,BV,t> Mul_rep;
  return (Series_rep*) new Mul_rep (f, g, nf, ng); }

/******************************************************************************
* Recursive monoblock abstractions -- require monoblock-recursivity!
******************************************************************************/

template<typename Op, typename M, typename V, nat s,
         typename BV, nat t, typename X>
class binary_scalar_recursive_monoblock_series_rep: public Series_rep {
  Monoblock_transformer blocker;
  Series f;
  X x;
  bool h1_init;
  Series h0, h1;  

public:
  binary_scalar_recursive_monoblock_series_rep (const Series& f2, const X& x2):
    Series_rep (CF(f2)), f(f2), x(x2), h1_init (false) {
    h0= Op::op (f, x); }
  binary_scalar_recursive_monoblock_series_rep (const Series& f2, const X& x2,
                                                const M& init):
    Series_rep (CF(f2)), f(f2), x(x2), h1_init (false) {
    h0= Op::op (f, x, init); }
  syntactic expression (const syntactic& z) const {
    return apply (Op::name (), flatten (f, z), flatten (x)); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    if (l <= t)
      increase_order (h0, l);
    else {
      if (! h1_init) {
        Monoblock_modular init= as<Lift_type(M)> (truncate (h0, s));
        h1= from_monoblock (Op::op_init (to_monoblock (f, blocker), x, init),
                            blocker);
        h1_init= true;
      }
      increase_order (h1, l); } }
  M next () {
    return this->n < t ? h0[this->n] : h1[this->n]; }
};

template<typename Op, typename M, typename V, nat s,
         typename BV, nat t, typename X> inline Series
binary_scalar_recursive_monoblock_series (const Series& f, const X& x) {
  typedef binary_scalar_recursive_monoblock_series_rep<Op,M,V,s,BV,t,X> Op_rep;
  return (Series_rep*) new Op_rep (f, x); }

#undef Monoblock_modular
#undef Monoblock_series_rep
#undef Monoblock_series

/******************************************************************************
* Monoblock multiplication
******************************************************************************/

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_multiply,U,series_carry_monoblock<W,s,BV,t> >:
  public implementation<series_abstractions,U>
{

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return as<Series> (binary_monoblock_series<mul_op,M,W,s,BV,t>
                     (as<series<M,W> > (f), as<series<M,W> > (g))); }

TMPL static inline Series
ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
    return Series (CF(f));
  return as<Series> (binary_monoblock_series<mul_op,M,W,s,BV,t>
                     (as<series<M,W> > (f), as<series<M,W> > (g))); }
}; // implementation<series_multiply,U,series_carry_monoblock<W> >:

/******************************************************************************
* Monoblock Division
******************************************************************************/

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_divide,U,series_carry_monoblock<W,s,BV,t> >:
  public implementation<series_multiply,U>
{

TMPL static inline Series
ser_rdiv_sc (const Series& f, const M& c) {
  typedef implementation<series_divide,W> Ser;
  return Ser::ser_rdiv_sc (f, c); }

TMPL static inline Series
ser_div (const Series& f, const Series& g) {
  if (is_exact_zero (f))
    return Series (CF(f));
  return as<Series> (binary_monoblock_series<div_op,M,W,s,BV,t>
                     (as<series<M,W> > (f), as<series<M,W> > (g))); }

TMPL static inline Series
ser_quo (const Series& f, const M& c) {
  return ser_div_sc (f, c); }
  
TMPL static inline Series
ser_quo (const Series& f, const Series& g) {
  return ser_div (f, g); }

TMPL static inline Series
ser_rrem_sc (const Series& f, const M& c) {
  return Series (CF(f)); }

}; // implementation<series_divide,U,series_carry_monoblock<W> >

/******************************************************************************
* separable roots
******************************************************************************/

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_separable_root,U,
                      series_carry_monoblock<W,s,BV,t> >:
    public implementation<series_multiply,U>
{

TMPL static inline Series
sep_root (const Series& f, nat r) {
  ASSERT (M(r) != 0, "not implemented");
  return as<Series>
    (binary_scalar_recursive_monoblock_series
     <ser_separable_root_op,M,W,s,BV,t,nat> (as<series<M,W> > (f), r)); }

TMPL static inline Series
sep_root_init (const Series& f, nat r, const M& init) {
  ASSERT (M(r) != 0, "not implemented");
  return as<Series>
    (binary_scalar_recursive_monoblock_series
     <ser_separable_root_op,M,W,s,BV,t,nat> (as<series<M,W> > (f), r, init)); }
};

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

template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_multiply,U,series_carry_blocks<W,s,BV,t> >:
  public implementation<series_abstractions,U>
{

TMPL
class mul_series_rep: public Series_rep {
  series_carry_monoblock_transformer<M,W,s,BV> blocker;
  Series f, g, h;
  
public:
  mul_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f(f2), g(g2)
  {
    VERIFY (t >= s-1, "wrong threshold");
    series<M,W> f0 = as<series<M,W> > (f);
    series<M,W> g0 = as<series<M,W> > (g);
    series<M,W> f1 = rshiftz (f0, (int) t);
    series<M,W> g1 = rshiftz (g0, (int) t);
    h= as<Series> (truncate_mul (f0, g0, t, t)
                   + lshiftz (truncate_mul (f0, g1, t, Maximal (nat)) +
                              truncate_mul (f1, g0, Maximal (nat), t),
                              (int) t)
                   + lshiftz (from_monoblock (to_monoblock (f1, blocker) *
                                              to_monoblock (g1, blocker),
                                              blocker),((int) t) << 1)); }
  syntactic expression (const syntactic& z) const {
    return flatten (f, z) * flatten (g, z); }
  virtual void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (h, l); }
  M next () {
    return h[this->n]; }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return (Series_rep*) new mul_series_rep<M,V> (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_blocks<W> >:

#undef Series
#undef Series_rep
#undef C
#undef Modulus
#undef TMPL
#undef Monoblock_transformer
} // namespace mmx
#endif // __MMX_SERIES_CARRY_BLOCKS_HPP