From e8ca6b583ea60c32809f75c4526fc44ac759afd0 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 9 Jun 2021 22:59:19 +0200 Subject: [PATCH] starting implementation of polynomials --- Makefile.am | 4 +- la.h | 1 + polynomial.cc | 44 ++++++++++++++++++++++ polynomial.h | 102 ++++++++++++++++++++++++++++++++++++++++++++++++++ t.cc | 28 ++++++++++++-- vec.h | 13 ++++--- 6 files changed, 181 insertions(+), 11 deletions(-) create mode 100644 polynomial.cc create mode 100644 polynomial.h diff --git a/Makefile.am b/Makefile.am index 13b42b7..483f126 100644 --- a/Makefile.am +++ b/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 -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 +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 polynomial.cc check_PROGRAMS = t test t_SOURCES = t.cc t2.cc test_SOURCES = test.cc diff --git a/la.h b/la.h index e15c857..512f009 100644 --- a/la.h +++ b/la.h @@ -44,6 +44,7 @@ #include "sparsesmat.h" #include "csrmat.h" #include "vec.h" +#include "polynomial.h" using namespace LA; typedef NRMat NRIMat; diff --git a/polynomial.cc b/polynomial.cc new file mode 100644 index 0000000..c9c09e5 --- /dev/null +++ b/polynomial.cc @@ -0,0 +1,44 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2021 Jiri Pittner or + + 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 . +*/ + +#include "polynomial.h" +#include +#include + + +namespace LA { + + + +/***************************************************************************//** + * forced instantization in the corresponding object file + ******************************************************************************/ +template class Polynomial; +template class Polynomial; +template class Polynomial; + +#define INSTANTIZE(T) \ + + +//INSTANTIZE(float) + +//INSTANTIZE(double) + + + +}//namespace diff --git a/polynomial.h b/polynomial.h new file mode 100644 index 0000000..fbdca35 --- /dev/null +++ b/polynomial.h @@ -0,0 +1,102 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2021 Jiri Pittner or + + 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 . +*/ + + +#ifndef _POLYNOMIAL_H +#define _POLYNOMIAL_H + +#include "la_traits.h" +#include "vec.h" + +namespace LA { + +template +class Polynomial : public NRVec { +public: + Polynomial(): NRVec() {}; + Polynomial(const int n) : NRVec(n+1) {}; + Polynomial(const T &a, const int n) : NRVec(n+1) {NRVec::clear(); (*this)[0]=a;}; + + int degree() const {return NRVec::size()-1;}; + void resize(const int n, const bool preserve=true) {NRVec::resize(n+1,preserve);} + + Polynomial& operator+=(const T &a) {NOT_GPU(*this); NRVec::copyonwrite(); (*this)[0]+=a; return *this;} + Polynomial& operator-=(const T &a) {NOT_GPU(*this); NRVec::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::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::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 +C value(const Polynomial &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 +inline Polynomial operator+(const T &a, const Polynomial &rhs) {return Polynomial(rhs)+=a;} +template +inline Polynomial operator-(const T &a, const Polynomial &rhs) {return Polynomial(rhs)-=a;} + + +}//namespace +#endif diff --git a/t.cc b/t.cc index 7fb9227..45b4bdb 100644 --- a/t.cc +++ b/t.cc @@ -2146,9 +2146,9 @@ int tot=p.generate_all_lex(printme); cout <<"generated "<>n >>unitary_n; Partition p(n); space_dim=0; @@ -2167,7 +2167,29 @@ cout <>n >>m >>x; +Polynomial p(n),q(m); +p.randomize(1.); +q.randomize(1.); +p*=10.; +Polynomial r=p*q; +Polynomial z=value(p,q); //p(q(x)) +Polynomial y=value(q,p); +cout < u(5,5); +u.randomize(1.); +cout << (value(p,u)*value(q,u) -value(r,u)).norm()<::resize(const int n, const bool preserve) #ifdef DEBUG if(n < 0) laerror("illegal dimension in NRVec::resize"); #endif -if(preserve && n::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(nnnnold) 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