include/lacunaryx/linear_factors_bivariate.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : linear_factors_bivariate.hpp
* DESCRIPTION: Linear factors of bivariate 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 <lacunaryx/linear_factors_univariate.hpp>
#include <factorix/irreducible_polynomial_polynomial.hpp>

namespace mmx {

#define Monomial              monomial< vector<E> >
#define Sparse_polynomial     sparse_polynomial<C, Monomial, V>
#define Univariate_polynomial polynomial<C>
#define Bivariate_polynomial  polynomial<polynomial<C> >
#define Bivariate_factor      irreducible_factor< polynomial <polynomial<rational> > >
#define Bivariate_factors     vector<Bivariate_factor>
#define Sparse_factor         irreducible_factor<Sparse_polynomial>
#define Sparse_factors        vector<Sparse_factor>

namespace lacunaryx {
  template<typename C, typename V1, typename V2, typename M, typename W>
  void set_as_bivariate_polynomial (polynomial<polynomial<C,V2>,V1>& p,
                                    const sparse_polynomial<C,M,W>& q) {
  nat n1= 1 + as<nat> (deg (q,0));
  vector<nat> n2 (0u, n1);
  for (nat i= 0; i < N(q); i++) {
    nat e1= as<nat> (deg (get_monomial (q, i),0));
    nat e2= as<nat> (deg (get_monomial (q, i),1));
    n2[e1]= max (n2[e1], e2);
  }
  n2= n2+1;
  vector<C*> buf ((C*) NULL,n1);
  vector<nat> l2 (0u, n1);
  for (nat i= 0; i < n1; i++) {
    l2[i]= aligned_size<C,V2> (n2[i]);
    buf[i]= mmx_new<C> (l2[i], promote (0,CF(q)));
  }

  for (nat i= 0; i < N(q); i++) {
    nat e1= as<nat> (deg (get_monomial (q, i),0));
    nat e2= as<nat> (deg (get_monomial (q, i),1));
    buf[e1][e2]= get_coefficient (q, i);
  }
  nat l1= aligned_size<C,V2> (n1);
  polynomial<C,V2>* c1= mmx_new<polynomial<C,V2> > (l1);
  for (nat i= 0; i < n1; i++)
    c1[i]= polynomial<C,V2> (buf[i],n2[i],l2[i],CF(q));
  polynomial<C,V2> sample (CF(q));  
  p= polynomial<polynomial<C,V2>,V1> (c1, n1, l1, get_format (sample));
    
  }

  template<typename C, typename M, typename W>
  polynomial<polynomial<C> > as_bivariate_polynomial (const sparse_polynomial<C,M,W> &q) {
    polynomial<C> sample (CF(q));  
    polynomial<polynomial<C> > p (get_format (sample));
    set_as_bivariate_polynomial(p,q);
    return p;
  }

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

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

  template<typename C, typename E, typename V> 
  Sparse_polynomial homogenize(const Sparse_polynomial& p) {
    //ASSERT (dim(p) == 1, "wrong dimension"); 
    E d= deg(p);
    vector<Monomial> v= supp(p);
    for (nat i=0; i<N(v); i++) {
      E dm= deg (v[i],0);
      v[i]= Monomial (vec (dm,d-dm));
    }
    return Sparse_polynomial ( coefficients(p), v, 2);
  }
  
  template<typename E, typename C, typename V>
  E w_deg (const Sparse_polynomial& p, vector<E> w) {
    return dot (multi_deg(p), w);
  }
  
  template<typename E>
  E w_deg (const Monomial m, vector<E> w) {
    return dot (multi_deg(m),w);
  }
  
  // Sort the monomials according the their weighted degree
  template<typename E, typename C, typename V>
  vector<nat> sort_poly (const Sparse_polynomial& q, vector<E> w) {
    vector<E> degrees= fill<E>(0,N(q));
    vector<nat> order= fill<nat>(0,N(q));
    
    for (nat i=0; i<N(q); i++) {
      Monomial m= get_monomial(q,i);
      degrees[i]= w_deg(m,w);
    }
  
    sort (degrees, order);
    return order;
  }
  
  // Sort the monomials according to their degree in one specific variable
  template<typename E, typename C, typename V>
  vector<nat> sort_poly (const Sparse_polynomial& q, nat var) {
    vector<E> degrees= fill<E>(0,N(q));
    vector<nat> order= fill<nat>(0,N(q));
    
    for (nat i=0; i<N(q); i++) {
      Monomial m= get_monomial(q,i);
      degrees[i]= deg(m,var);
    }
  
    sort (degrees, order);
    return order;
  }
  
