include/lacunaryx/linear_factors_univariate.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : linear_factors_univariate.hpp
* DESCRIPTION: Linear factors of univariate sparse 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 <multimix/sparse_polynomial.hpp>
#include <multimix/multivariate.hpp>
#include <factorix/irreducible_polynomial_integer.hpp>
#include <factorix/irreducible_polynomial_rational.hpp>

namespace mmx {

#define Monomial              monomial< vector<E> >
#define Sparse_polynomial     sparse_polynomial<C, Monomial, V>
#define Univariate_polynomial polynomial<C>
#define Univariate_factor     irreducible_factor<polynomial <rational> >
#define Univariate_factors    vector<Univariate_factor>
#define Sparse_factor         irreducible_factor<Sparse_polynomial>
#define Sparse_factors        vector<Sparse_factor>


namespace lacunaryx {
  template<typename C, typename V, typename M, typename W>
  void set_as_polynomial (polynomial<C,V> &p, const sparse_polynomial<C,M,W> &q) {
    //ASSERT (dim(q) == 1, "wrong dimension");
    nat n= as<nat> (deg(q)+1);
    nat l= aligned_size<C,V>(n);
    C* c= mmx_new<C>(l,promote(0,CF(q)));
    for (nat i=0; i<N(q); i++) 
      c[as<nat> (deg(get_monomial(q,i)))]= get_coefficient(q,i);
    p= polynomial<C,V>(c,n,l,CF(q));
  }

  template<typename C, typename M, typename W>
  polynomial<C> as_polynomial (const sparse_polynomial<C,M,W> &q) {
    polynomial<C> p (CF(q));
    set_as_polynomial(p,q);
    return p;
  }

  // multiplicity
  template<typename C>
  nat multiplicity (const Univariate_polynomial& f, 
                    const Univariate_polynomial& p) {
    nat m=0;
    Univariate_polynomial q= copy(p);
    while (divides (f,q)) {
      m+=1;
      q=derive (q);
    }
    return m;
  } //

  template<typename C, typename E, typename V>
  Sparse_polynomial 
  as_sparse_polynomial (const polynomial<rational>& q) {
    vector<C> coefficients;
    vector<Monomial> monomials;
    C lcm_den=1;
    for (nat i=0; i<N(q); i++) {
      rational c= q[i];
      if (c != 0) lcm_den= lcm(lcm_den,as<C> (denominator (c)));
    }
    for (nat i=0; i<N(q); i++) {
      rational c= q[i];
      if (c != 0) {
        C num= as<C> (numerator (c));
        C den= as<C> (denominator (c));
        coefficients << num*(lcm_den/den);
        monomials << Monomial (vec(i));
      }
    }
    return Sparse_polynomial (coefficients, monomials, 1);
  }

  template<typename C, typename E, typename V>
  Sparse_factor as_sparse_factor (const Univariate_factor& f) {
    return Sparse_factor (as_sparse_polynomial<C,E,V> (f.f), f.m);
  }

  template<typename C, typename E, typename V>
  Sparse_factors 
  intersect_factors (const Sparse_factors& f1, 
                     const Sparse_factors& f2) {
    Sparse_factors f;
    for (nat i=0; i<N(f1); i++) {
      for (nat j=0; j<N(f2); j++) {
        if ((f1[i].f == f2[j].f) || (f1[i].f == -f2[j].f)) {
          if ( f1[i].m < f2[j].m ) f << f1[i]; 
          else f << f2[j];
          break;
        }
      }
    }
    return f;
  }
  
  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);
  }
  
  /* tentative: sparse_derive based on derive
  template<typename E, typename C>
  Sparse_polynomial new_sparse_derive (const Sparse_polynomial& p) {
    
    Monomial m= get_monomial(p,0);
    Sparse_polynomial p_derive= range(p,1,N(p))/m;
    p_derive= derive (p,0);
  }
  */
  
} // lacunaryx namespace

using namespace lacunaryx;

// Computation of the gap (cf [Lenstra'99])
template<typename E, typename C, typename V>
C gap_univariate (const Sparse_polynomial& p) {
  C g=0;
  Sparse_polynomial q= copy(p);

  for (nat i=0; i<N(p); i++) {
    C max_coeff= abs (get_coefficient (q,0));
    for (nat j=1; j<N(q); j++) {
      C curr_coeff = abs (get_coefficient (q,j));
      if (curr_coeff > max_coeff) max_coeff= curr_coeff;
    }
    C length= as<C> (N(q)-1);
    g= max (g, as<C> (bit_size (length*max_coeff)));
    q= sparse_derive (q);
  }
  
  return g;
}

template<typename E, typename C, typename V>
vector<Sparse_polynomial> 
partition_univariate (const Sparse_polynomial& q, E g) {
    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 > g ) {
        tmp= range (q,index_min,i);
        tmp /= Monomial (vec (val(tmp))); 
        polys << tmp; 
        prev_exponent = curr_exponent;
        index_min= i;
      } 
      else prev_exponent= curr_exponent; 
    }
    tmp= range(q, index_min, N(q));
    tmp /= Monomial (vec (val(tmp))); 
    polys << tmp; 

    return polys;
}

