synaps/upol/gcd_approx.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): B. Mourrain, GALAAD, INRIA
 *   $Id: gcd_approx.h,v 1.1 2005/07/11 08:12:26 mourrain Exp $
 ********************************************************************/
#ifndef SYNAPS_UPOL_GCD_APPROX_H
#define SYNAPS_UPOL_GCD_APPROX_H
/********************************************************************/


//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/linalg/lapack.h>
#include <synaps/upol.h>
#include <synaps/upol/resultant.h>
#include <synaps/linalg/Svd.h>
#include <synaps/linalg/Decomp.h>


__BEGIN_NAMESPACE_SYNAPS

//--------------------------------------------------------------------
template<class C>
UPolDse<C> Gcd(const UPolDse<C> & a, const UPolDse<C> & b,
                    double eps)
{
  MatrDse<lapack<double> > M,R;
  M = SylvesterD(a,b, lapack<double>());
  M.transpose();
  int m = M.nbrow();
  VectDse<double> S;
  S = Svd(M);

#ifdef VERBATIM
  // cout <<"\nTranspose Sylvester matrix:\n";  cout << M <<endl;
  std::cout <<"Singular values : "<<S/S[0]<<std::endl;
#endif

  int i; for(i=m-1; i>0 && S[i] < eps*S[0];i--);
  if (i != m-1)
    {
      int d= m-1-i;
      R=Decomp(M,QR());
#ifdef DEBUG
          std::cout<<"QR Decomposition:\n"<<R <<std::endl;
#endif
      UPolDse<C> x(d+1,AsSize());
      for(int i=0;i<=d;i++) x[d-i]=R(m-d-1,m-d-1+i)/R(m-d-1,m-d-1);
      return x;
    }
  else
    {
      return UPolDse<C>(1);
    }
}

__END_NAMESPACE_SYNAPS

/********************************************************************/
#endif // SYNAPS_UPOL_GCD_APPROX_H