synaps/mpol/elim.h

/*********************************************************************
*      This file is part of the source code of SYNAPS.               *
*      V.t TRAN , GALAAD, INRIA                                      *
**********************************************************************
$Id: Elim.h,v 1.1 2005/07/11 08:15:51 mourrain Exp $
\file 
**********************************************************************/
#ifndef SYNAPS_RESULTANT_ELIM_H
#define SYNAPS_RESULTANT_ELIM_H

#include <synaps/init.h>
#include <synaps/linalg.h>
#include <synaps/linalg/SousMatMax.h>
#include <synaps/mpol/macaulay.h>
#include <synaps/mpol/bezoutian.h>

__BEGIN_NAMESPACE_SYNAPS


template<class LL>
int numvar(const LL l)
{
  int nv=0, n=l.size();
  for (int i=0;i<n;i++) {
    int lv = MPOLDST::lvar(l[i]);
    nv = (nv > lv) ? nv : lv;
  }
  return ++nv;
};

template<class LL>
LL transformType(LL h) {
  typedef typename LL::value_type POL;
  typedef typename POL::monom_t MON;

  LL f;
  int i, n = h.size()-1;
  char xi[20];// = new char(10);

  //number of variables
  int nv=numvar(h);

  //recharging values of the variables in x
  /*  for (i=0;i<n;i++) {
    sprintf(xi,"x%d",i); 
    strcat(xi,itoa(i));
    POL::var(i+1,xi);
  }
  */
  //  delete xi;

  for (i=0;i<=n;i++)
    for (int j=0;j<nv;j++)
      h[i] = MPOLDST::swap(h[i],j,j+n);

  POL h0 = h[0];
  for (i=0; i<n; i++) {
    POL p = h0*MON(i,1) - h[i+1];//pn = h0*xi - hi
    f.push_back(p);
  }

  return f;
};


template <class LL>
LL hiddenVar(const LL & f,  const int  l0, const int l1,  LL& f2)
{
  typedef typename LL::value_type POL;
  typedef typename POL::monom_t MON;
  typedef typename POL::coeff_t TT;
 
  LL f1, f3, f4;
  int N=1; 
  int i,j, flag;
  typename LL::const_iterator ll,ll1;
  typename POL::const_iterator pp,pp1;
  
  for (ll=f.begin();ll!=f.end();ll++) {
    POL tmp1 = (*ll);
    for (j=0;j<l0;j++)
      tmp1 = MPOLDST::subs(tmp1,j,1); //Replace the @l0@ first variables by 1

    for (pp=(*ll).begin();pp!=(*ll).end();pp++) {
      POL ct(*pp);
      if (MPOLDST::degree(ct)==0) {
                int perd =1;
                for (pp1=tmp1.begin();pp1!=tmp1.end();pp1++) {
                  POL rc(*pp1);
                  if (MPOLDST::degree(rc)==0){
                        perd = 0;
                        break;
                  }
                }
                if (perd == 1)
                  tmp1+=POL("1");
                break;
      }
    }
    
    POL tmp2("0");
    for (pp=tmp1.begin(); pp!=tmp1.end();pp++) {
      MON mm = (*pp);
      mm.setcoeff(TT(1));  //set coefficients to 1
      tmp2 +=mm;
    }
    f1.push_back(tmp2);
  }
 
  f2.push_back(POL("0"));//value NUL of Macaulay matrix if it exists

  for (ll=f.begin(),ll1=f1.begin();ll!=f.end();ll++,ll1++) {
    POL pol("0");
    for (pp=(*ll1).begin(); pp!=(*ll1).end();pp++) {
      POL ini("0");//Init
      for (typename POL::const_iterator kk=(*ll).begin(); kk!=(*ll).end();kk++) {
        POL tmp(*kk), div(*pp);
        tmp = tmp/div;
        if (tmp!=0) {//divisible
          flag=1;
          for (j=l0;j<l0+l1;j++) {
            POL con(MON(j,1)); 
            con = tmp/con;
            if (con != 0) 
              flag=0;
          }
          if (flag==1)//yes
            ini = ini + tmp;//insert
        }
      }
      f2.push_back(ini); 
      
      MON mm =(*pp);
      mm.setcoeff(TT(N));//Identify the coefficient of monomial mm by an integer N
      pol +=mm;
      N++;
    }
    f3.push_back(pol);
  }
 
  //Shift the variables has x0,..,xl1
  for (ll=f3.begin();ll!=f3.end();ll++) {
    POL tmp = (*ll);
    for (j=0;j<l1;j++)
      tmp=MPOLDST::swap(tmp,j,j+l0);
    f4.push_back(tmp);
  }

  return f4;
};

template<class MAT, class LL>
MAT Macaulay(const LL & f, const int  l0, const int l1)
{
  typedef typename LL::value_type POL;
  LL f1,f2;
  f1 = hiddenVar<LL>(f,l0,l1,f2); //Hide the l0 firsts variables in f2
 
  //Compute the Macaulay matrix with f1.
  MAT M=macaulay<MAT>(f1);
  
  //Find the coefficients according to the hidden variables stored in f2
  int mc;
  for (int i=0;i<M.nbrow(); i++)
    for (int j=0;j<M.nbcol();j++) {
      POL pm = M(i,j);
      if (pm==0)
                mc =0;
      else
                mc = int((pm.begin()->coeff()).get_d());
      matrixof::assigncoeff(M,i,j,f2[mc]);
    }
  
  return M;  
};