  template<typename E, typename C, typename V>
  vector<Sparse_polynomial> 
  w_homogeneous_parts (const Sparse_polynomial& p, vector<E> w) {
    vector<nat> o= sort_poly(p,w);
    vector<Sparse_polynomial> v;
    
    E d= w_deg (get_monomial (p,o[0]), w);
    nat imin=0;
    for (nat i=1; i<N(p); i++) {
      if (w_deg (get_monomial (p,o[i]), w) > d) {
        v << range_permute(p, imin, i, o);
        imin=i;
        d= w_deg (get_monomial (p,o[i]), w);
      }
    }
    v << range_permute(p, imin, N(p), o);
    return v;
  }

} // lacunaryx namespace

using namespace lacunaryx;

// Partition with respect to one specific variable
template<typename E, typename C, typename V>
vector<Sparse_polynomial> 
partition_bivariate_helper (const Sparse_polynomial &q, nat var ) {

  vector<nat> order= sort_poly( q, var );
  
  E min_exponent= deg(get_monomial(q,order[0]),var);
  E curr_exponent;
  E accu=0;
  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,order[i]),var);
  
    // The bound is curr_exponent - min_exponent > binom{number_of_terms}{2}
    // It is iteratively computed in accu
    if ( (curr_exponent - min_exponent) > accu ) {
      tmp= range_permute(q,index_min,i,order);
      tmp/= gcd_supp (tmp);
      polys << tmp;
      accu= 0;
      min_exponent= curr_exponent;
      index_min= i;
    } else {
      accu+=i;
    }
  }
  tmp= range_permute(q, index_min, N(q), order);
  tmp/= gcd_supp (tmp);
  polys << tmp;
  return polys;
}

// Partition with respect to both variables
template<typename E, typename C, typename V>
vector<Sparse_polynomial> 
partition_bivariate(const Sparse_polynomial &p ) {

  vector<Sparse_polynomial> current= vec(p);
  nat prev_length=0;
  nat curr_length= N( current );

  while ( prev_length < curr_length ) {
    
    prev_length= curr_length;

    // Partition according to both variables
    for (nat i=0; i<2; i++) {
      vector<Sparse_polynomial> next;
      for (nat j=0; j<curr_length; j++) {
        next << partition_bivariate_helper( current[j], i );
      }
      current= next;
      curr_length= N(current);
    }
  }

  return current;
}

template<typename E, typename C, typename V>
Sparse_factors
truly_bivariate_linear_factors ( const Sparse_polynomial &p ) {
  vector<Sparse_polynomial> v= partition_bivariate (p);
  
  Bivariate_polynomial q= as_bivariate_polynomial (v[0]);
  for (nat i=1; i<N(v); i++) {
    q= gcd (q, as_bivariate_polynomial (v[i]) );
  }

  polynomial <polynomial<rational> > q_rat= 
    as< polynomial <polynomial<rational> > > (q);
  Bivariate_factors f= // irreducible_factorization(q_rat);
        bounded_irreducible_factorization(q_rat,2);
  Sparse_factors ret;
  for (nat i=0; i<N(f); i++) {
    if (deg(f[i].f)==1 && deg(f[i].f[0])==1 && val(f[i].f[0])==0)
      ret << as_sparse_factor<C,E,V> (f[i]);
  }
  return ret;
}

template<typename E, typename C, typename V>
Sparse_factors
linear_factors_bivariate_naive (const Sparse_polynomial &p) {
  Sparse_factors factors;
  vector<Sparse_polynomial> hom_comp;
  Sparse_factors lin_fact;
  
  // Univariate factors, first variable
  hom_comp= w_homogeneous_parts(p, vec((E)0, (E)1));
  lin_fact= linear_factors_univariate (project (hom_comp[0],vec(0)));
  for (nat i=1; i<N(hom_comp); i++) {
    Sparse_factors tmp_fact= 
      linear_factors_univariate (project (hom_comp[i],vec(0))); 
    lin_fact= intersect_factors (tmp_fact, lin_fact);
  }
  for (nat i=0; i<N(lin_fact); i++)
    factors << Sparse_factor (inject (lin_fact[i].f, vec(0), 2), lin_fact[i].m);

  // Univariate factors, second variable
  hom_comp= w_homogeneous_parts(p,vec((E)1, (E)0));
  lin_fact= linear_factors_univariate (project (hom_comp[0],vec(1)));
  for (nat i=1; i<N(hom_comp); i++) {
    Sparse_factors tmp_fact= 
      linear_factors_univariate (project (hom_comp[i],vec(1)));
    lin_fact= intersect_factors (tmp_fact, lin_fact);
  }
  for (nat i=0; i<N(lin_fact); i++) 
    factors << Sparse_factor (inject (lin_fact[i].f, vec(1), 2), lin_fact[i].m);

  // Bivariate factors, binomials
  hom_comp= w_homogeneous_parts(p,vec((E)1, (E)1));
  lin_fact= linear_factors_univariate (project (hom_comp[0],vec(0)));
  for (nat i=1; i<N(hom_comp); i++) {
    Sparse_factors tmp_fact= 
      linear_factors_univariate (project (hom_comp[i],vec(0)));
    lin_fact= intersect_factors (tmp_fact, lin_fact);
  }
  for (nat i=0; i<N(lin_fact); i++)
    if (val (lin_fact[i].f) == 0) // remove monomials 
      factors << Sparse_factor (homogenize (lin_fact[i].f), lin_fact[i].m);
  
  // Truly bivariate factors
  Sparse_factors tbf= truly_bivariate_linear_factors(p);
  for (nat i=0; i<N(tbf); i++) factors << tbf[i];

  return factors;
}

