synaps/upol/resultant.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
*********************************************************************/
#ifndef SYNAPS_RESULTANT_UNIRESLT_H
#define SYNAPS_RESULTANT_UNIRESLT_H
//====================================================================
#include <complex>
#include <synaps/init.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/Det.h>
#include <synaps/linalg/hankel.h>
__BEGIN_NAMESPACE_SYNAPS
//====================================================================
//====================================================================

template <class MATR, class R>
MATR bezout(const UPolDse<R> & P, const UPolDse<R> & Q) 
{
  unsigned int l = std::max(degree(P),degree(Q));
  MATR M(l,l);
  if(degree(P)<degree(Q)){
    for(int i=0; i<degree(P)+1; i++){
      for( int j=0; j<degree(Q)-i; j++)
        for( int l=0; l<=j; l++)
          M(i+j-l,i+l) += P[i]*Q[j+i+1];
      for( int j=0; j<i; j++)
        for( int l=0; l<=j; l++)
          M(i-1-l,i-1-j+l) -= P[i]*Q[i-1-j];
    }
  } else {
    for( int i=0; i<degree(Q)+1; i++){
      for( int j=0; j<degree(P)-i; j++)
        for( int l=0; l<=j; l++)
          M(i+j-l,i+l) -= Q[i]*P[j+i+1];
      for( int j=0; j<i; j++)
        for( int l=0; l<=j; l++)
          M(i-1-l,i-1-j+l) += Q[i]*P[i-1-j];
    }
  }
  return (M);
} 
//------------------------------------------------------------------------
template <class R>
MatrStr<linalg::hankel<typename R::value_type> > bezout(const R & P)
{
  typedef typename R::value_type C;
  typedef MatrStr<linalg::hankel<C> >  Mat;
  Mat H(degree(P),degree(P));
  for(int i=0; i<degree(P); i++) {
    (H.rep())[i]=P[i+1];
    if (i!=degree(P)-1) (H.rep())[degree(P)+i]=(C)0;
  }
  return H;
}
//-------------------------------------------------------------------------

template <class MATR, class R>
MATR sylvester(const UPolDse<R> & P, const UPolDse<R> & Q) 
{
  int l = degree(P)+degree(Q); 
  MATR M(l,l);
  for( int i=0; i<l; i++)
    for( int j=0; j<l; j++) M(i,j)=0;
  for( int i=0; i< degree(Q); i++)
    for( int j=0; j< degree(P)+1; j++)
      M(j+i,i)=P[j];
  for( int i=0; i< degree(P); i++)
    for( int j=0; j< degree(Q)+1; j++)
      M(j+i,degree(Q)+i)=Q[j];

  return (M);
}
//-------------------------------------------------------------------------
template <class C, class R, class MATR>
MATR sylvester(const UPolDse<C,R> & P, const UPolDse<C,R> & Q, char c)
{
  int dp = degree(P), dq = degree(Q);
  int l = dp+dq;
  MATR M(l,l);
  for( int i=0; i<l; i++)
    for( int j=0; j<l; j++) M(i,j)=0;
  for( int i=0; i< dq; i++)
    for( int j=0; j< dp+1; j++)
      M(j+i,i)=P[dp-j];
  for( int i=0; i< dp; i++)
    for( int j=0; j< dq+1; j++)
      M(j+i,dq+i)=Q[dq-j];

  return (M);
}

