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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : p_adic.hpp
* DESCRIPTION: p-adic integers
* 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_P_ADIC_HPP
#define __MMX_P_ADIC_HPP
#include <basix/pair.hpp>
#include <algebramix/p_expansion_modular_int.hpp>
#include <algebramix/p_expansion_modular_integer.hpp>
#include <algebramix/series_int.hpp>
#include <algebramix/series_integer.hpp>
#include <algebramix/series_modular_int.hpp>
#include <algebramix/series_integer.hpp>
#include <algebramix/series_modular_integer.hpp>
#include <algebramix/series_carry.hpp>
#include <algebramix/series_carry_relaxed.hpp>
#include <algebramix/series_carry_blocks.hpp>

namespace mmx {

#define P_adic_variant(M) Series_carry_variant (M)
#define p_adic(M, V) series<M,V>
#define default_p_adic(M) p_adic(M, typename P_adic_variant(M))

#define Series series<M,V>
#define Series_rep series_rep<M,V>
#define C typename M::C
#define TMPL template<typename M,typename V>

/******************************************************************************
* Specific variant for p-adic integer
******************************************************************************/

template<typename V>
struct series_carry_p_adic {};

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

/******************************************************************************
* Naive product optimized for small hardware integer
******************************************************************************/

template<typename V>
struct series_carry_modular_int_naive {};

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

template<typename U,typename W>
struct implementation<series_multiply,U,series_carry_modular_int_naive<W> >:
  public implementation<series_abstractions,U>
{
  typedef implementation<series_abstractions,U> Ser;

template<typename M,typename V>
class mul_series_rep: public Ser::template binary_series_rep<mul_op,M,V> {
  typedef typename M::modulus Modulus;
  typedef integer L;
  typedef unsigned long int N;
  nat nf, ng;
  L carry;
  xnat free_bits;
  
public:
  inline mul_series_rep (const Series& f, const Series& g,
                         nat nf2= Maximal (nat),
                         nat ng2= Maximal (nat)):
    Ser::template binary_series_rep<mul_op,M,V> (f, g),
    nf (nf2), ng (ng2), carry (0) {
    Modulus p= M::get_modulus ();
    xnat p_size= bit_size (* p);
    xnat N_size= 8 * sizeof (N);
    free_bits= N_size <= 2 * p_size ? 0 : (N_size - 2 * p_size); }
  M next () {
    nat n= this->n;
    Modulus p= M::get_modulus ();
    if (n < 8 || free_bits <= 3) {
      C acc= as<C> (rem (carry, * p, carry)), tmp_a, tmp_c;
      for (nat i= n >= ng ? n - ng + 1 : 0; i < min (1 + n, nf); i++) {
        tmp_c= 0;
        mul_mod (tmp_a, (this->f[i]).rep, (this->g[n-i]).rep, p, tmp_c);
        carry += tmp_c; tmp_c= 0;
        add_mod (acc, tmp_a, p, tmp_c);
        carry += tmp_c;
      }
      return acc;
    }
    C acc= as<C> (rem (carry, * p, carry)), tmp_a, tmp_c;
    nat n_max= min (1 + n, nf);
    for (nat i= n >= ng ? n - ng + 1 : 0; i < n_max;) {
      N t0= 0, t1= 0;
      nat i_max= i + min (((nat) 1 << free_bits), (nat) (n_max - i));
      for (; i < i_max; i++)
        t0 += N((this->f[i]).rep) * N((this->g[n-i]).rep);
      tmp_a= rem (t0, * p, t1);
      tmp_c= 0; add_mod (acc, tmp_a, p, tmp_c);
      carry += tmp_c + t1;
    }
    return acc; }
};

template<typename M,typename V> static inline Series
ser_mul (const Series& f, const Series& g) {
  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); }

template<typename M,typename V> 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 (f, g, nf, ng); }

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

/******************************************************************************
* p-th root
******************************************************************************/
#define Pair pair<integer,integer>

struct ser_carry_pth_root_reg_op {
  // For the only computation of the regularized part for p >= 3 only

  template<typename M, typename V> static Series
  binpow_no_tangent_normalized (const Series& x0, const Series& x1,
                                const integer& r, const integer& p) {
    // ((x0 + p * x1)^r - r * p * x1 * x0^(r-1) - x0^r) / p^2
    // x1 has valuation at least 1
    ASSERT (p >= 3, "bug");
    if (r <= 1) return Series (0);

    Series x1_2= lshiftz (square (rshiftz (x1, 1)), 2);
    if (r == 2) return x1_2;

    Series s_p   (coefficients (as<default_p_expansion(M)> (p)));

    if ((r & 1) == 1) {
      integer h= r - 1;
      Series w= binpow_no_tangent_normalized (x0, x1, h, p);
      Series s_h (coefficients (as<default_p_expansion(M)> (h)));
      Series x0_h_1= binpow (x0, h - 1);
      return lshiftz ((x0 + s_p * x1) * rshiftz (w, 1), 1)
        + s_h * x0_h_1 * x1_2;
    }

    integer h= r >> 1; // r is now even
    Series s_h   (coefficients (as<default_p_expansion(M)> (h)));
    Series s_p_2 (coefficients (as<default_p_expansion(M)> (square (p))));
    Series x0_h_1= binpow (x0, h - 1);
    Series x0_h= x0 * x0_h_1;

    Series d= binpow_no_tangent_normalized (x0, x1, h, p);
    Series l= s_h * s_p * x0_h_1 * x1 + x0_h;
    return lshiftz (rshiftz (d, 1) * (s_p_2 * d + l + l), 1)
      + square (s_h * x0_h_1) * x1_2;
  }

