starting implementation of polynomials

Makefile.am

@@ -1,6 +1,6 @@
 lib_LTLIBRARIES = libla.la
-include_HEADERS = vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h  davidson.h   laerror.h    mat.h          qsort.h             vec.h bisection.h   diis.h       la.h         noncblas.h     smat.h bitvector.h   fourindex.h  la_traits.h  nonclass.h     sparsemat.h sparsesmat.h csrmat.h conjgrad.h    gmres.h      matexp.h     permutation.h  
+include_HEADERS = vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h  davidson.h   laerror.h    mat.h          qsort.h             vec.h bisection.h   diis.h       la.h         noncblas.h     smat.h bitvector.h   fourindex.h  la_traits.h  nonclass.h     sparsemat.h sparsesmat.h csrmat.h conjgrad.h    gmres.h      matexp.h     permutation.h   polynomial.h
-libla_la_SOURCES = quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc  bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc
+libla_la_SOURCES = quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc  bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc
 check_PROGRAMS = t test
 t_SOURCES = t.cc t2.cc 
 test_SOURCES = test.cc

la.h

@@ -44,6 +44,7 @@
 #include "sparsesmat.h"
 #include "csrmat.h"
 #include "vec.h"
+#include "polynomial.h"
 
 using namespace LA;
 typedef NRMat<int> NRIMat;

polynomial.cc

@@ -0,0 +1,44 @@
+/*
+    LA: linear algebra C++ interface library
+    Copyright (C) 2021 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "polynomial.h"
+#include <stdio.h>
+#include <string.h>
+
+
+namespace LA {
+
+
+
+/***************************************************************************//**
+ * forced instantization in the corresponding object file
+ ******************************************************************************/
+template class Polynomial<int>;
+template class Polynomial<float>;
+template class Polynomial<double>;
+
+#define INSTANTIZE(T) \
+
+
+//INSTANTIZE(float)
+
+//INSTANTIZE(double)
+
+
+
+}//namespace

polynomial.h

@@ -0,0 +1,102 @@
+/*
+    LA: linear algebra C++ interface library
+    Copyright (C) 2021 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef _POLYNOMIAL_H
+#define _POLYNOMIAL_H
+
+#include "la_traits.h"
+#include "vec.h"
+
+namespace LA {
+
+template <typename T>
+class Polynomial : public NRVec<T> {
+public:
+	Polynomial():  NRVec<T>() {};
+	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;};
+
+	int degree() const {return NRVec<T>::size()-1;};
+	void resize(const int n, const bool preserve=true) {NRVec<T>::resize(n+1,preserve);}
+	
+	Polynomial& operator+=(const T &a) {NOT_GPU(*this); NRVec<T>::copyonwrite(); (*this)[0]+=a; return *this;}
+	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 
+	//unary- inherited
+	
+	Polynomial& operator+=(const Polynomial &rhs) 
+		{
+		NOT_GPU(*this); NRVec<T>::copyonwrite(); 
+		if(rhs.degree()>degree()) resize(rhs.degree());
+		for(int i=0; i<=rhs.degree(); ++i) (*this)[i] += rhs[i];
+		return *this;
+		}
+
+	Polynomial& operator-=(const Polynomial &rhs) 
+                {
+                NOT_GPU(*this); NRVec<T>::copyonwrite(); 
+                if(rhs.degree()>degree()) resize(rhs.degree());
+                for(int i=0; i<=rhs.degree(); ++i) (*this)[i] -= rhs[i];
+                return *this;
+                }
+	Polynomial operator+(const Polynomial &rhs) const {return Polynomial(*this) += rhs;};
+	Polynomial operator-(const Polynomial &rhs) const {return Polynomial(*this) -= rhs;};
+	Polynomial operator*(const Polynomial &rhs) const //for very long polynomials FFT should be used
+		{
+		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;
+		};
+
+	//@@@@
+	//simplify(threshold)
+	//derivative,integral
+	//division remainder 
+	//gcd, lcm
+	//roots, interpolation ... special only for real->complex - declare only and implent only template specialization  in .cc
+
+};
+
+//this is very general, can be used also for nesting polynomials
+template <typename T, typename C>
+C value(const Polynomial<T> &p, const C &x)
+{
+C sum(x); 
+sum=0; //get matrix dimension if C is a matrix
+for(int i=p.degree(); i>0; --i)
+  		 {
+       		 sum+= p[i];
+       		 sum= sum*x; //not *= for matrices
+       		 }
+sum += p[0];
+return sum;
+}
+
+//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;}
+
+
+}//namespace
+#endif

t.cc

@@ -2146,9 +2146,9 @@ int tot=p.generate_all_lex(printme);
 cout <<"generated "<<tot<<endl;
 }
 
-if(1)
+if(0)
 {
-int n;
+int n,unitary_n;;
 cin >>n >>unitary_n;
 Partition<int> p(n);
 space_dim=0;
@@ -2167,7 +2167,29 @@ cout <<Sn;
 if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
 }
 
-
+if(1)
+{
+int n,m;
+double x;
+cin >>n >>m >>x;
+Polynomial<double> p(n),q(m);
+p.randomize(1.);
+q.randomize(1.);
+p*=10.;
+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 <<r;
+cout <<z;
+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.h

@@ -970,7 +970,6 @@ void NRVec<T>::resize(const int n, const bool preserve)
 #ifdef DEBUG
 	if(n < 0) laerror("illegal dimension in NRVec::resize");
 #endif
-if(preserve && n<nn) laerror("cannot resize to smaller vector and preserve data");
 T *vold=0;
 int nnold=0;
 bool preserved=false;
@@ -1054,24 +1053,26 @@ do_preserve:
 if(!preserve || !preserved) laerror("assertion failed in NRVec::resize");
 // omit this check since we would need to have traits for presently unknown user defined classes
 // if(!LA_traits<T>::is_plaindata()) laerror("do not know how to preserve non-plain data");
-if(nnold>=nn)  laerror("assertion2 failed in NRVec::resize");
+
+int nnmin=nnold;
+if(nn<nnmin) nnmin=nn;
 
 #ifdef CUDALA
 if(location == cpu)
         {
 #endif       
-	for(int i=0; i<nnold; ++i) v[i]=vold[i]; //preserve even non-plain data classes
+	for(int i=0; i<nnmin; ++i) v[i]=vold[i]; //preserve even non-plain data classes
-	memset(v+nnold,0,(nn-nnold)*sizeof(T)); //just zero the new memory
+	if(nn>nnold) memset(v+nnold,0,(nn-nnold)*sizeof(T)); //just zero the new memory
 	if(do_delete) delete[] vold;
 #ifdef CUDALA
         }
 else
         {
 	//!!!works only with plain data
-	cublasSetVector(nnold, sizeof(T), vold, 1, v, 1);
+	cublasSetVector(nnmin, sizeof(T), vold, 1, v, 1);
 	TEST_CUBLAS("cublasSetVector"); 
 	T a(0);
-	smart_gpu_set(nn-nnold, a, v+nnold);
+	if(nn>nnold) smart_gpu_set(nn-nnold, a, v+nnold);
 	if(do_delete) gpufree(vold);
         }
 #endif