//-------------------------------------------------------------------------
#if 0
template<class UPOL, class MPOL>
MatrDse<UPOL> u_sylvester(const MPOL & p1, const MPOL & p2)
{
   int d1 = degree(p1,1), d2=degree(p2,1);
   MatrDse<UPOL> A(d1+d2,d1+d2);

   for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
   {
      UPOL m(it->coeff(),(*it)[0]);
      for(int k=0;k<d2;k++) A((*it)[1]+k,k) +=m;
   }
   for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
   {
      UPOL m(it->coeff(),(*it)[0]);
      for(int k=0;k<d1;k++) A((*it)[1]+k,k+d2) +=m;
   }

   return A;
}
#endif
//#if 0
template<class UPOL, class MPOL>
MatrDse<UPOL> u_sylvester(const MPOL & p1, const MPOL & p2,
        typename MPOL::monom_t::index_t i = 1)
{
  typename MPOL::monom_t::index_t j = (i == 0 ? 1 : 0);

  int d1 = degree(p1, i), d2 = degree(p2, i);
  MatrDse<UPOL> A(d1+d2, d1+d2);

  for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
  {
        UPOL m(it->coeff(),(*it)[j]);
//      cout << *it << " " << m << endl;
        //for(int k=0;k<d2;k++) A((*it)[0]+k,k) +=m;
        for(int k=0;k<d2;k++) A(k, k+ (*it)[i]) +=m;
  }
//  cout << "---" << endl;
  for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
  {
        UPOL m(it->coeff(),(*it)[j]);
//      cout << m << endl;
        //for(int k=0;k<d1;k++) A((*it)[0]+k,k+d2) +=m;
        for(int k=0;k<d1;k++) A(d2 + k, k+ (*it)[i]) +=m;
  }

  return A;
}
//#endif
//-------------------------------------------------------------------------
template <class MATR, class R>
MATR companion(const UPolDse<R> & P)
{
  typedef typename UPolDse<R>::size_type size_type;
  size_type n = degree(P);
  MATR M(n,n);
  for(size_type i=0; i<n; i++)
    for(size_type j=0; j<n-1; j++)
      if(i==j+1) M(i,j)=1; else M(i,j)=0;
  for(size_type i=0; i<n; i++) M(i,n-1)=-(P[i]/P[n]);

  return M;
}
//--------------------------------------------------------------------
template<class MPOL>
MatrDse<MPOL>  u_sylvester(const MPOL & p1, const MPOL & p2,int I)
{
  typedef typename MPOL::monom_t monom_t;
  int d1 = degree(p1,I), d2=degree(p2,I);
  MatrDse<MPOL> A(d1+d2,d1+d2);
  for(typename MPOL::const_iterator it=p1.begin();it!=p1.end();it++)
    {
      int d=(*it).nvars();
      MPOL P("1;");
      if(d<0)
        {
          P=(*it);
        }
      else
        {
          for(int j=0;j<=d;j++)
            {
              MPOL P1;
              if(j!=I)
                {
                  int power=(*it)[j];
                  if(power!=0)
                    {
                      P1= monom_t(1,j,power);
                      P*=P1; 
                    }
                }
              else if(j==I && (*it)[I]>0)
              {
                P1=monom_t(1,0,0);
                P*=P1;
              }
            }
          P*=(it->coeff());
        }
      for(int k=0;k<d2;k++)
        {
          A((*it)[I]+k,k) +=P;
        }
    }
  for(typename MPOL::const_iterator it=p2.begin();it!=p2.end();it++)
    {
      int d=(*it).nvars();
      MPOL P("1;");
      if(d<0)
        {
          P=(*it);
        }
      else
        {
          for(int j=0;j<=d;j++)
            {
              MPOL P1;
              if(j!=I)
                {
                  int power=(*it)[j];
                  if(power!=0)
                    { 
                      P1 = monom_t(1,j,power);
                      P*=P1;
                    }
                }
              else if(j==I && (*it)[I]>0)
              {
                P1=monom_t(1,0,0);
                P*=P1;
              }
            }
          P*=(it->coeff());
        }
      for(int k=0;k<d1;k++) 
        {
          A((*it)[I]+k,k+d2) = A((*it)[I]+k,k+d2)+P;
        }
    }
  return A;
}

//--------------------------------------------------------------------

//pour l'instant elle ne calcule le resultant que losque p1 et p2 sont trois variable (a generaliser)