  static generic name () { return "pth_root"; }

  template<typename M> static inline M
  op (const M& a, const Pair& p) {
    (void) p; return a; }

  template<typename M, typename V> static inline Series
  op (const Series& a, const Pair& p) {
    typedef implementation<series_pth_root_reg,V> Ser;
    return Ser::unsep_root_reg (a, p); }

  template<typename M, typename V> static inline Series
  op_init (const Series& a, const Pair& p, const M& init) {
    typedef implementation<series_pth_root_reg,V> Ser;
    return Ser::unsep_root_reg (a, p, init); }

  static inline syntactic
  op (const syntactic& a, const syntactic& p) {
    return apply ("pth_root_reg", a, p); }

  static inline nat nr_init () { return 1; }

  template<typename M, typename V> static inline Series
  def (const Series& me, const Series& a, const Pair& p_b0) {
    integer p= first (p_b0); integer b0= second (p_b0);
    integer q= as<integer> (* M::get_modulus ());
    Series ser_p (coefficients (as<default_p_expansion(M)> (p)));
    if (p == 2)
      return (a - me[0] - ser_p
              * (square (Series (me[0]))
                 + lshiftz (square (rshiftz (me, 1)), 2)))
        / (M(1) + square (ser_p) * me[0]);

    Series tmp= binpow (Series (M(b0) + ser_p * me[0]), p)
      - binpow (Series (M(b0)), p);
    tmp= q == p ? rshiftz (tmp, 2): rshiftz (M (q / square (p)) * tmp, 1);
    Series ap= a - tmp;
    Series c= M (b0) + ser_p * me[0];
    return (ap - binpow_no_tangent_normalized (c, me - me[0], p, p))
      / binpow (c, p - 1); }

}; // struct ser_carry_pth_root_reg_op

template<typename U, typename W>
struct implementation<series_pth_root_reg,U,
                      series_carry_p_adic<W> > {

TMPL static inline Series
unsep_root_reg (const Series& a, const Pair& p) {
  return binary_scalar_recursive_series
    <ser_carry_pth_root_reg_op> (a, p); }

TMPL static inline Series
unsep_root_reg (const Series& a, const Pair& p, const M& init) {
  return binary_scalar_recursive_series
    <ser_carry_pth_root_reg_op> (a, p, init); }
};

template<typename U, typename W>
struct implementation<series_pth_root,U,series_carry_p_adic<W> > {

typedef implementation<series_pth_root_reg,U> Reg;

TMPL static inline Series
unsep_root (const Series& a) {
  if (is_exact_zero (a)) return Series (CF(a));
  C p= * M::get_modulus ();
  if (p == 2) {
    ASSERT (a[0] == 1, "positive valuation is not allowed");
    ASSERT (a[1] == 0 && a[2] == 0, "no root");
    return 1 + lshiftz (Reg::unsep_root_reg (rshiftz (a, 3), Pair (p, 1)), 2);
  }
  ASSERT (! is_exact_zero (a[0]), "positive valuation is not allowed");
  M b0= pth_root (a[0]);
  Series ap= a - binpow (Series (b0), p);
  ASSERT (ap[0] == 0 && ap[1] == 0, "no root");
  return b0 + lshiftz (Reg::unsep_root_reg
                       (rshiftz (ap, 2), Pair (p, * b0)), 1); }

}; // implementation<series_pth_root,U,series_carry_p_adic<W> >

// wrapper for block operations
 
template<typename U,typename W,nat s,typename BV,nat t>
struct implementation<series_pth_root_reg,U,
                      series_carry_monoblock<W,s,BV,t> > {

TMPL static inline Series
unsep_root_reg (const Series& a, const Pair& p) {
  return as<Series>
      (binary_scalar_recursive_monoblock_series
       <ser_carry_pth_root_reg_op,M,W,s,BV,t,Pair>
       (as<series<M,W> > (a), p)); }

TMPL static inline Series
unsep_root_reg (const Series& a, const Pair& p, const M& init) {
  return as<Series>
      (binary_scalar_recursive_monoblock_series
       <ser_carry_pth_root_reg_op,M,W,s,BV,t,Pair>
       (as<series<M,W> > (a), p, init)); }
};
#undef Pair

#undef Series
#undef Series_rep
#undef C
#undef TMPL

/******************************************************************************
* Default variants
******************************************************************************/

#define IMPLEMENTATION_HELPER(C)                                        \
  template<typename V,typename W>                                       \
  struct series_carry_variant_helper<modular<modulus<C,V>,W> > {        \
    typedef series_carry_p_adic<                                        \
      series_carry_relaxed<series_carry_naive> > SV; };

IMPLEMENTATION_HELPER(signed char)
IMPLEMENTATION_HELPER(unsigned char)
IMPLEMENTATION_HELPER(short int)
IMPLEMENTATION_HELPER(unsigned short)
IMPLEMENTATION_HELPER(int)
IMPLEMENTATION_HELPER(unsigned int)
IMPLEMENTATION_HELPER(long int)
IMPLEMENTATION_HELPER(unsigned long int)
IMPLEMENTATION_HELPER(long long int)
IMPLEMENTATION_HELPER(unsigned long long int)
#undef IMPLEMENTATION_HELPER

template<typename V,typename W>
struct series_carry_variant_helper<modular<modulus<integer,V>,W> > {
  typedef series_carry_p_adic<
    series_carry_relaxed<series_carry_naive> > SV;
};

} // namespace mmx
#endif // __MMX_P_ADIC_HPP