diff --git a/Makefile.am b/Makefile.am index 35d386a..c2ce214 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 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 reg.h regsurf.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 reg.cc regsurf.cc +include_HEADERS = version.h miscfunc.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 reg.h regsurf.h +libla_la_SOURCES = simple.cc miscfunc.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 reg.cc regsurf.cc nodist_libla_la_SOURCES = version.cc check_PROGRAMS = t test tX test_reg test_regsurf t_SOURCES = t.cc t2.cc diff --git a/la.h b/la.h index 12353e6..2b85a74 100644 --- a/la.h +++ b/la.h @@ -35,6 +35,7 @@ #include "matexp.h" #include "noncblas.h" #include "nonclass.h" +#include "miscfunc.h" #include "permutation.h" #include "qsort.h" #include "smat.h" diff --git a/miscfunc.cc b/miscfunc.cc new file mode 100644 index 0000000..8ee3e8a --- /dev/null +++ b/miscfunc.cc @@ -0,0 +1,113 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2024 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 "miscfunc.h" +#include "laerror.h" +#include +#include + + +namespace LA { + + +#define MAXFACT 20 +unsigned long long factorial(const int n) +{ +static int ntop=20; +static unsigned long long a[MAXFACT+1]={1,1,2,6,24,120,720,5040, +40320, +362880, +3628800, +39916800ULL, +479001600ULL, +6227020800ULL, +87178291200ULL, +1307674368000ULL, +20922789888000ULL, +355687428096000ULL, +6402373705728000ULL, +121645100408832000ULL, +2432902008176640000ULL}; + +int j; +if (n < 0) laerror("negative argument of factorial"); +if (n > MAXFACT) laerror("overflow in factorial"); + while (ntopn||k<0) return(0); +if(k>n/2) k=n-k; +if(k==0) return(1); +if(k==1) return(n); +int ind=0; +if(n<=MAXBINOM) + { + ind=ibidx(n,k); + if (ibitab[ind]) return ibitab[ind]; + } +/* nonrecurent method used anyway */ +d=n-k; +p=1; +for(;n>d;n--) p *= n; +value=p/factorial(k); +if(n<=MAXBINOM) ibitab[ind]=value; +return value; +} +#undef ibidx +#undef ibidxmaxeven +#undef MAXBINOM + + +unsigned long long longpow(unsigned long long x, int i) +{ +if(i<0) return 0; +unsigned long long y=1; +do + { + if(i&1) y *= x; + i >>= 1; + if(i) x *= x; + } +while(i); +return y ; +} + +}//namespace diff --git a/miscfunc.h b/miscfunc.h new file mode 100644 index 0000000..72f233d --- /dev/null +++ b/miscfunc.h @@ -0,0 +1,35 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2024 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 _MISCFUNC_H +#define _MISCFUNC_H + + + +namespace LA { + + + +extern unsigned long long factorial(const int n); +extern unsigned long long binom(int n, int k); +extern unsigned long long longpow(unsigned long long x, int i); + + +}//namespace +#endif diff --git a/permutation.cc b/permutation.cc index f8f5322..5707eae 100644 --- a/permutation.cc +++ b/permutation.cc @@ -18,6 +18,7 @@ #include "vec.h" #include "permutation.h" +#include "miscfunc.h" #include #include #include @@ -729,75 +730,6 @@ for(k=n-1; k>=0; --k) *this = NRPerm(3,inv); } -#define MAXFACT 20 -PERM_RANK_TYPE factorial(const int n) -{ -static int ntop=20; -static PERM_RANK_TYPE a[MAXFACT+1]={1,1,2,6,24,120,720,5040, -40320, -362880, -3628800, -39916800ULL, -479001600ULL, -6227020800ULL, -87178291200ULL, -1307674368000ULL, -20922789888000ULL, -355687428096000ULL, -6402373705728000ULL, -121645100408832000ULL, -2432902008176640000ULL}; - -int j; -if (n < 0) laerror("negative argument of factorial"); -if (n > MAXFACT) laerror("overflow in factorial"); - while (ntopn||k<0) return(0); -if(k>n/2) k=n-k; -if(k==0) return(1); -if(k==1) return(n); -int ind=0; -if(n<=MAXBINOM) - { - ind=ibidx(n,k); - if (ibitab[ind]) return ibitab[ind]; - } -/* nonrecurent method used anyway */ -d=n-k; -p=1; -for(;n>d;n--) p *= n; -value=p/factorial(k); -if(n<=MAXBINOM) ibitab[ind]=value; -return value; -} -#undef ibidx -#undef ibidxmaxeven -#undef MAXBINOM template @@ -1290,20 +1222,6 @@ this->clear(); for(int i=1; i<=nparts; ++i) (*this)[i]=x[i].size(); } -PERM_RANK_TYPE longpow(PERM_RANK_TYPE x, int i) -{ -if(i<0) return 0; -PERM_RANK_TYPE y=1; -do - { - if(i&1) y *= x; - i >>= 1; - if(i) x *= x; - } -while(i); -return y ; -} - //aux for the recursive procedure static PERM_RANK_TYPE partitioncount; diff --git a/permutation.h b/permutation.h index 6f67540..7408492 100644 --- a/permutation.h +++ b/permutation.h @@ -198,9 +198,6 @@ public: -extern PERM_RANK_TYPE factorial(const int n); -extern PERM_RANK_TYPE binom(int n, int k); -extern PERM_RANK_TYPE longpow(PERM_RANK_TYPE x, int i); //permutations represented in the cycle format template diff --git a/polynomial.cc b/polynomial.cc index 220a424..00a5d40 100644 --- a/polynomial.cc +++ b/polynomial.cc @@ -17,6 +17,7 @@ */ #include "polynomial.h" +#include "miscfunc.h" #include #include #include