template<typename C, typename E, typename V>
nat multiplicity (const Sparse_polynomial& f, const Sparse_polynomial& p) {
  Univariate_polynomial f_dense= as_polynomial (f);
  if (f_dense[0] == 0) return as<nat> (deg (get_monomial (p,0)));
  if (f_dense[0] == -f_dense[1]) {
    Sparse_polynomial q= copy(p);
    for (nat i=0; i<N(p); i++) {
      C s=0;
      for (nat j=0; j<N(q); j++) s+= get_coefficient (q,j); 
      if ( s != 0 ) return i; 
      q= sparse_derive (q);
    }
    return N(p)-1; // should never occur
  }
  if (f_dense[0] == f_dense[1]) {
    Sparse_polynomial q= copy(p);
    for (nat i=0; i<N(p); i++) {
      C 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 ) return i; 
      q= sparse_derive (q);
    }
    return N(p)-1; // should never occur
  }
  E gap= as<E> (gap_univariate (p));
  vector<Sparse_polynomial> polys= partition_univariate(p, gap);
  nat m= multiplicity (f_dense, as_polynomial (polys[0]));
  for ( nat i=0; i < N( polys ); i++ ) {
    nat mm= multiplicity (f_dense, as_polynomial (polys[i]));
    if (mm == 0) return mm;
    if (mm < m) m=mm;
  }
  return m;
}

template<typename E, typename C, typename V>
Sparse_factors
non_cyclotomic_factors (const Sparse_polynomial &q, E gap) {
    vector<Sparse_polynomial> polys= partition_univariate(q, gap);
    vector<Univariate_polynomial> vp (as <Univariate_polynomial> (0),N(polys));
    
    for ( nat i=0; i < N( polys ); i++ ) vp[i]= as_polynomial (polys[i]);

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

    Sparse_factors sparse_fact;
    for (nat i=0; i < N(fact); i++)
    {
      polynomial<rational> f=factor(fact[i]);
      if ( deg(f) == 1 && abs(f[1]*f[0]) != 1 ) {
        sparse_fact <<  as_sparse_factor<C,E,V> (fact[i]);
      }
    }

    return sparse_fact;
}

template<typename E, typename C, typename V>
Sparse_factors
cyclotomic_factors (const Sparse_polynomial &p ){
  
  Sparse_factors f;
  Sparse_polynomial x ((C)1, Monomial (vec(1)));
  Sparse_polynomial q= copy(p);
  nat mul1=0;
  nat mul2=0;
  C s;

  for (nat i=0; i<N(p); i++){
    
    // Factor x-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;
    }

    // Factor x+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) f << Sparse_factor (x-1,mul1);
  if (mul2 > 0) f << Sparse_factor (x+1,mul2);

  return f;
}

template<typename E, typename C, typename V>
pair<C,vector <Sparse_factor> >
gap_and_cyclotomic_factors (const Sparse_polynomial &p ){

  Sparse_polynomial q= copy(p);
  
  // For the factors
  vector <Sparse_factor> f;
  Sparse_polynomial x ((C)1, Monomial (vec(1)));
  nat mul1=0;
  nat mul2=0;
  C s1, s2;

  // For the gap
  C g=0;

  for (nat i=0; i<N(p); i++){
    
    if (mul1 == i) s1=get_coefficient(q,0); 
    if (mul2 == i)
      s2=get_coefficient(q,0) * binpow(-1, deg (get_monomial(q,0)) % 2);
    C max_coeff= abs (get_coefficient(q,0));

    for (nat j=1; j<N(q); j++) {
   
      if (mul1 == i) s1+= get_coefficient (q,j); 
      if (mul2 == i) 
        s2+= get_coefficient (q,j) * binpow(-1, deg (get_monomial(q,j)) % 2);

      C curr_coeff = abs (get_coefficient (q,j));
      if (curr_coeff > max_coeff) max_coeff= curr_coeff;
    }

    if (mul1 == i && s1 == 0) mul1=i+1;
    if (mul2 == i && s2 == 0) mul2=i+1;
    
    C length= as<C> (N(q)-1);
    g= max (g, as<C> (bit_size(length*max_coeff)) );
    q= sparse_derive (q);
  }

  if (mul1 > 0) f << Sparse_factor (x-1,mul1);
  if (mul2 > 0) f << Sparse_factor (x+1,mul2);

  pair<C, vector <Sparse_factor> > result(g,f);

  return result;
}

template<typename E, typename C, typename V>
Sparse_factors
linear_factors_univariate ( const Sparse_polynomial &q ){
  Sparse_polynomial x ((C)1, Monomial (vec(1)));
  Sparse_factors f;
  nat v= as<nat> (val(q));
  
  pair<C, vector <Sparse_factor> > gap_cyclo= gap_and_cyclotomic_factors (q);
  E gap= as<E> (car (gap_cyclo));

  Sparse_factors cyclotomic_fact= cdr (gap_cyclo);

  if ( v > 0 ) f << Sparse_factor( x , v );
  f << cyclotomic_fact; 
  f << non_cyclotomic_factors (q,gap);
  return f;
}

#undef Monomial
#undef Sparse_polynomial
#undef Univariate_polynomial
#undef Univariate_factor
#undef Univariate_factors
#undef Sparse_factor
#undef Sparse_factors

} // mmx namespace

---