continueing on polynomials, fix of NRVec unary minus

polynomial.cc

@@ -23,19 +23,49 @@
 
 namespace LA {
 
+template <typename T>
+void Polynomial<T>::polydiv(const Polynomial &rhs, Polynomial &q, Polynomial &r) const
+{
+if(rhs[rhs.degree()]==0) laerror("division by a polynomial with zero leading coefficient - simplify it first");
+if(rhs.degree()==0) //scalar division
+	{
+	q= *this/rhs[0];
+	r.resize(0,false);
+	r[0]=0;
+	return;
+	}
+int rdegree= rhs.degree();
+int qdegree= degree()-rdegree;
+if(qdegree<0)
+	{
+	q.resize(0,false);
+	q[0]=0;
+	r= *this;
+	return;
+	}
+//general case
+q.resize(qdegree,false);
+r= *this; r.copyonwrite();
+for(int i=degree(); i>=rdegree; --i)
+	{
+	T tmp= r[i]/rhs[rdegree];
+	q[i-rdegree]= tmp;
+	r -= rhs.shifted(i-rdegree)*tmp;
+	}
+r.resize(rhs.degree()-1,true);
+}
 
 
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
 template class Polynomial<int>;
-template class Polynomial<float>;
 template class Polynomial<double>;
+template class Polynomial<std::complex<double> >;
 
 #define INSTANTIZE(T) \
 
 
-//INSTANTIZE(float)
 
 //INSTANTIZE(double)
 

polynomial.h

@@ -29,6 +29,7 @@ template <typename T>
 class Polynomial : public NRVec<T> {
 public:
 	Polynomial():  NRVec<T>() {};
+	Polynomial(const NRVec<T> &v) : NRVec<T>(v) {}; //allow implicit conversion from NRVec
 	Polynomial(const int n) : NRVec<T>(n+1) {};
 	Polynomial(const T &a, const int n) : NRVec<T>(n+1) {NRVec<T>::clear(); (*this)[0]=a;};
 
@@ -39,13 +40,16 @@ public:
 	Polynomial& operator-=(const T &a) {NOT_GPU(*this); NRVec<T>::copyonwrite(); (*this)[0]-=a; return *this;}
 	Polynomial operator+(const T &a) const {return Polynomial(*this) += a;};
 	Polynomial operator-(const T &a) const {return Polynomial(*this) -= a;};
-	//operator *= and * by a scalar inherited 
+	Polynomial operator-() const {return NRVec<T>::operator-();}
-	//unary- inherited
+	Polynomial operator*(const T &a) const {return NRVec<T>::operator*(a);}
+	Polynomial operator/(const T &a) const {return NRVec<T>::operator/(a);}
+	Polynomial& operator*=(const T &a) {NRVec<T>::operator*=(a); return *this;}
+	Polynomial& operator/=(const T &a) {NRVec<T>::operator/=(a); return *this;}
 	
 	Polynomial& operator+=(const Polynomial &rhs) 
 		{
 		NOT_GPU(*this); NRVec<T>::copyonwrite(); 
-		if(rhs.degree()>degree()) resize(rhs.degree());
+		if(rhs.degree()>degree()) resize(rhs.degree(),true);
 		for(int i=0; i<=rhs.degree(); ++i) (*this)[i] += rhs[i];
 		return *this;
 		}
@@ -53,7 +57,7 @@ public:
 	Polynomial& operator-=(const Polynomial &rhs) 
                 {
                 NOT_GPU(*this); NRVec<T>::copyonwrite(); 
-                if(rhs.degree()>degree()) resize(rhs.degree());
+                if(rhs.degree()>degree()) resize(rhs.degree(),true);
                 for(int i=0; i<=rhs.degree(); ++i) (*this)[i] -= rhs[i];
                 return *this;
                 }
@@ -61,16 +65,69 @@ public:
 	Polynomial operator-(const Polynomial &rhs) const {return Polynomial(*this) -= rhs;};
 	Polynomial operator*(const Polynomial &rhs) const //for very long polynomials FFT should be used
 		{
+		NOT_GPU(*this); 
 		Polynomial r(degree()+rhs.degree());
 		r.clear();
 		for(int i=0; i<=rhs.degree(); ++i) for(int j=0; j<=degree(); ++j) r[i+j] += rhs[i]*(*this)[j];
 		return r;
 		};
+	void simplify(const typename LA_traits<T>::normtype thr)
+		{
+		NOT_GPU(*this); 
+		int n=degree();
+		while(n>0 && abs((*this)[n])<thr) --n;
+		resize(n,true);
+		};
+	Polynomial shifted(const int shift) const
+		{
+		if(shift==0) return *this;
+		if(shift>0)
+			{
+			Polynomial r(degree()+shift);
+			for(int i=0; i<shift; ++i) r[i]=0;
+			for(int i=0; i<=degree(); ++i) r[shift+i] = (*this)[i];
+			return r;
+			}
+		else
+			{
+			if(shift+degree()<0)
+				{
+				Polynomial r(0);
+				r[0]=0;
+				return r;
+				}
+			Polynomial r(shift+degree());
+			for(int i= -shift; i<=degree(); ++i) r[shift+i] = (*this)[i];
+			return r;
+			}
+		}
+	Polynomial derivative() const
+		{
+		NOT_GPU(*this); 
+		int n=degree();
+		if(n==0)
+			{
+			Polynomial r(0);
+			r[0]=0;
+			return r;
+			}
+		Polynomial r(n-1);
+		for(int i=1; i<=n; ++i) r[i-1] = (*this)[i]* ((T)i);
+		return r;
+		};
+	Polynomial integral() const
+		{
+		NOT_GPU(*this); 
+		int n=degree();
+		Polynomial r(n+1);
+		r[0]=0;
+		for(int i=0; i<=n; ++i) r[i+1] = (*this)[i]/((T)(i+1));
+		return r;
+		}
+	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;};
 
-	//@@@@
-	//simplify(threshold)
-	//derivative,integral
-	//division remainder 
 	//gcd, lcm
 	//roots, interpolation ... special only for real->complex - declare only and implent only template specialization  in .cc
 

t.cc

@@ -2172,23 +2172,33 @@ if(1)
 int n,m;
 double x;
 cin >>n >>m >>x;
+if(n<m) {int t=n; n=m; m=t;}
 Polynomial<double> p(n),q(m);
 p.randomize(1.);
 q.randomize(1.);
+NRVec<double> qq(q);
+Polynomial<double> qqq(qq);
+Polynomial<double> pp(n); pp.randomize(1.); 
 p*=10.;
+Polynomial<double> a,b;
+p.polydiv(q,a,b);
 Polynomial<double> r=p*q;
 Polynomial<double> z=value(p,q); //p(q(x))
 Polynomial<double> y=value(q,p); 
 cout <<p;
 cout <<q;
+cout <<a;
+cout <<b;
 cout <<r;
 cout <<z;
+cout << "polydiv test "<<(a*q+b -p).norm()<<endl;
 cout << value(p,x)*value(q,x) -value(r,x)<<endl;
 cout << value(p,value(q,x)) -value(z,x)<<endl;
 cout << value(q,value(p,x)) -value(y,x)<<endl;
 NRMat<double> u(5,5);
 u.randomize(1.);
 cout << (value(p,u)*value(q,u) -value(r,u)).norm()<<endl;
+
 }
 
 