template<class MAT, class LL>
MAT Minor(const MAT& mat, const LL & f)
{
  typedef typename LL::value_type POL;
  typedef typename POL::monom_t MON;
  typedef typename POL::order_t O;
  int d[f.size()], n = f.size()-1, nu = 0;
  POL F;

  typename LL::const_iterator fi;int i;
  for(i=0, fi=f.begin(); i<=n; i++, fi++) {
    d[i] = std::max(MPOLDST::degree(*fi),1);
    nu += d[i];
  }
  nu -= n; 

  MON* xdi = new MON[n+1];
  POL* B   = new POL[n+1];

  xdi[0]= MON(n+1,d[0]);
  B[0]  = POL(1);
  for(int j=1; j<=n; j++){
    xdi[j] = MON(j-1,d[j]);
    B[j]=0;
    for(int i=0;i<d[j];i++) 
      B[j]+=POL(MON(j-1,i));
  }

  POL* E= new POL[n+1];
  E[0]=1;

  POL U(1), V(1);
  int mu;
    
  for(int i=n;i>0;i--) {
    mu = nu-d[i];
    for (typename POL::const_iterator j = E[0].begin(); j!=E[0].end(); j++) {
      choose(mu- degree(*j),i,V);
      E[i] += V*(*j);
    }
    E[0]*=B[i];
  } 
  for (int i=1;i<=n;i++)
    E[i] *= xdi[i];

  //Index of the monomial
  std::map<MON,int,O> index;
  int c=0;
  MON m;
  for(int i=0; i<=n; i++)
    for (typename POL::reverse_iterator j = E[i].rbegin(); j!=E[i].rend();j++){
      m = (*j);
      index[m]=c;
      c++;
    }

  //Construction of the dodus monomials
  POL DO("0"); 
  for (int i=0;i<n;i++) {
    for (int j=i+1;j<=n;j++) {
      POL tmp = xdi[i];
      tmp = tmp*xdi[j];
      DO += tmp;//xdi[i]*xdi[j];
    }
  }

  int v=0; 
  std::vector<int> ind;

  for (int i=1;i<=n;i++) //for each E[i]
    for (typename POL::const_iterator ee=E[i].begin();ee!=E[i].end();ee++) {//for each monomial of E[i]
      POL tmp1(*ee);// tmp1*=xdi[i];
      tmp1 *= MON(n+1,nu-degree(*ee));
      for (typename POL::const_iterator pp=DO.begin();pp!=DO.end();pp++) {
                POL tmp2(*pp), tmp3 = tmp1/tmp2;
                if (tmp3!=0) {//pp in tmp1
                  ind.push_back(index[(*ee)]);
                  v++;
                  break;
                }
      }
    }

 //Construction of maximum minor of Macaulay matrix
  MAT Res;
  if (v<2) {//Not exists
    POL p("1");
    matrixof::reserve(Res,1,1,1);
    matrixof::assigncoeff(Res,0,0,p);
  }
  else {
    matrixof::reserve(Res,v,v,1);
    for (int i=0;i<v;i++)
      for (int j=0;j<v;j++)
                matrixof::assigncoeff(Res,i,j,mat(ind[i],ind[j]));
  }

  return Res;
};


//Applicable with QQ number
template<class LL>
typename LL::value_type Elim_Bezout(const LL & f) 
{
  typedef typename LL::value_type POL;
  typedef typename POL::coeff_t TT;
  typedef MatrDse<POL, linalg::rep2d<POL> > MAT;

  int n = f.size();
  
  MAT B = mbezout<MAT>(f,n);
#ifdef DEBUG
  std::cout<<"\nBezout matrix B :\n"<<B;
#endif 
  VectDse<TT> v(n); 
  for (int i=0; i<n;i++)
    v[i] = rand();

  MatrDse<TT, linalg::rep2d<TT> > BS(B.nbrow(), B.nbcol());
  for (int i=0; i<B.nbrow(); i++) 
    for (int j=0; j<B.nbcol(); j++) {
      POL p = B(i,j);
      TT tmp = p(v.rep());
      BS(i,j) = tmp;
    }

  MAT SM = SousMatMax(BS, B);
#ifdef DEBUG
  std::cout<<"\nMinor of B : \n"<<SM<<std::endl;
#endif
  POL result = Det(SM);//Berko(),SM);

  return result;
};

template<class LL>
typename LL::value_type Elim_Macaulay(const LL & f)
{
  typedef typename LL::value_type POL;
  typedef MatrDse<POL, linalg::rep2d<POL> > MAT;

  int nv =numvar(f)-f.size()+1;
  //  LL f = transformType(h);

  int n = f.size();
  MAT M= Macaulay<MAT,LL>(f,n,nv);
#ifdef DEBUG
  std::cout<<"Macaulay matrix M :\n"<<M<<std::endl;
#endif
  LL f1,f2;
  f1 = hiddenVar(f,n,nv,f2);
  MAT m = Minor<MAT>(M,f1);
#ifdef DEBUG
  std::cout<<"Minor matrix m :\n"<<m<<std::endl;
#endif
  POL d = Det(M);
  POL dm = Det(m);
#ifdef DEBUG
  std::cout<<"Determinant of M : "<<d<<"\nDeterminant of m : "<<dm<<std::endl;
#endif
  if (dm!=0) d /=dm;

  return d;
};

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_RESULTANT_ELIM_H