synaps/upol/SturmSeq.h

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s):
 *       Initial version:
 *          B. Mourrain, INRIA, GALAAD, <mourrain@sophia.inria.fr>
 *       Implementatation of Sturm Habicht by
 *         F. Carreras, U. of Santander, <fernando.carreras@unican.es>
 *         C. Tansescu, U. of Santander, <claudiu.tanasescu@unican.es>
 *         L. Gonzales-Vega, U. of Santander, <laureano.gonzalez@unican.es>
 *         O. Ruatta, INRIA, GALAAD, <oruatta@sophia.inria.fr>
 *       Extended list of SturmSeq constructors by
 *         Elias P. TSIGARIDAS, U. of Athenes,  <et@di.uoa.gr>
 *
 *   $Id: SturmSeq.h,v 1.1 2005/07/11 08:15:52 mourrain Exp $
 ********************************************************************/
#ifndef SYNAPS_RESULTANT_STURM_H
#define SYNAPS_RESULTANT_STURM_H
/********************************************************************/
#include <synaps/init.h>
#include <synaps/arithm/sign.h>
#include <synaps/arithm/pow.h>
#include <synaps/arithm/Infinity.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/UPOLDAR.m>
#include <synaps/upol/EUCLIDIAN.h>
#include <synaps/upol/SturmNames.h>
#include <synaps/arithm/gmp.h>
#include <vector>
/********************************************************************/
__BEGIN_NAMESPACE_SYNAPS
//====================================================================
//--------------------------------------------------------------------
template<class UPoly> inline
bool isNull(const UPoly& b) {return b==0; }


template<class UPoly> inline
int parity(const UPoly& p)
{
  return (p.get_degree() % 2 ? -1 : 1);
}

template <class Pol> inline
Pol Prem(const Pol & a, const Pol & b)
{
  Pol r(a),q;
  EUCLIDIAN::div_pseudorem0(q,r,b);
  return r;
}


namespace sturmseq {};


//--------------------------------------------------------------------
using UPOLDAR::degree;
using UPOLDAR::lcoeff;
//--------------------------------------------------------------------




// /** Construct the Sturm sequence of two polynomials @a@ and @b@.
//  * The Sturm sequence, is actually a signed pseudo remainder sequence
//  * and can be computed by various methods.
//  */
// template<class C, class UPoly, class SEQ>
// SturmSeq<C, UPoly, SEQ>::SturmSeq(const UPoly & a, const UPoly& b, PRSmethod method)
//   : mth(method), tab_(std::min(degree(a), degree(b)) + 3)
// {
//   switch (method) {
//   case PRS_STANDARD:
//     compute_Standard(a, b);
//     break;
//   case PRS_EUCLIDEAN:
//     compute_Euclidean(a, b);
//     break;
//   case PRS_REDUCED:
//     compute_Reduced(a, b);
//     break;
//   case PRS_SUBRESULTANT:
//     compute_Subres(a, b);
//     break;
//   case PRS_STHA:
//     compute_StHa(a,b);
//     break;
//   default:
//     compute_Standard(a, b);
//     break;
//   }
// }