vec.cc

@@ -153,10 +153,10 @@ const NRVec<double> NRVec<double>::operator-() const {
 #ifdef CUDALA
 	if(location == cpu){
 #endif
-		cblas_dscal(nn, -1.0, v, 1);
+		cblas_dscal(nn, -1.0, result.v, 1);
 #ifdef CUDALA
 	}else{
-		cublasDscal(nn, -1.0, v, 1);
+		cublasDscal(nn, -1.0, result.v, 1);
 		TEST_CUBLAS("cublasDscal");
 	}
 #endif
@@ -174,10 +174,10 @@ const NRVec<std::complex<double> > NRVec<std::complex<double> >::operator-() con
 #ifdef CUDALA
 	if(location == cpu){
 #endif
-		cblas_zdscal(nn, -1.0, v, 1);
+		cblas_zdscal(nn, -1.0, result.v, 1);
 #ifdef CUDALA
 	}else{
-		cublasZdscal(nn, -1.0, (cuDoubleComplex*)v, 1);
+		cublasZdscal(nn, -1.0, (cuDoubleComplex*)result.v, 1);
 		TEST_CUBLAS("cublasZdscal");
 	}
 #endif

vec.h

@@ -193,10 +193,12 @@ public:
 	inline NRVec& operator+=(const T &a);
 	inline NRVec& operator-=(const T &a);
 	inline NRVec& operator*=(const T &a);
+	inline NRVec& operator/=(const T &a);
 
 	inline const NRVec operator+(const T &a) const;
 	inline const NRVec operator-(const T &a) const;
 	inline const NRVec operator*(const T &a) const;
+	inline const NRVec operator/(const T &a) const;
 
 
 	//! determine the actual value of the reference counter 
@@ -753,6 +755,16 @@ inline NRVec<T> & NRVec<T>::operator*=(const T &a) {
 	return *this;
 }
 
+template <typename T>
+inline NRVec<T> & NRVec<T>::operator/=(const T &a) {
+        NOT_GPU(*this);
+        copyonwrite();
+
+        for(register int i=0; i<nn; ++i) v[i] /= a;
+        return *this;
+}
+
+
 /***************************************************************************//**
  * compute scalar product \f$d\f$ of this vector \f$\vec{x}\f$ of general type <code>T</code> 
  * with given vector \f$\vec{y}\f$ of type <code>T</code> and order \f$N\f$
@@ -858,6 +870,7 @@ inline const NRVec<T> NRVec<T>::unitvector() const {
 NRVECMAT_OPER(Vec,+)
 NRVECMAT_OPER(Vec,-)
 NRVECMAT_OPER(Vec,*)
+NRVECMAT_OPER(Vec,/)
 
 /***************************************************************************//**
  * generate operators involving vector and vector
@@ -1035,13 +1048,13 @@ nn = n;
 		if(location == cpu)
 			{
 #endif
-			if(preserve) {vold=v; do_delete=true;} else delete[] v;
+			if(preserve) {vold=v; preserved= do_delete=true;} else delete[] v;
 			v = new T[nn];
 #ifdef CUDALA
 			}
 		else
 			{
-			 if(preserve) {vold=v; do_delete=true;} else gpufree(v);
+			 if(preserve) {vold=v; d preserved= o_delete=true;} else gpufree(v);
 			v = (T*) gpualloc(nn*sizeof(T));
 			}
 #endif
@@ -1374,6 +1387,9 @@ inline NRVec<double>& NRVec<double>::operator*=(const double &a) {
         return *this;
 }
 
+template<>
+inline NRVec<double>& NRVec<double>::operator/=(const double &a) {return *this *= (1./a);}
+
 /***************************************************************************//**
  * multiplies this complex vector \f$\vec{x}\f$ by a complex scalar value \f$\alpha\f$
  * \f[\vec{x}_i\leftarrow\alpha\vec{x}_i\f]
@@ -1397,6 +1413,9 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator*=(co
         return *this;
 }
 
+template<>
+inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator/=(const std::complex<double> &a) {return *this *= (1./a);}
+
 /***************************************************************************//**
  * computes the inner product of this real vector \f$\vec{x}\f$ with given real vector \f$\vec{y]\f$
  * @param[in] rhs real vector \f$\vec{y}\f$