create numbers.h and numbers.cc

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
+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  bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc la_random.cc
+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 

bitvector.cc

@@ -1,6 +1,6 @@
 /*
     LA: linear algebra C++ interface library
-    Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+    Copyright (C) 2008-2023 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

bitvector.h

@@ -1,6 +1,6 @@
 /*
     LA: linear algebra C++ interface library
-    Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+    Copyright (C) 2008-2023 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
@@ -20,6 +20,7 @@
 #define _BITVECTOR_H_
 
 #include "vec.h"
+#include "numbers.h"
 #include <stdint.h>
 
 

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"

numbers.cc

@@ -0,0 +1,137 @@
+/*
+    LA: linear algebra C++ interface library
+    Copyright (C) 2008-2023 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 "numbers.h"
+
+
+namespace LA {
+
+template<typename N>
+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<typename N>
+FACTORIZATION<N>  factorization(const N &x)
+{
+FACTORIZATION<N> f;
+N y=x;
+N last=0;
+while(y>1)
+	{
+	N z=primefactor(y);
+	if(z!=last)
+		{	
+		std::pair<N,N> p;
+		p.first=z;
+		p.second=1;
+		f.push_back(p);
+		}
+	else
+		++f.back().second;
+	last=z;
+	y/=z;
+	}
+return f;
+}
+
+
+template<typename N>
+N nextprime(N x)
+{
+if ( x < 2 ) return 2;
+if ( !(x & 1) ) x--;
+while ( !isprime(x += 2) );
+return x;
+}
+
+
+template <typename N>
+std::ostream & operator<<(std::ostream &s, const FACTORIZATION<N>  &x) 
+{
+for(auto p=x.begin(); p!=x.end(); ++p) s<<"("<<p->first<<","<<p->second<<") ";
+return s;
+}
+
+
+template <typename N>
+N eulerphi(const FACTORIZATION<N> &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 <typename N>
+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<N>  factorization(const N &x); \
+template N nextprime(N x); \
+template std::ostream & operator<<(std::ostream &s, const FACTORIZATION<N>  &x); \
+template N pow(const N &x, N i); \
+template N eulerphi(const FACTORIZATION<N> &f); \
+
+
+
+
+
+INSTANTIZE(uint64_t)
+
+}//namespace

numbers.h

@@ -0,0 +1,67 @@
+/*
+    LA: linear algebra C++ interface library
+    Copyright (C) 2008-2023 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/>.
+*/
+
+//some number-theoretical utilities using naive algorithms suitable for small numbers
+
+#ifndef _NUMBERS_H_
+#define _NUMBERS_H_
+
+#include <stdint.h>
+#include <list>
+#include <tuple>
+#include <iostream>
+
+
+
+namespace LA {
+
+template<typename N>
+N primefactor(const N &x);
+
+template<typename N>
+bool isprime(const N &x) {return x>1 && x==primefactor(x);}
+
+template<typename N>
+bool isoddprime(const N &x) {return x>2 && x==primefactor(x);}
+
+template<typename N>
+class FACTORIZATION : public  std::list<std::pair<N,N> >
+{
+};
+
+template<typename N>
+FACTORIZATION<N> factorization(const N &x);
+
+template<typename N>
+N nextprime(N x);
+
+template <typename N>
+std::ostream & operator<<(std::ostream &s, const FACTORIZATION<N>  &x);
+
+template <typename N>
+N eulerphi(const FACTORIZATION<N> &f);
+
+template <typename N>
+N eulerphi(const N &x) {return eulerphi(factorization(x));}
+
+template <typename N>
+N pow(const N &x, N i);
+
+
+}//namespace
+#endif

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<v.size(); ++i) cout <<(v<<i)<<endl;
 }
 
-if(1)
+if(0)
 {
 int seed;
 int f=open("/dev/random",O_RDONLY);
@@ -2936,7 +2937,13 @@ if(!(l%u).is_zero()) laerror("error in gcd");
 if(!(l%v).is_zero()) laerror("error in gcd");
 if(!(u%g).is_zero()) laerror("error in gcd");
 if(!(v%g).is_zero()) laerror("error in gcd");
+}
 
+if(1)
+{
+uint64_t n;
+cin >>n;
+cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl;
 }
 
 }