synaps/msolve/Bezout.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): B. Mourrain, GALAAD, INRIA.
 *   $Id: Bezout.h,v 1.1 2005/07/11 10:03:50 jppavone Exp $
 ********************************************************************/
#ifndef SYNAPS_SOLVE_HIDDEN_H
#define SYNAPS_SOLVE_HIDDEN_H
/********************************************************************/
#include <stdio.h>
#include <synaps/init.h>
#include <synaps/linalg.h>
#include <synaps/lapack.h>
#include <synaps/mpol/resultant/MBezout.h>



__BEGIN_NAMESPACE_SYNAPS

//------------------------------------------------------------------
template <class MAT, class MLAP, class SCAL> inline 
void Expand(const MAT & M, MLAP & ML, typename MAT::size_type ind, SCAL ty)
{
  typedef typename MAT::size_type size_type;
  typedef typename MAT::value_type coeff_t;

  for (size_type i = 0; i < M.nbrow(); i++)
    for (size_type j = 0; j < M.nbcol(); j++)
      {
        if (M(i,j)!=0) ML(i,j) = SCAL(0);
        else
          {
            bool find=false;
            typename coeff_t::const_iterator m=M(i,j).begin();
            for (int l=0; l < M(i,j).size(); l++) {
              if ((m->nvars() == ind) && (m->GetDegree(ind) > 0)) {
                ML(i,j)=m->GetCoeff();
                find=true;
              }
              m++;
            }
            if (!find) ML(i,j)=SCAL(0);
          }
      }
}
//---------------------------------------------------------------------------
template <class POL>
MatrDse<array2d<typename POL::coeff_t> >  MPSolveURes(const list<POL> & f, bool badcond=true )
{
  typedef typename POL::coeff_t scal_t;
  typedef MatrDse<array2d<POL> > MatRes;
  typedef typename MatRes::value_type value_type;
  typedef typename MatRes::size_type size_type;
  typedef typename value_type::monom_t monom_t;
  
  size_type nbrow, nbcol, rank;
  size_type nvar=f.size();
  MatRes ResMat;
  POL s;
  list<POL> L;
  POL f0(1); 

  monom_t m0(1,0,1);

  // We construct the linear form, f0, which will be added to the system
  f0*= m0;
  for (size_type i=1; i <=f.size(); i++) {
    monom_t mi(1,i,1), mj(1,i+nvar,1);
    mi*=mj;
    f0+=mi;
    }

  L.push_back(f0);

  // Here we shift the index of variables of the original system, in order
  // to consider the coefficients of the linear form as variables.

  for (list<POL>::const_iterator p=f.begin(); p!=f.end(); p++)
    {
      POL fi(0);
      for (typename POL::const_iterator m=p->begin(); m!= p->end(); m++) {
        int n = m->nvars()+1;
        if (n!=0) {
          typename monom_t::index_t * t= new (typename monom_t::index_t)[n+nvar];
          for(int j=0; j<n;j++) {t[j]=0; t[j+nvar]=(*m)[j];}
          scal_t coef(m->GetCoeff());
          fi +=POL(monom_t(coef,n+nvar,t));
          delete [] t;
        }
        else 
          fi+=POL(m->GetCoeff());
      }
      L.push_back(fi);
    }

  // Construction of the resultant matrix of the system of type
  // defined by MatRes.
  // The first variables, representing the coefficients of the linear form,
  // are hidden.
 
  POL bzt;
  bzt=MBezout(L,nvar);
  ResMat=DeltaOf(bzt,L.size()-1,L.size()-1,L.size()-1,MatRes());
  cout <<ResMat<<endl;
  nbrow=ResMat.nbrow();  nbcol=ResMat.nbcol();
  
  cout <<nbrow<<" "<<nbcol<<endl;
  // Evaluation of the entries of the matrix in a random point.
  // These evaluations are put in another matrix.

  array1d<scal_t> ui(f.size()+1);
  for (size_type i=0; i <= f.size(); i++) ui[i]=rand()%((i+1)*100);

  MatrDse<array2d<scal_t> > Meval(nbrow, nbcol);
  for (size_type i=0; i<nbrow; i++)
    for (size_type j=0; j<nbcol; j++)
      {
        Meval(i,j)=ResMat(i,j)(ui);
      }
  
  array1d<size_type> permut_row(Meval.nbrow()), 
                     permut_col(Meval.nbcol());

  // Computation of the rank of the evaluation matrix.
  // This rank is also the rank of the Bezoutian matrix.
  // It is firts computed with QR decomposition.
  cout <<Meval<<endl;
  
  rank = Rank(Meval,permut_row,permut_col);
 
  // But if we know in advance that the Bezoutian matrix has a
  // bad condition number, we verify the previous value found by
  // a computation with SVD.
  if (badcond) rank=Rank(Meval);

  cout << " Rang de la Bezoutienne : " << rank << endl;

  // If it is not of maximal rank, extract a maximal minor
  bool perm=(rank!=Max(nbrow,nbcol));
  if (perm) 
    ResMat=Extract(ResMat,permut_row, permut_col, rank);

  MatrDse<lapack<scal_t> > D0(rank,rank),Di(rank,rank),VE(rank,rank);
  VectStd<array1d<scal_t> > VE_V(rank), D0V(rank), DiV(rank);
  array2d<scal_t>* res = new array2d<scal_t>(rank,f.size()-1);

  // Take the coefficient matrices of the variable x0 and x1
  Expand(ResMat, D0, 0, scal_t());
  Expand(ResMat, Di, 1, scal_t());
 
  // Compute generalized eigenvector VE such, (Di- v D0)=0
  VE=Eigenvct(D0,Di, Real(),VE_V);

  // Apply each eigenvector to each coefficient matrix of all hidden variables.
  for (size_type l=0;l<rank;l++) {
    VE_V=Eval(Col(VE,l));
    DiV = Di*VE_V ;
    D0V = D0*VE_V ;
    
    for (size_type i=0; i<rank; i++)
      if (DiV[i]!=0) {
        for (size_type j=0;j<rank;j++)
          if (D0V[j]!=0) (*res)(i,0)=DiV[i]/D0V[j];
      }
  }
  for (size_type k=2; k<=f.size(); k++) {
    Expand(ResMat, Di, k,scal_t());
    
    for (size_type l=0;l<rank;l++) {
      VE_V=Eval(Col(VE,l));
      DiV=Di*VE_V;
      D0V = D0*VE_V;
      for (size_type i=0; i<rank; i++)
        if (DiV[i]!=0) {
          for (size_type j=0;j<rank;j++)
            if (D0V[j]!=0) 
              (*res)(i,k-1)=DiV[i]/D0V[j];
        }
    }
  }
  return MatrDse<array2d<scal_t> >(*res);
}

__END_NAMESPACE_SYNAPS

/********************************************************************/
#endif // SYNAPS_SOLVE_HIDDEN_H