continueing on polynomials

polynomial.cc

@@ -56,6 +56,94 @@ r.resize(rhs.degree()-1,true);
 }
 
 
+template <typename T>
+NRMat<T> Sylvester(const Polynomial<T> &p, const Polynomial<T> &q)
+{
+int nm=p.degree()+q.degree();
+NRMat<T> a(nm,nm);
+a.clear();
+for(int i=0; i<q.degree(); ++i)
+	for(int j=p.degree(); j>=0; --j)
+		a(i,i+p.degree()-j)=p[j];
+for(int i=0; i<p.degree(); ++i)
+	for(int j=q.degree(); j>=0; --j)
+		a(q.degree()+i,i+q.degree()-j)=q[j];
+return a;
+}
+
+template <typename T>
+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()];
+return a;
+}
+
+template<>
+NRVec<typename LA_traits<int>::complextype> Polynomial<int>::roots() const
+{
+laerror("roots not implemented for integer polynomials");
+return NRVec<typename LA_traits<int>::complextype>(1);
+}
+
+template <typename T>
+NRVec<typename LA_traits<T>::complextype> Polynomial<T>::roots() const
+{
+NRMat<T> a=this->companion();
+NRVec<typename LA_traits<T>::complextype> r(degree());
+gdiagonalize(a,r,NULL,NULL);
+return r;
+}
+
+template <typename T>
+NRVec<T> Polynomial<T>::realroots(const typename LA_traits<T>::normtype thr) const
+{
+NRVec<typename LA_traits<T>::complextype> r = roots();
+NRVec<T> rr(degree());
+int nr=0;
+for(int i=0; i<degree(); ++i)
+	{
+	if(abs(r[i].imag())<thr)
+		{
+		rr[nr++] = r[i].real();
+		}
+	}
+rr.resize(nr,true);
+rr.sort();
+return rr;
+}
+
+template <typename T>
+Polynomial<T> lagrange_interpolation(const NRVec<T> &x, const NRVec<T> &y)
+{
+if(x.size()!=y.size()) laerror("vectors of different length in lagrange_interpolation");
+if(x.size()<1) laerror("empty vector in lagrange_interpolation");
+if(x.size()==1)
+        {
+        Polynomial<T> p(0);
+        p[0]=y[0];
+        return p;
+        }
+int n=x.size()-1;
+Polynomial<T> p(n);
+p.clear();
+for(int i=0; i<=n; ++i)
+        {
+	T prod=(T)1;
+	for(int j=0; j<=n; ++j) if(j!=i) prod *= (x[i]-x[j]);
+	if(prod==(T)0) laerror("repeated x-value in lagrange_interpolation");
+	Polynomial<T> tmp=polyfromroots(x,i);
+	p += tmp * y[i] / prod;
+        }
+return p;
+}
+
+
+
+
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
@@ -64,10 +152,13 @@ template class Polynomial<double>;
 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); \
 
 
-
-//INSTANTIZE(double)
+INSTANTIZE(int)
+INSTANTIZE(double)
+INSTANTIZE(std::complex<double>)
 
 
 

polynomial.h

@@ -22,6 +22,7 @@
 
 #include "la_traits.h"
 #include "vec.h"
+#include "nonclass.h"
 
 namespace LA {
 
@@ -78,6 +79,7 @@ public:
 		while(n>0 && abs((*this)[n])<thr) --n;
 		resize(n,true);
 		};
+	void normalize() {if((*this)[degree()]==(T)0) laerror("zero coefficient at highest degree - simplify first"); *this /= (*this)[degree()];};
 	Polynomial shifted(const int shift) const
 		{
 		if(shift==0) return *this;
@@ -127,9 +129,11 @@ public:
 	void polydiv(const Polynomial &rhs, Polynomial &q, Polynomial &r) const;
 	Polynomial operator/(const Polynomial &rhs) const {Polynomial q,r; polydiv(rhs,q,r); return q;};
 	Polynomial operator%(const Polynomial &rhs) const {Polynomial q,r; polydiv(rhs,q,r); return r;};
+	NRMat<T> companion() const; //matrix which has this characteristic polynomial
+	NRVec<typename LA_traits<T>::complextype> roots() const; //implemented for complex<double> and double only
+	NRVec<T> realroots(const typename LA_traits<T>::normtype thr) const;
 
-	//gcd, lcm
-	//roots, interpolation ... special only for real->complex - declare only and implent only template specialization  in .cc
+	//@@@gcd, lcm  euler and svd  
 
 };
 
@@ -148,12 +152,43 @@ sum += p[0];
 return sum;
 }
 
+template <typename T, typename C>
+NRVec<C> values(const Polynomial<T> &p, const NRVec<C> &x)
+{
+NRVec<C> r(x.size());
+for(int i=0; i<x.size(); ++i) r[i]=value(p,x[i]);
+return r;
+}
+
+
 //scalar+-polynomial
 template <typename T>
 inline Polynomial<T> operator+(const T &a, const Polynomial<T> &rhs) {return Polynomial<T>(rhs)+=a;}
 template <typename T>
 inline Polynomial<T> operator-(const T &a, const Polynomial<T> &rhs) {return Polynomial<T>(rhs)-=a;}
 
+//Sylvester matrix
+template <typename T>
+extern NRMat<T> Sylvester(const Polynomial<T> &p, const Polynomial<T> &q);
+
+//polynomial from given roots
+template <typename T>
+Polynomial<T> polyfromroots(const NRVec<T> &roots, const int skip= -1)
+{
+Polynomial<T> p(0);
+p[0]=(T)1;
+for(int i=0; i<roots.size(); ++i) 
+    if(i!=skip)
+	p = p.shifted(1) - p * roots[i];
+return p;
+}
+
+
+template <typename T>
+extern Polynomial<T> lagrange_interpolation(const NRVec<T> &x, const NRVec<T> &y);
+
+
+
 
 }//namespace
 #endif

t.cc

@@ -2167,7 +2167,7 @@ cout <<Sn;
 if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
 }
 
-if(1)
+if(0)
 {
 int n,m;
 double x;
@@ -2187,6 +2187,8 @@ Polynomial<double> z=value(p,q); //p(q(x))
 Polynomial<double> y=value(q,p); 
 cout <<p;
 cout <<q;
+cout <<q.companion();
+cout<<Sylvester(p,q);
 cout <<a;
 cout <<b;
 cout <<r;
@@ -2201,5 +2203,30 @@ cout << (value(p,u)*value(q,u) -value(r,u)).norm()<<endl;
 
 }
 
+if(0)
+{
+int n;
+cin >>n ;
+NRVec<double> r(n);
+r.randomize(1.);
+r.sort(0);
+Polynomial<double> p=polyfromroots(r);
+cout <<p;
+cout <<r;
+cout <<p.realroots(1e-10);
+}
+
+if(1)
+{
+int n;
+cin >>n ;
+NRVec<double> x(n+1),y(n+1);
+x.randomize(1.);
+y.randomize(1.);
+Polynomial<double> p=lagrange_interpolation(x,y);
+cout <<x<<y<<p;
+NRVec<double> yy=values(p,x);
+cout<<"interpolation error= "<<(y-yy).norm()<<endl;
+}
 
 }