From c6a9a8e456c879a368c6e6d071f03849b6616c29 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 28 Dec 2023 23:48:05 +0100 Subject: [PATCH] create numbers.h and numbers.cc --- Makefile.am | 4 +- bitvector.cc | 2 +- bitvector.h | 3 +- la.h | 1 + numbers.cc | 137 +++++++++++++++++++++++++++++++++++++++++++++++++++ numbers.h | 67 +++++++++++++++++++++++++ t.cc | 9 +++- 7 files changed, 218 insertions(+), 5 deletions(-) create mode 100644 numbers.cc create mode 100644 numbers.h diff --git a/Makefile.am b/Makefile.am index e86bb16..293c321 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,6 +1,6 @@ lib_LTLIBRARIES = libla.la -include_HEADERS = version.h simple.h 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 la_random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h -libla_la_SOURCES = simple.cc 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 contfrac.cc graph.cc la_random.cc +include_HEADERS = version.h simple.h 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 numbers.h bitvector.h fourindex.h la_traits.h la_random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h +libla_la_SOURCES = simple.cc quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc numbers.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc la_random.cc nodist_libla_la_SOURCES = version.cc check_PROGRAMS = t test tX t_SOURCES = t.cc t2.cc diff --git a/bitvector.cc b/bitvector.cc index afd8d7f..ca630ee 100644 --- a/bitvector.cc +++ b/bitvector.cc @@ -1,6 +1,6 @@ /* LA: linear algebra C++ interface library - Copyright (C) 2008 Jiri Pittner or + Copyright (C) 2008-2023 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 diff --git a/bitvector.h b/bitvector.h index 7edf0b6..ae911a0 100644 --- a/bitvector.h +++ b/bitvector.h @@ -1,6 +1,6 @@ /* LA: linear algebra C++ interface library - Copyright (C) 2008 Jiri Pittner or + Copyright (C) 2008-2023 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 @@ -20,6 +20,7 @@ #define _BITVECTOR_H_ #include "vec.h" +#include "numbers.h" #include diff --git a/la.h b/la.h index fc21e1c..12353e6 100644 --- a/la.h +++ b/la.h @@ -24,6 +24,7 @@ #include "laerror.h" #include "auxstorage.h" #include "bisection.h" +#include "numbers.h" #include "bitvector.h" #include "conjgrad.h" #include "davidson.h" diff --git a/numbers.cc b/numbers.cc new file mode 100644 index 0000000..ac8c217 --- /dev/null +++ b/numbers.cc @@ -0,0 +1,137 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2008-2023 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 "numbers.h" + + +namespace LA { + +template +N primefactor(const N &x) +{ +N i,t; + if ( x <= 2 ) return x; + if ( !(x & 1) ) return 2; + i = 3; + for (;;) { + t = x / i; + if ( t < i ) break; + if ( t * i == x ) + return i; + i += 2; + } + return x; +} + + +template +FACTORIZATION factorization(const N &x) +{ +FACTORIZATION f; +N y=x; +N last=0; +while(y>1) + { + N z=primefactor(y); + if(z!=last) + { + std::pair p; + p.first=z; + p.second=1; + f.push_back(p); + } + else + ++f.back().second; + last=z; + y/=z; + } +return f; +} + + +template +N nextprime(N x) +{ +if ( x < 2 ) return 2; +if ( !(x & 1) ) x--; +while ( !isprime(x += 2) ); +return x; +} + + +template +std::ostream & operator<<(std::ostream &s, const FACTORIZATION &x) +{ +for(auto p=x.begin(); p!=x.end(); ++p) s<<"("<first<<","<second<<") "; +return s; +} + + +template +N eulerphi(const FACTORIZATION &x) +{ +N e=1; +for(auto p=x.begin(); p!=x.end(); ++p) + { + e *= (p->first-1); + if(p->second>1) e*= pow(p->first, p->second-1); + } +return e; +} + + +template +N pow(const N &x, N i) +{ +if(i==0) return 1; +if(i==1) return x; +N y,z; +z=x; +while(!(i&1)) + { + z = z*z; + i >>= 1; + } +y=z; +while((i >>= 1)/*!=0*/) + { + z = z*z; + if(i&1) y = y*z; + } +return y; +} + + + +//force instantization + +#define INSTANTIZE(N) \ +template N primefactor(const N &x); \ +template FACTORIZATION factorization(const N &x); \ +template N nextprime(N x); \ +template std::ostream & operator<<(std::ostream &s, const FACTORIZATION &x); \ +template N pow(const N &x, N i); \ +template N eulerphi(const FACTORIZATION &f); \ + + + + + +INSTANTIZE(uint64_t) + +}//namespace diff --git a/numbers.h b/numbers.h new file mode 100644 index 0000000..da3a9d7 --- /dev/null +++ b/numbers.h @@ -0,0 +1,67 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2008-2023 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 . +*/ + +//some number-theoretical utilities using naive algorithms suitable for small numbers + +#ifndef _NUMBERS_H_ +#define _NUMBERS_H_ + +#include +#include +#include +#include + + + +namespace LA { + +template +N primefactor(const N &x); + +template +bool isprime(const N &x) {return x>1 && x==primefactor(x);} + +template +bool isoddprime(const N &x) {return x>2 && x==primefactor(x);} + +template +class FACTORIZATION : public std::list > +{ +}; + +template +FACTORIZATION factorization(const N &x); + +template +N nextprime(N x); + +template +std::ostream & operator<<(std::ostream &s, const FACTORIZATION &x); + +template +N eulerphi(const FACTORIZATION &f); + +template +N eulerphi(const N &x) {return eulerphi(factorization(x));} + +template +N pow(const N &x, N i); + + +}//namespace +#endif diff --git a/t.cc b/t.cc index 8ac0a84..a7308c2 100644 --- a/t.cc +++ b/t.cc @@ -27,6 +27,7 @@ #include "contfrac.h" #include "simple.h" #include "graph.h" +#include "numbers.h" using namespace std; @@ -2883,7 +2884,7 @@ v.randomize(); for(int i=0; i>n; +cout <