implemented polynomial gcd

polynomial.cc

@@ -19,6 +19,7 @@
 #include "polynomial.h"
 #include <stdio.h>
 #include <string.h>
+#include <math.h>
 
 
 namespace LA {
@@ -77,8 +78,8 @@ NRMat<T> Polynomial<T>::companion() const
 if((*this)[degree()]==(T)0) laerror("zero coefficient at highest degree - simplify first");
 NRMat<T> a(degree(),degree());
 a.clear();
-for(int i=0; i<degree(); ++i) a(degree()-1,i) = -(*this)[i];
-for(int i=0; i<degree()-1; ++i) a(i,i+1) = (*this)[degree()];
+for(int i=0; i<degree(); ++i) a(degree()-1,i) = -(*this)[i]/(*this)[degree()];
+for(int i=0; i<degree()-1; ++i) a(i,i+1) = (T)1; 
 return a;
 }
 
@@ -158,6 +159,49 @@ return x;
 }
 
 
+template <typename T>
+Polynomial<T> poly_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr, const int d)
+{
+Polynomial<T> big,small;
+if(p.degree() < q.degree()) {big=q; small=p;} else {big=p; small=q;}
+small.simplify(thr);
+if(small.degree()==0) return small;
+Polynomial<T> help;
+do      {
+        help=small;
+        small= big%small;
+        big=help;
+	small.simplify(thr);
+        }
+while((d<0 && small.degree() != 0) || (d>=0 && small.degree()>=d));
+return big;
+}
+
+
+template <>
+Polynomial<int> svd_gcd(const Polynomial<int> &p, const Polynomial<int> &q, const typename LA_traits<int>::normtype thr)
+{
+laerror("SVD gcd only for floating point numbers");
+return p;
+}
+
+template <typename T>
+Polynomial<T> svd_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr)
+{
+Polynomial<T> big,small;
+if(p.degree() < q.degree()) {big=q; small=p;} else {big=p; small=q;}
+small.simplify(thr);
+if(small.degree()==0) return small;
+NRMat<T> s = Sylvester(p,q);
+NRMat<T> u(s.nrows(),s.nrows()),v(s.ncols(),s.ncols());
+NRVec<typename LA_traits<T>::normtype> w(s.nrows());
+singular_decomposition(s,&u,w,&v,0);
+int rank=0;
+for(int i=0; i<w.size(); ++i) if(w[i]>thr*::sqrt((double)big.degree()*small.degree())) ++rank;
+int d=big.degree()+small.degree()-rank;
+return poly_gcd(big,small,thr,d);
+}
+
 
 
 /***************************************************************************//**
@@ -170,6 +214,9 @@ template class Polynomial<std::complex<double> >;
 #define INSTANTIZE(T) \
 template NRMat<T> Sylvester(const Polynomial<T> &p, const Polynomial<T> &q); \
 template Polynomial<T> lagrange_interpolation(const NRVec<T> &x, const NRVec<T> &y); \
+template Polynomial<T> poly_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr,const int d); \
+template Polynomial<T> svd_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr); \
+
 
 
 INSTANTIZE(int)

polynomial.h

@@ -75,6 +75,7 @@ public:
 	void simplify(const typename LA_traits<T>::normtype thr=0)
 		{
 		NOT_GPU(*this); 
+		this->copyonwrite();
 		int n=degree();
 		while(n>0 && abs((*this)[n])<=thr) --n;
 		resize(n,true);
@@ -148,8 +149,6 @@ public:
 	NRVec<T> realroots(const typename LA_traits<T>::normtype thr) const;
 	T newton(const T x0, const typename LA_traits<T>::normtype thr=1e-14, const int maxit=1000) const; //solve root from the guess
 
-	//@@@gcd, lcm  euler and svd  
-
 };
 
 //this is very general, can be used also for nesting polynomials
@@ -203,6 +202,17 @@ template <typename T>
 extern Polynomial<T> lagrange_interpolation(const NRVec<T> &x, const NRVec<T> &y);
 
 
+template <typename T>
+extern Polynomial<T> poly_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr=0, const int d= -1);
+
+template <typename T>
+extern Polynomial<T> svd_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr=0);
+
+template <typename T>
+Polynomial<T> poly_lcm(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr=0)
+{
+return p*q/poly_gcd(p,q,thr);
+}
 
 
 }//namespace

t.cc

@@ -2203,18 +2203,23 @@ cout << (value(p,u)*value(q,u) -value(r,u)).norm()<<endl;
 
 }
 
-if(1)
+if(0)
 {
 int n;
 cin >>n ;
 NRVec<double> r(n);
-r.randomize(1.);
+//r.randomize(1.);
+//wilkinson's ill-conditionel polynomial
+for(int i=0; i<n;++i) r[i]=i+1;
 double x0=r[0]*0.8+r[1]*0.2;
 r.sort(0);
 Polynomial<double> p=polyfromroots(r);
 cout <<p;
 cout <<r;
-cout <<p.realroots(1e-10);
+NRVec<double> rr= p.realroots(1e-10);
+rr.resize(n,true);
+cout <<rr;
+cout <<"root error = "<<(r-rr).norm()<<endl;
 double x=p.newton(x0);
 double xdif=1e10;
 for(int i=0; i<n; ++i) if(abs(x-r[i])<xdif) xdif=abs(x-r[i]);
@@ -2237,4 +2242,31 @@ Polynomial<double>pp=q.derivative(2);
 cout<<"test deriv. "<<(pp-p).norm()<<endl;
 }
 
+if(1)
+{
+int n;
+cin >>n;
+NRVec<double> rr(n);
+rr.randomize(1.);
+rr.sort(0);
+if(rr.size()>2) rr[1]=rr[0];//make a degenerate root
+NRVec<double> pr(2*n);
+NRVec<double> qr(2*n);
+pr.randomize(1.);
+qr.randomize(1.);
+for(int i=0; i<n; ++i) {pr[i]=qr[i]=rr[i];}
+Polynomial<double> p=polyfromroots(pr);
+Polynomial<double> q=polyfromroots(qr);
+Polynomial<double> g=poly_gcd(p,q,1e-8);
+cout <<"GCD ="<<g;
+Polynomial<double> gg=svd_gcd(p,q,1e-13);
+cout <<"SVDGCD ="<<gg;
+NRVec<double> rrr=g.realroots(1e-5);
+NRVec<double> rrrr=gg.realroots(1e-5);
+cout <<rr<<rrr<<rrrr;
+cout <<"test gcd "<<(rr-rrr).norm()<<endl;
+cout <<"test svdgcd "<<(rr-rrrr).norm()<<endl;
+}
+
+
 }