include/algebramix/polynomial_series.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_series.hpp
* DESCRIPTION: Polynomials over power series numbers
* COPYRIGHT  : (C) 2009  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_POLYNOMIAL_SERIES_HPP
#define __MMX_POLYNOMIAL_SERIES_HPP
#include <algebramix/series.hpp>
#include <algebramix/series_vector.hpp>
#include <algebramix/polynomial.hpp>
#include <algebramix/polynomial_ring_dicho.hpp>

namespace mmx {
#define TMPL   template<typename C, typename U, typename V, typename W>
#define Series series<C,V>
#define Series_vector     series <vector<C,W> >
#define Series_vector_rep series_rep <vector<C,W> >
#define Recursive_series_vector_rep recursive_series_rep<vector<C,W> >
#define Polynomial polynomial<Series,U>

/******************************************************************************
* Variants
******************************************************************************/

template<typename V>
struct polynomial_series: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_series<typename V::Positive> Positive;
  typedef polynomial_series<typename V::No_simd> No_simd;
  typedef polynomial_series<typename V::No_thread> No_thread;
  typedef polynomial_series<typename V::No_scaled> No_scaled;
};

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

DEFINE_VARIANT (polynomial_series_dicho,
                  polynomial_series<
                    polynomial_dicho<
                      polynomial_naive> >)

template<typename C, typename V>
struct polynomial_variant_helper<Series> {
  typedef polynomial_series_dicho PV;
};

/******************************************************************************
* Modular inversion
******************************************************************************/

template<typename X, typename BV>
struct implementation<polynomial_gcd,X,polynomial_series<BV> >:
  public implementation<polynomial_gcd,X,BV> {

typedef implementation<polynomial_gcd,X,BV> Pol;
  
TMPL
class inv_mod_polynomial_series_rep: public Recursive_series_vector_rep {
  
public:
  const Polynomial a, b; // inverse of a mod b

  inv_mod_polynomial_series_rep (const Polynomial& a2,
                                 const Polynomial& b2) : 
    Recursive_series_vector_rep (CF (b2)), a (a2), b (b2) {}
  virtual void Increase_order (nat l) {
    Recursive_series_vector_rep::Increase_order (l); }
  syntactic expression (const syntactic& z) const {
    return flatten (1) / flatten (a); }
  Series_vector initialize () {
    VERIFY (N(b) > 0, "unexpected zero polynomial");
    polynomial<C> a0, b0;
    vector<C> coeff (C (), N (a));
    for (nat i = 0; i < N (a); i++)
      coeff[i] = a[i][0];
    a0 = polynomial<C> (coeff);
    coeff = vector<C> (C (), N (b));
    for (nat i = 0; i < N (b); i++)
      coeff[i] = b[i][0];
    b0 = polynomial<C> (coeff);
    polynomial<C> u= invert_modulo (a0, b0);
    if (u == 0) return Series_vector (); 
    nat n= N(b) - 1;
    Polynomial c (rem (a * as<Polynomial> (u) - 1, b));
    vector<Series> w= coefficients (c, 0, n);
    w= as_vector<C,V,W> (rshiftz (as_series (w),1),n);
    c= rem (Polynomial (w)
            * Polynomial (as_vector<C,V,W> (this->me (), n)), b);
    w= coefficients (c, 0, n);
    w= as_vector<C,V,W> (lshiftz (as_series (w),1), n);
    c= as<Polynomial> (u) - Polynomial (w);
    return as_series (coefficients (c, 0, n)); }
};

TMPL static inline Series_vector
inv_mod_polynomial_series (const Polynomial& a, const Polynomial &b) {
  Series_vector_rep* rep= new inv_mod_polynomial_series_rep<C,U,V,W> (a, b);
  return recursive (Series_vector (rep)); }

template<typename C, typename U, typename V> static Polynomial
invert_mod (const Polynomial& a, const Polynomial& b) {
  typedef typename Vector_variant(C) W;
  Series_vector v= inv_mod_polynomial_series<C,U,V,W> (a, b);
  if (b == 0) {
    if (deg (a) != 0) return 0;
    return 1 / a[0];
  }
  vector<Series> tmp= as_vector (v, N(b) - 1);
  return Polynomial (tmp); }

}; // implementation<polynomial_gcd,X,polynomial_series<BV> >

#undef TMPL
#undef Series
#undef Series_vector
#undef Series_vector_rep
#undef Recursive_series_vector_rep
#undef Polynomial
} // namespace mmx
#endif // __MMX_POLYNOMIAL_SERIES_HPP

---