template<class MPOL>
MPOL  Resultant(const MPOL & p1, const MPOL & p2,int I)
{
  typedef typename MPOL::monom_t monom_t;
  MatrDse<MPOL> M=u_sylvester(p1,p2,I);
  return Det(M);
}

//--------------------------------------------------------------------
template<class C, class O, class R> struct MPol;
//--------------------------------------------------------------------
template<class C, class O, class R>
C  Subresultant(const MPol<C,O,R> & p1, const MPol<C,O,R> & p2,int I, int Ord,int Deg, C x,C y) 
{
  int d1= degree(p1,I), d2=degree(p2,I);
  int taille = d1+d2-2*Ord;
  if (taille <1)
    {
      std::cerr<<"unable computing Subresultant"<<std::endl; //a voir proprement
      return 0.;
    }
  else
    {
      MatrDse<MPol<C,O,R> > A(d1+d2-2*Ord,d1+d2-Ord);
      typedef typename MPol<C,O,R>::monom_t monom_t;
      for(typename MPol<C,O,R>::const_iterator it=p1.begin();it!=p1.end();it++)
        {
          int d=(*it).nvars();
          MPol<C,O,R> P;
          if(d<0)
            {
              P=(*it);
            }
          for(int j=0;j<=d;j++)
            {
              MPol<C,O,R> P1;
              if(j!=I)
                {
                  int power=(*it)[j];
                  if(power>0)
                    {
                      P1= monom_t(1,j,power);
                      P+=P1; 
                    }
                }
              else if(j==I && (*it)[I]>0)
                {
                  P1=monom_t(1,0,0);
                  P+=P1;
                }
              P*=(it->coeff());
            }
          for(int k=0;k<d2-Ord;k++)
            {
              A(k,(*it)[I]+k) +=P;
            }
        }
      for(typename MPol<C,O,R>::const_iterator it=p2.begin();it!=p2.end();it++)
        {
          int d=(*it).nvars();
          MPol<C,O,R> P;
          if(d<0)
            {
              P=(*it);
            }
          else
            {
              for(int j=0;j<=d;j++)
                {
                  MPol<C,O,R> P1;
                  if(j!=I)
                    {
                      int power=(*it)[j];
                      if(power!=0)
                        {
                          P1 = monom_t(1,j,power);
                          P+=P1;
                        }
                    }
                  else if(j==I && (*it)[I]>0)
                    {
                      P1=monom_t(1,0,0);
                      P+=P1;
                    }
                }
              P*=(it->coeff());
            }
          for(int k=d2-Ord;k<d1+d2-2*Ord;k++)
            {
              A(k,k+(*it)[I]-d2+Ord) +=P;
            }
        }
      MatrDse<MPol<C,O,R> > Minor(d1+d2-2*Ord,d1+d2-2*Ord);
      for(int j=0;j<d1+d2-2*Ord-1;j++)
        {
          for(int i=0;i<d1+d2-2*Ord;i++)
            {
              Minor(i,j)=A(i,j);
            }
        }
      for(int i=0;i<d1+d2-2*Ord;i++)
        {
          Minor(i,d1+d2-2*Ord-1)=A(i,d1+d2-Ord-Deg-1);
        }
      MPol<C,O,R> subres=(Det(Berko(),Minor));
      C z=0.;
      int var1=0;
      for(int j=0;j<=2;j++)
        {
          if (j!=I)
            {
              var1=j;break;
            }
        }
   int var2=0;
      for(int j=0;j<=2;j++)
        {
          if (j!=I && j!=var1)
            {
              var2=j;break;
            }
        }
      for(typename MPol<C,O,R>::const_iterator it=subres.begin();it!=subres.end();it++)
        {
          int dx=(*it)[var1];
          int dy=(*it)[var2];
          C z1=pow(x,dx)*pow(y,dy);
          z1*=(*it).coeff();
          z+=z1;
        }
      return(z);
    }
}
//----------------------------------------------------------------------


__END_NAMESPACE_SYNAPS

#endif // SYNAPS_RESULTANT_UNIRESLT_H