From 50b2447535ced5e826a44c4a40084bd79f0c7823 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 2 Jan 2024 14:05:28 +0100
Subject: [PATCH] finishing bitvector

---
 bitvector.h | 41 ++++++++++++++++++++++++++++++++++++++---
 numbers.cc  | 24 +++++++++++++++++++++---
 numbers.h   |  4 ++++
 t.cc        |  6 ++++++
 4 files changed, 69 insertions(+), 6 deletions(-)

diff --git a/bitvector.h b/bitvector.h
index e3abdcf..cc838c9 100644
--- a/bitvector.h
+++ b/bitvector.h
@@ -21,8 +21,12 @@
 
 #include "vec.h"
 #include "numbers.h"
+#include "laerror.h"
 #include <stdint.h>
 
+//TODO: if efficiency is requires, make also a monic_bitvector, which will not store the leading 1 explicitly
+//and then the field operations will be done without any resize
+//To avoid confusion this class must NOT be derived from bitvector and have only explicit constructor conversion
 
 namespace LA {
 //compressed storage of large bit vectors
@@ -46,6 +50,19 @@ public:
 	explicit bitvector (const unsigned int n):NRVec<bitvector_block>((n+blockbits-1)/blockbits) {modulo=n%blockbits; memset(v,0,nn*sizeof(bitvector_block));};
 	bitvector (const bitvector_block a, const unsigned int n):NRVec<bitvector_block>(a,(n+blockbits-1)/blockbits) {modulo=n%blockbits;};
 	bitvector(const bitvector &rhs) : NRVec<bitvector_block>(rhs) {modulo=rhs.modulo;};
+	explicit bitvector(const uint8_t *data, const unsigned int n): NRVec<bitvector_block>((n+blockbits-1)/blockbits) 
+		{
+		modulo=n%blockbits;
+		if(endianity()) laerror("not portable to big endian");
+		else memcpy(&v[0],data,(n+7)/8);
+		zero_padding();
+		};      
+	void getdata(uint8_t *data)
+		{
+		if(endianity()) laerror("not portable to big endian");
+                else memcpy(data,&v[0],(size()+7)/8);
+		}
+
 	//operator= seems to be correctly synthetized by the compiler
 	//override dereferencing to address single bits, is however possible
 	//only in the const context (otherwise we would have to define a type which, when assigned to, changes a single bit - possible but probably inefficient)
@@ -57,9 +74,27 @@ public:
 	bitvector_block getblock(const unsigned int i) const {return v[i];}; //integer interpretation
 	void setblock(const unsigned int i, const bitvector_block b) {v[i]=b;};
 	int getblocksize() const {return 8*sizeof(bitvector_block);};
-	void set(const unsigned int i) {v[i/blockbits] |= (1UL<<(i%blockbits));};
-	void reset(const unsigned int i) {v[i/blockbits] &= ~(1UL<<(i%blockbits));};
-	void flip(const unsigned int i) {v[i/blockbits] ^= (1UL<<(i%blockbits));};
+	void set(const unsigned int i) 
+		{
+#ifdef DEBUG
+		if(i>=size()) laerror("bitvector index out of range in");
+#endif
+		v[i/blockbits] |= (1UL<<(i%blockbits));
+		};
+	void reset(const unsigned int i) 
+		{
+#ifdef DEBUG    
+                if(i>=size()) laerror("bitvector index out of range in");
+#endif
+		v[i/blockbits] &= ~(1UL<<(i%blockbits));
+		};
+	void flip(const unsigned int i) 
+		{
+#ifdef DEBUG    
+                if(i>=size()) laerror("bitvector index out of range in");
+#endif
+		v[i/blockbits] ^= (1UL<<(i%blockbits));
+		};
 	const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;};
 	void clear() {copyonwrite(true); memset(v,0,nn*sizeof(bitvector_block));};
 	void fill() {memset(v,0xff,nn*sizeof(bitvector_block));};
diff --git a/numbers.cc b/numbers.cc
index 4d954dc..2d389e6 100644
--- a/numbers.cc
+++ b/numbers.cc
@@ -117,6 +117,23 @@ return y;
 }
 
 
+//avoiding overflow which would occur very soon in (x*y)%m
+template <typename N>
+N multmod(N x, N y, const N &m)
+{
+N sum=0;
+if(y==0) return 0;
+while(x)
+	{
+	if(x&1) sum= (sum+y)%m;
+	x>>=1;
+	y = (y<<1)%m; //still can overflow here but for much larger numbers
+	}
+return sum;
+}
+
+
+
 template <typename N>
 N powmod(const N &x, N i, const N &m)
 {
@@ -126,14 +143,14 @@ N y,z;
 z=x%m;            
 while(!(i&1))   
         {       
-        z = (z*z)%m;
+        z = multmod(z,z,m);
         i >>= 1;
         }       
 y=z;    
 while((i >>= 1)/*!=0*/)
                 {
-                z = (z*z)%m;
-                if(i&1) y = (y*z)%m;
+		z = multmod(z,z,m);
+                if(i&1) y = multmod(y,z,m);
                 }
 return y;
 }
@@ -151,6 +168,7 @@ 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 powmod(const N &x, N i, const N &m); \
+template N multmod(N x, N i, const N &m); \
 template N eulerphi(const FACTORIZATION<N> &f); \
 
 
diff --git a/numbers.h b/numbers.h
index c7249e8..46de580 100644
--- a/numbers.h
+++ b/numbers.h
@@ -65,6 +65,10 @@ N pow(const N &x, N i);
 template <typename N>
 N powmod(const N &x, N i, const N &m);
 
+template <typename N>
+N multmod(N x, N y, const N &m);
+
+
 
 }//namespace
 #endif
diff --git a/t.cc b/t.cc
index c53aa97..a93a2a5 100644
--- a/t.cc
+++ b/t.cc
@@ -2944,6 +2944,12 @@ if(0)
 uint64_t n;
 cin >>n;
 cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl;
+uint64_t a;
+do {
+	a=RANDINT64()%n;
+
+	}while(gcd(a,n)!=1);
+cout <<"E-F test "<<powmod(a,eulerphi(n),n)<<endl;
 }
 
 if(0)