include/lacunaryx/integer_roots.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : integer_roots.hpp
* DESCRIPTION: Integer roots computation of lacunary polynomials
* COPYRIGHT  : (C) 2014  Bruno Grenet
*******************************************************************************
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see <http://www.gnu.org/licenses/>.
******************************************************************************/

#include <algebramix/polynomial_integer.hpp>
#include <factorix/irreducible_polynomial_rational.hpp>
#include <multimix/sparse_polynomial.hpp>
#include <multimix/multivariate.hpp>
#include <lacunaryx/lpolynomial.hpp>

namespace mmx {

#define Monomial monomial<vector<E> >
#define Sparse_polynomial sparse_polynomial<C, Monomial, V>
#define Dense_polynomial polynomial<C>
#define Root pair<C,E>


template<typename E, typename C, typename V>
Sparse_polynomial sparse_derive (const Sparse_polynomial& p) {
  E d= deg(get_monomial(p,0));
  vector<Monomial> monomials= fill<Monomial> (vec(0),N(p)-1);
  vector<C> coefficients= fill<C> (0,N(p)-1);

  for (nat i=1; i<N(p); i++) {
    pair<C,Monomial> term= p[i];
    monomials[i-1]= cdr(term)/Monomial(vec(d+1));
    coefficients[i-1]= car(term)*(as<C> (deg (cdr (term)) - d));
  }
  
  return Sparse_polynomial(coefficients,monomials,1);
}

// Partition the input polynomial into a vector
// Integer roots of q equal the common integer roots of the 
//  polynomials in the vector
template<typename E, typename C, typename V>
vector<Sparse_polynomial> partition (const Sparse_polynomial& q) {
    C max_coeff= abs(get_coefficient(q,0));
    E prev_exponent= deg(get_monomial(q,0));
    E curr_exponent;
    vector<Sparse_polynomial> polys; 
    Sparse_polynomial tmp;
    nat index_min=0;

    for (nat i=1; i<N(q); i++) {
      
      curr_exponent= deg(get_monomial(q,i));

      if ( curr_exponent - prev_exponent > bit_size( max_coeff ) ) {
        tmp= range (q,index_min,i);
        tmp /= Monomial (vec (val(tmp))); 
        polys << tmp; 
        prev_exponent = curr_exponent;
        max_coeff= abs(get_coefficient(q,i));
        index_min= i;
      } else {
        prev_exponent= curr_exponent;
        max_coeff= max(abs(get_coefficient(q,i)),max_coeff);
      }
    }
    tmp= range(q, index_min, N(q));
    tmp /= Monomial (vec (val(tmp))); 
    polys << tmp; 

    return polys;
}


template<typename E, typename C, typename V>
vector<Root> 
large_integer_roots (const Sparse_polynomial &q) {
    vector<Sparse_polynomial> polys= partition(q);
    vector<Dense_polynomial> vp (as<Dense_polynomial> (0),N(polys));

    for ( nat i=0; i < N( polys ); i++ ) {
      Sparse_polynomial p= polys[i];
      vector<C> c(0, as<nat> (deg(p)+1));
      for ( nat j=0; j < N( p ); j++) 
        c[as<nat> (deg(get_monomial(p,j)))]= get_coefficient(p,j);
      vp[i]= Dense_polynomial (c);
    }

    Dense_polynomial g= primitive_part ( big<gcd_op> ( vp ));
    polynomial<rational> g_rat= as< polynomial<rational> > (g);
    vector <irreducible_factor <polynomial <rational> > > fact=
      bounded_irreducible_factorization (g_rat, 2);
      //irreducible_factorization (g_rat);

    vector<Root> r;

    for (nat i=0; i < N(fact); i++)
    {
      polynomial<rational> f=factor(fact[i]);
      if ( deg(f) == 1 && is_integer(f[0]/f[1]) && abs (f[0]/f[1]) != 1 ) 
        r << Root (as<C> (numerator (-f[0]/f[1])), fact[i].m);
    }

    return r;
}

template<typename E, typename C, typename V>
vector<Root> special_integer_roots (const Sparse_polynomial &p ){
  
  vector<Root> r;
  C s;
  
  // Root O
  E v= val(p);
  if ( v > 0 ) r << Root( 0, v );

  // Root 1 & -1
  Sparse_polynomial q= copy(p); 
  E mul1=0;
  E mul2=0;

  for (nat i=0; i<N(p); i++){
    
    // Root 1
    if (mul1 == i) {
      s=0;
      for (nat j=0; j<N(q); j++) s+=get_coefficient (q,j);
      if ( s == 0 ) mul1=i+1;
    }

    // Root -1
    if (mul2 == i) {
      s=0;
      for (nat j=0; j<N(q); j++) 
        s+= get_coefficient (q,j) * binpow(-1, deg (get_monomial (q,j)) % 2 );
      if ( s == 0 ) mul2=i+1;
    }

    if ( mul1 < i+1 && mul2 < i+1 ) break; // Useless to go further

    q= sparse_derive(q);
  }

  if (mul1 > 0) r << Root (1,mul1);
  if (mul2 > 0) r << Root (-1,mul2);

  return r;
}

template<typename E, typename C, typename V>
vector<Root> integer_roots ( const Sparse_polynomial &q ){
  vector<Root> r= large_integer_roots (q);
  r << special_integer_roots (q);
  return r;
}

// Wrapper for multivariate
template<typename E, typename C, typename V>
vector<generic> integer_roots ( const multivariate<Sparse_polynomial> &p ){
  vector<generic> ret;
  vector<Root> roots= integer_roots (p.rep);
  for (nat i=0; i<N(roots); i++) 
    ret << as<generic> (vec (car(roots[i]), as<integer> (cdr(roots[i]))));
  return ret; //integer_roots (p.rep);
}

// Wrapper for lpolynomial
template<typename E, typename C, typename W>
vector<generic> integer_roots ( const lpolynomial<C,E,W> &lp ){
  typedef typename Sparse_polynomial_variant(C) V;
  // set as Sparse_polynomial
  vector<E> w (entries (*lp));
  vector<C> c;
  vector<Monomial> m;
  for (nat i= 0; i < N(w); i++)
    if (lp[w[i]] != 0) {
      c << lp[w[i]];
      m << Monomial (vec(w[i]));
    }
  Sparse_polynomial p (c, m, 1);

  // Roots
  vector<Root> roots= integer_roots (p);
  vector<generic> ret;
  for (nat i=0; i<N(roots); i++)
    ret << as<generic> (vec (car(roots[i]), as<integer> (cdr(roots[i]))));
  return ret;
}

#undef Monomial
#undef Sparse_polynomial
#undef Dense_polynomial
#undef Root


} // mmx namespace

---