00001
00002
00003
00004
00005
00006 #ifndef SYNAPS_UPOL_GCD_APPROX_H
00007 #define SYNAPS_UPOL_GCD_APPROX_H
00008
00009
00010
00011
00012 #include <synaps/init.h>
00013 #include <synaps/linalg/lapack.h>
00014 #include <synaps/upol.h>
00015 #include <synaps/upol/resultant.h>
00016 #include <synaps/linalg/Svd.h>
00017 #include <synaps/linalg/Decomp.h>
00018
00019
00020 __BEGIN_NAMESPACE_SYNAPS
00021
00022
00024 template<class C>
00025 UPolDse<C> Gcd(const UPolDse<C> & a, const UPolDse<C> & b,
00026 double eps)
00027 {
00028 MatrDse<lapack<double> > M,R;
00029 M = SylvesterD(a,b, lapack<double>());
00030 M.transpose();
00031 int m = M.nbrow();
00032 VectDse<double> S;
00033 S = Svd(M);
00034
00035 #ifdef VERBATIM
00036
00037 std::cout <<"Singular values : "<<S/S[0]<<std::endl;
00038 #endif
00039
00040 int i; for(i=m-1; i>0 && S[i] < eps*S[0];i--);
00041 if (i != m-1)
00042 {
00043 int d= m-1-i;
00044 R=Decomp(M,QR());
00045 #ifdef DEBUG
00046 std::cout<<"QR Decomposition:\n"<<R <<std::endl;
00047 #endif
00048 UPolDse<C> x(d+1,AsSize());
00049 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);
00050 return x;
00051 }
00052 else
00053 {
00054 return UPolDse<C>(1);
00055 }
00056 }
00057
00058 __END_NAMESPACE_SYNAPS
00059
00060
00061 #endif // SYNAPS_UPOL_GCD_APPROX_H