//====================================================================
template<class C, class UPoly, class SEQ> struct  SturmSeq;
//--------------------------------------------------------------------
namespace sturmseq
{
template<class C, class UPoly, class SEQ>
void init(STD mthd, SturmSeq<C,UPoly,SEQ> & tab,
                   const UPoly & a, const UPoly& b)
{
  //  cout << __FUNCTION__ << endl;
  tab[0] = a;
  tab[1] = b;
  unsigned i = 1;
  while (i < tab.size()-1 && !isNull(tab[i]) )
    {
      tab[i+1] = - (tab[i-1] % tab[i]);
      i++;
    }
  if ( isNull(tab[i]) ) tab.resize(i);
}



template< class C, class UPoly, class SEQ >
void init(EUCLIDEAN mth, SturmSeq<C,UPoly,SEQ> & tab,
                   const UPoly & a, const UPoly& b)
{
  //  cout << __FUNCTION__ << endl;
  tab[0] = a;
  tab[1] = b;

  unsigned i;
  for (i = 1; degree(tab[i]) > 0; ++i) {
        tab[i+1] = EUCLIDIAN::prem(tab[i-1], tab[i])
          *(-1)*pow( sign(EUCLIDIAN::lcoeff(tab[i])), degree(tab[i-1]) - degree(tab[i]) + 1);
  }
  tab.resize(i+1);
  return;
}


template< class C, class UPoly, class SEQ >
void init(REDUCED mth, SturmSeq<C,UPoly,SEQ> & tab,
          const UPoly & a, const UPoly& b)
{
  //  cout << __FUNCTION__ << endl;
  unsigned i = 1;
  tab[0] = a;
  tab[1] = b;
  C alpha = tab[1][tab[1].get_degree()];
  int delta = tab[0].get_degree() - tab[1].get_degree();
//  C psi = C(1);
  C beta = C(1);
  while ( (i < tab.size()-1) && (tab[i].get_degree() > 0) )
    {
      tab[i+1] = EUCLIDIAN::prem(tab[i-1], tab[i]);

      C aa = pow( sign(lcoeff(tab[i])), tab[i-1].get_degree()-tab[i].get_degree()+1);
        int sgn_aa = sign(aa);
        int sgn = (sgn_aa  == 1
                ? (sign(beta) == 1 ? -1 :  1 )
                : (sign(beta) == 1 ?  1 : -1 )
                );
        tab[i+1] /= beta*sgn;

        delta = tab[i-1].get_degree() - tab[i].get_degree();
        alpha = tab[i-1+1][tab[i-1+1].get_degree()];
        beta = pow(alpha, delta + 1);
        ++i;
    }
  if (i+1 < tab.size()) tab.resize(i+1);
}


template< class C, class UPoly, class SEQ >
void init(SUBRES, SturmSeq<C, UPoly, SEQ> & tab,
          const UPoly & a, const UPoly& b)
{
  int du = a.get_degree();
  int dv = b.get_degree();
  C g(1);
  C h(1);
  int i = 1;

  tab[0] = a;
  tab[1] = b;
  if (du < dv) {
        std::swap(du, dv);
        tab[++i] = -a;
  }

//  int sgn = 1;
  C aa = 1;
  UPoly r;
  while (dv > 0) {
        int delta = du - dv;
        aa = pow( sign(lcoeff(tab[i])), delta+1 );

        tab[i+1] = EUCLIDIAN::prem(tab[i-1], tab[i]);
        du = dv;
        dv = tab[i+1].get_degree();

        C beta = g * pow(h, delta);
        int sgn_aa = sign(aa);
        int sgn = (sgn_aa  == 1
                ? (sign(beta) == 1 ? -1 :  1 )
                : (sign(beta) == 1 ?  1 : -1 )
                );

        if (beta != 0) tab[i+1] /= beta*sgn;
        else tab[i+1] = tab[i+1] * sgn;
        
        g = lcoeff(tab[i]);
        if (delta == 1) {
          h = g;
        } else {
//        h = pow(g, delta) / pow(h, delta - 1);
          h = pow( C(g / h), delta)*h;
        }
        ++i;
  }

  if (i+1 < tab.size()) tab.resize(i+1);
  return;
}



template<class C, class Pol, class TAB_TYPE>
void init(HABICHT mth,SturmSeq<C,Pol,TAB_TYPE> & tab,
          const Pol & a, const Pol & b)
{
        typedef typename Pol::value_type coeff_type;
        assert(degree(a) >=0 && degree(b) >= 0);
        Pol p=a,
                        q=b;
        int n1 = degree(p), n2= degree(q);
        coeff_type lp = p[n1], lq = q[n2];

        tab.resize(n1+n2+3);

        tab[0]=p;
        tab[1]=q;
        int j,n,d,m,nS1,nS0,r;
        //      int i, k;
        Pol S0, S1, NS, Sr;
        int c=2;
        //m==0;
        if (n1>n2)
        {
          n=n1-1; d=n-n2;
          if (d!=0)
            {
              S1=q*pow(-1,d*(d+1)/2)*pow(coeff_type(lp*lq),d);
              S0=Prem(p,q)*pow<int>(-1,(d+2)*(d+3)/2)*pow(lp,d);
              if (S0!=Pol(0))
                {tab[c++]=S1; tab[c++]=S0;}
              else
                {tab[c++]=S0;}
            }
          else
            {
              S1=q;S0=S0-Prem(p,q);
              if (S0!=Pol(0)) {tab[c++]=S0;}
            }
          j=n2-1;
        }
        else
        {
                n=n2; S1=q; S0=S0-((p*(lq*lq))%q); j=n-1;
                if (S0!=Pol(0)) {tab[c++]=S0;}
        }
        nS1=degree(S1);
        coeff_type R=S1[nS1];
        if (S0==Pol(0)) r=-1; else r=degree(S0);

        while (r>=0)
        {
                if (r>0)    {
                        if (r==j)
                        {
                          coeff_type denom=pow<coeff_type>(-R,j-r+2);
                          NS=Prem(S1,S0); NS/= denom;
                          NS*=pow<coeff_type>(-1,((j-r+2)*(j-r+3))/2);
                          if (NS!=Pol(0)) {tab[c++] =NS;}
                          S1=S0; S0=NS; nS1=degree(S1);R=S1[nS1]; j=r-1;
                        }
                        else
                        {
                                for (m=r+1;m<j;m++)
                                {
                                        tab[c++]= Pol(0);
                                }
                                Sr=Prem(S1,S0); Sr/=pow(R,j-r+2);
                                nS0=degree(S0);
                                S1=S0*(pow(S0[nS0],j-r)); S1/=pow(R,j-r);
                                S1*=coeff_type(pow<int>(-1,((j-r)*(j-r+1)/2)));
                                S0=Sr*(pow<int>(-1,((j-r+2)*(j-r+3)/2)));
                                if (S0!=Pol(0)) {
                                        tab[c++]=(S1);
                                        tab[c++]=(S0);}
                                else {
                                        tab[c++]=(S1);}
                                nS1=degree(S1);
                                R=S1[nS1]; j= r-1;
                        }
                }
                else
                {
                        if (j>0)
                        {
                                nS0=degree(S0); S1=S0*pow(S0[nS0],j); S1 /=pow(R,j);
                                S1*=coeff_type(pow<int>(-1,(j*(j+1)/2))); tab[c++]=(S1);
                        }
                        S0=Pol(0);
                }
                if (S0!=Pol(0)) r=degree(S0); else r=-1;
        }
        tab.resize(c);
}
//--------------------------------------------------------------------
}//sturmseq
//====================================================================
template< class C,
  class UPoly = UPolDse<C>,
  class SEQ = std::vector<UPoly> >