template<typename E, typename C, typename V>
Sparse_factors
linear_factors_bivariate (const Sparse_polynomial &p) {
  Sparse_factors factors;
  vector<Sparse_polynomial> hom_comp;
  Sparse_factors lin_fact;
  
  // Univariate factors, first variable
  hom_comp= w_homogeneous_parts(p, vec((E)0, (E)1));
  hom_comp[0]= project (hom_comp[0],vec(0));
  nat Nmin= N(hom_comp[0]); 
  nat index=0;
  for (nat i=1; i<N(hom_comp); i++) {
    hom_comp[i]= project (hom_comp[i],vec(0));
    if (N(hom_comp[i]) < Nmin) { 
      Nmin= N(hom_comp[i]); 
      index=i;
    }
  }
  lin_fact= linear_factors_univariate (hom_comp[index]);
  for (nat i=0; i<N(hom_comp); i++) {
    if (i != index) {
      Sparse_factors tmp_fact;
      for (nat j=0; j<N(lin_fact); j++) {
        nat m= multiplicity (lin_fact[j].f, hom_comp[i]);
        if (m > 0) 
          tmp_fact << Sparse_factor (lin_fact[j].f, min(lin_fact[j].m,m));
      }
      lin_fact= tmp_fact;
    }
  }
  for (nat i=0; i<N(lin_fact); i++)
    factors << Sparse_factor (inject (lin_fact[i].f, vec(0), 2), lin_fact[i].m);

  // Univariate factors, second variable
  hom_comp= w_homogeneous_parts(p,vec((E)1, (E)0));
  hom_comp[0]= project (hom_comp[0],vec(1));
  Nmin= N(hom_comp[0]); 
  index=0;
  for (nat i=1; i<N(hom_comp); i++) {
    hom_comp[i]= project (hom_comp[i],vec(1));
    if (N(hom_comp[i]) < Nmin) { 
      Nmin= N(hom_comp[i]); 
      index=i;
    }
  }
  lin_fact= linear_factors_univariate (hom_comp[index]);
  for (nat i=0; i<N(hom_comp); i++) {
    if (i != index) {
      Sparse_factors tmp_fact;
      for (nat j=0; j<N(lin_fact); j++) {
        nat m= multiplicity (lin_fact[j].f,hom_comp[i]);
        if (m > 0) 
          tmp_fact << Sparse_factor (lin_fact[j].f, min(lin_fact[j].m,m));
      }
      lin_fact= tmp_fact;
    }
  }
  for (nat i=0; i<N(lin_fact); i++) 
    factors << Sparse_factor (inject (lin_fact[i].f, vec(1), 2), lin_fact[i].m);

  // Bivariate factors, binomials
  hom_comp= w_homogeneous_parts(p,vec((E)1, (E)1));
  hom_comp[0]= project (hom_comp[0],vec(0));
  Nmin= N(hom_comp[0]); 
  index=0;
  for (nat i=1; i<N(hom_comp); i++) {
    hom_comp[i]= project (hom_comp[i],vec(0));
    if ( N(hom_comp[i]) < Nmin) { 
      Nmin= N(hom_comp[i]); 
      index=i;
    }
  }
  lin_fact= linear_factors_univariate (hom_comp[index]);
  for (nat i=0; i<N(hom_comp); i++) {
    if (i != index) {
      Sparse_factors tmp_fact;
      for (nat j=0; j<N(lin_fact); j++) {
        nat m= multiplicity (lin_fact[j].f, hom_comp[i]);
        if (m > 0) 
          tmp_fact << Sparse_factor (lin_fact[j].f, min(lin_fact[j].m,m));
      }
      lin_fact= tmp_fact;
    }
  }
  for (nat i=0; i<N(lin_fact); i++)
    if (val (lin_fact[i].f) == 0) // remove monomials 
      factors << Sparse_factor (homogenize (lin_fact[i].f), lin_fact[i].m);
  
  // Truly bivariate factors
  Sparse_factors tbf= truly_bivariate_linear_factors(p);
  for (nat i=0; i<N(tbf); i++) factors << tbf[i];

  return factors;
}

#undef Monomial
#undef Sparse_polynomial
#undef Univariate_polynomial
#undef Bivariate_polynomial
#undef Bivariate_factor
#undef Bivariate_factors
#undef Sparse_factor
#undef Sparse_factors

} // mmx namespace

---