/Users/mourrain/Devel/mmx/realroot/include/realroot/subresultant.hpp

Go to the documentation of this file.

#ifndef subdivix_subresultant_hpp
#define subdivix_subresultant_hpp
/********************************************************************/
#include <realroot/Seq.hpp>
#include <realroot/polynomial_univariate.hpp>
namespace mmx {
struct euclidean {


  template<class R> static
  void pseudo_div(R& q, R& a, const R& b)
  {
    typedef typename R::iterator iterator;
    R t;//(degree(b)+1, AsSize());
    typename R::value_type lb=b[degree(b)];
    q=R(0);
    int c=0, da = degree(a), db = degree(b), dab =da-db+1;
    while (da>=db) {
      R m(a[da],da-db);
      shift(t,b,da-db); t*=a[da];   // t=b*m;
      c++;
      a*=lb;
      t[da]=a[da];
      a-=t;
      q+=m;
      da=degree(a);
    }
    a *= pow(lb, dab-c);
  }

  template<class Polynomial> static inline 
  Polynomial prem(const Polynomial& a, const Polynomial& b)
  {
    Polynomial r=a,q;
    pseudo_div(q, r, b);
    return r;
  }

  template<class Polynomial> static inline 
  Polynomial pquo(const Polynomial& a, const Polynomial& b)
  {
    Polynomial r=a, q;
    pseudo_div(q, r, b);
    return q;
  }

};



template<class PREM>
struct subresultant {

  template<class Pol> 
  static Seq<Pol> sequence(const Pol & p, const Pol & q) {

    typedef typename Pol::value_type coeff_type;
    assert(degree(p) >=0 && degree(q) >= 0);
    int n1 = degree(p), n2= degree(q);
    coeff_type lp = p[n1], lq = q[n2];
    
    Seq<Pol> tab;
    unsigned N = std::min(n1,n2);
    tab.resize(N);
    for(unsigned i=0;i<N;i++) tab[i]=Pol(0);
    int d,da,db;
    coeff_type la;
    Pol a, b, Sr;
    
    if (n1>=n2) {
        a=q;
        b=PREM::prem(p,q);
        d=n1-n2-1;
        la=pow(lq,n1-n2);
    } else {
      a=p;
      b=PREM::prem(q,p);
      d=n2-n1-1;
      la=pow(lp,n2-n1);
    }
    
    if (d!=0) {
      a*=(coeff_type)pow<int>(-1,d*(d+1)/2);
      b*=(coeff_type)pow<int>(-1,(d+2)*(d+3)/2);
    }
    da=degree(a);
    db=degree(b);
    if(db>=0) tab[db]=b;
    
    while (db>=0) {
      Sr =PREM::prem(a,b);
      Sr/=(coeff_type)(pow(la,da-db)*a[da]);
      Sr*=(coeff_type)pow<int>(-1,((da-db+1)*(da-db+2))/2);
      a  = b;
      a *=(coeff_type)pow(b[db],da-db-1);
      a /=(coeff_type)pow(la,da-db-1); 
      da=degree(a); 
      b=Sr; 
      la=a[da]; 
      db=degree(b);
      if (db>=0) tab[db]=b;
    }
    return tab;
  }

  template<class Polynomial> 
  static Seq<Seq<Polynomial> > sequence(const Polynomial & p, const Polynomial & q, int v) {

    Seq<Polynomial> 
      sp= coefficients(p,v),
      sq= coefficients(q,v);

    typedef polynomial< Polynomial, with<Univariate> > upolmpol_t;
    upolmpol_t 
      P(1,sp.size()-1),
      Q(1,sq.size()-1);

    for (unsigned i=0;i<sp.size();i++) P[i]= sp[i];
    for (unsigned i=0;i<sq.size();i++) Q[i]= sq[i];

    Seq< upolmpol_t > s = subresultant<euclidean>::sequence(P,Q);
    
    Seq< Seq<Polynomial> > res(s.size());
    for(unsigned i=0;i<s.size();i++)
      for(int j=0; j<=degree(s[i]);j++)
        res[i]<<s[i][j];

    return res;
  }

  template<class Pol> static typename Pol::value_type 
  resultant(const Pol & p, const Pol & q) {
    Seq<Pol> s= sequence(p,q);
    if(s.size()>0) 
      return s[0][0];
    else return 0;
  }
    
  template<class Pol> static Pol 
  resultant(const Pol & p, const Pol & q, int v) {
    Seq< Seq<Pol> > s = sequence(p,q,v);
    if(s.size()>0) return s[0][0];
    else return 0;
  }

  template<class Pol> static bool 
  is_zero_seq(const Seq<Pol> & s) {
    for (unsigned i=0;i<s.size() ;i++) 
      if(s[i]!=0) return false;
    return true;
  }
    
  template<class Polynomial> static Polynomial 
  subres_gcd(const Polynomial & p, const Polynomial & q, int v) {
    typedef typename Polynomial::Monomial Monomial;    
    if(p==0) return q;
    if(q==0) return p;
    if(v==-1) return gcd(content(p),content(q));

    if(nbvar(p)==v+1 && nbvar(q)<= v) {
      Seq< Polynomial > s = coefficients(p,v);
      Polynomial d=q;
      for (unsigned i=0;i<s.size();i++) 
        d = subres_gcd(d,s[i]);
      return d;
    }

    if(nbvar(p)<=v && nbvar(q)== v+1) {
      Seq< Polynomial > s = coefficients(q,v);
      Polynomial d=p;
      for (unsigned i=0;i<s.size();i++) 
        d = subres_gcd(d,s[i]);
      return d;
    }


    // Remove content of p
    Polynomial dp=0;
    Seq<Polynomial> sp = coefficients(p,v);
    for(unsigned i=0;i< sp.size();i++) dp= subres_gcd(dp,sp[i]);
    Polynomial P= p/dp;

    // Remove content of q
    Polynomial dq=0;
    Seq<Polynomial> sq = coefficients(q,v);
    for(unsigned i=0;i< sq.size();i++) dq= subres_gcd(dq,sq[i]);
    Polynomial Q= q/dq;

    // Take gcd of contents
    Polynomial C= subres_gcd(dp,dq);

    // Subresultant sequence
    Seq< Seq<Polynomial> > s 
      = (degree(P)>=degree(Q)?sequence(P,Q,v):sequence(Q,P,v));

    unsigned i;
    for (i=0;i<s.size() && is_zero_seq(s[i]);i++) ;

    if(i==s.size()) return (degree(P)>=degree(Q)?Q*C:P*C);

    // Take last non-zero term
    Seq<Polynomial> D = s[i];

    // Remove content of last non-zero term.
    Polynomial d=0;
    for (i=0;i<D.size();i++) d = subres_gcd(d,D[i]);
    //    d/=content(d);
    for (i=0;i<D.size();i++) D[i]/=d;

    // Reconstruct the multivariate polynomial
    Polynomial r=0;
    for (i=0;i<D.size();i++) 
      r+= D[i]*Polynomial(Monomial(1,i,v));
    r/=content(r);
    r *= C;
    return r;
  }

  template<class Polynomial> static Polynomial 
  subres_gcd(const Polynomial & p, const Polynomial & q) {
    int v = std::max(nbvar(p)-1,nbvar(q)-1);
    return subres_gcd(p,q,v);
  }
};

  template<class C, class V> inline polynomial<C,V>
  gcd (const polynomial<C,V> & p, const polynomial<C,V> & q) {
    return subresultant<euclidean>::subres_gcd(p,q);
  }
}// namespace mmx
#endif //subdivix_subresultant_hpp