struct SturmSeq
{
  typedef UPoly value_type;
  //SEQ         tab_;
  shared_object<SEQ>  data;
  SEQ &       rep()       {return (*data);}
  const SEQ & rep() const {return (*data);}


  //definition of the types for iterator on Seq
  typedef typename SEQ::iterator       iterator;
  typedef typename SEQ::const_iterator const_iterator;
  typedef std::reverse_iterator<const_iterator>  const_reverse_iterator;
  typedef std::reverse_iterator<iterator>        reverse_iterator;

  //  void compute_Standard(const UPoly & a, const UPoly& b);
  //  void compute_Euclidean(const UPoly & a, const UPoly& b);
  //  void compute_Reduced(const UPoly & a, const UPoly& b);
  //  void compute_Subres(const UPoly & a, const UPoly& b);
  //  void compute_StHa(const UPoly & a, const UPoly& b);

  SturmSeq(): data() {}
  SturmSeq(const SEQ & r):data(r){}
  SturmSeq(const shared_object<SEQ> & r):data(r){}

  template<class MTH>
  SturmSeq(const UPoly & a, const UPoly& b, MTH  method =MTH())
      {
        rep().resize(std::min(degree(a), degree(b)) + 3);
        sturmseq::init(method,*this,a,b);
      }

  //  template<class R>  unsigned var(const R& x);
  //  unsigned var(const char * x);

  const UPoly& operator[](const int& i) const { return rep()[i]; }
        UPoly& operator[](const int& i) { return rep()[i]; }

  int  size() const        { return rep().size(); }
  void resize (unsigned i) {rep().resize(i);}
  void push_back(const UPoly & a) { return rep().push_back(a); }

};
//--------------------------------------------------------------------
template<class C, class UPoly, class SEQ, class R>
unsigned variation(const SturmSeq<C, UPoly, SEQ> & tab, const R& x)
{
  int i = 0, c = 0;
  int s0 = 0, s = 0;
  while (  i<tab.size() && (s0 = UPOLDAR::sign_at(tab[i],x)) == 0 )
        i++;
  while ( i < tab.size() ) {
    while( i<tab.size() && (s=UPOLDAR::sign_at(tab[i],x))==0)
          ++i;
        if(s !=0 && s!=s0) {
          c++;
          s0 = s;
        }
        ++i;
  }
  return c;
}



template<class C,class UPoly, class SEQ>
unsigned variation(const SturmSeq<C, UPoly, SEQ> & tab, 
                   const Infinity& x)
{
  unsigned i = 0, c = 0;
  int s0 = 0, s = 0;

  if (x.sign < 0) {
        while( (i < tab.size()) && ( (s0 = parity(tab[i])*sign(lcoeff(tab[i]))) == 0) )
          i++;
        while (i < tab.size()) {
          while ( (i < tab.size()) && ((s = parity(tab[i])*sign(lcoeff(tab[i]))) == 0) )
                ++i;
          if ((s != 0) && (s != s0)) {
                c++;
                s0 = s;
          }
          ++i;
        }
  } else {
        while( ( i < tab.size()) && ((s0 = sign(lcoeff(tab[i]))) == 0) )
          i++;
        while (i < tab.size()) {
          while ( (i < tab.size()) && ((s = sign(lcoeff(tab[i]))) == 0) )
                ++i;
          if ((s != 0) && (s != s0)) {
                c++;
                s0 = s;
          }
          ++i;
        }
  }
  return c;
}


template<class RET_TYPE, class C, class Pol, class TAB_TYPE >
RET_TYPE CoeffStHa(const SturmSeq<C,Pol,TAB_TYPE> & StHa)
{
  typename SturmSeq<C,Pol,TAB_TYPE>::const_iterator it = StHa.rep().begin();
  RET_TYPE RESULT;
  it++;

  for (; it != StHa.rep().end(); ++it) {
          RESULT.push_back(UPOLDAR::lcoeff(*it));
  }
  return RESULT;
}


template<class C, class UPoly, class SEQ>
std::ostream& operator<<(std::ostream &os, const SturmSeq<C,UPoly,SEQ> & tab)
{
  VECTOR::print(os,tab.rep());
  return os;
}

//====================================================================
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif // SYNAPS_RESULTANT_STURM_H