Starting to work on multithreaded code

2 files changed, 107 insertions, 104 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..04ab9e4
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+roots
+*.o
+*.swp
diff --git a/roots.cpp b/roots.cpp
index 617eea3..d0cffe7 100644
--- a/roots.cpp
+++ b/roots.cpp
@@ -3,6 +3,7 @@
 #include <vector>
 #include <tuple>
 #include <functional>
+#include <thread>
 #include <cstring>
 #include <cmath>
 #include <cstdint>
@@ -231,81 +232,62 @@ void check(const vector<Com> &v){
 	for(const Com &c : v)check(c);
 }
 
-template <Polynomial(*F)(i64 n)>
-class Generator{
-	i64 n=0;
-	bool isEnd=false;
-	Polynomial current;
-
-	Generator(bool isEnd):isEnd(isEnd){}
-
-	static Generator endGen;
-
-public:
-	Generator(){
-		operator++();
-	}
-
-	Generator begin() const {
-		return Generator(*this);
-	}
-
-	const Generator& end(){
-		return endGen;
-	}
-
-	bool operator!=(const Generator &g){
-		return isEnd!=g.isEnd;
-	}
-
-	Generator& operator++(){
-		if(!isEnd){
-			current=F(n++);
-			if(current.size()==0)isEnd=true;
-		}
-		return *this;
-	}
-
-	const Polynomial& operator*(){
-		return current;
-	}
-};
-
-template <Polynomial(*F)(i64)>
-Generator<F> Generator<F>::endGen = Generator<F>(true);
-
 namespace Gens {
+	class Generator{
+	public:
+		virtual Polynomial get(i64 n);
+		virtual i64 total();
+	};
+
 	template <int N,int D>
-	Polynomial christensen(i64 n){
-		int d=1;
-		while(d<=D){
-			i64 p=ipow(2*N+1,d);
-			if(n<p)break;
-			n-=p;
-			d++;
+	class Christensen : public Generator{
+	public:
+		Polynomial get(i64 n){
+			int d=1;
+			while(d<=D){
+				i64 p=ipow(2*N+1,d);
+				if(n<p)break;
+				n-=p;
+				d++;
+			}
+			if(d>D)return {};
+			Polynomial poly(d+1);
+			for(int dd=0;dd<d;dd++){
+				poly[dd]=n%(2*N+1)-N;
+				n/=2*N+1;
+			}
+			poly[d]=1;
+			return poly;
 		}
-		if(d>D)return {};
-		Polynomial poly(d+1);
-		for(int dd=0;dd<d;dd++){
-			poly[dd]=n%(2*N+1)-N;
-			n/=2*N+1;
+
+		i64 total(){
+			i64 n=0;
+			for(int d=1;d<=D;d++){
+				n+=ipow(2*N+1,d);
+			}
+			return n;
 		}
-		poly[d]=1;
-		return poly;
-	}
+	};
 
 	template <int D>
-	Polynomial derbyshire(i64 n){
-		n+=2;
-		const int d=ilog(n);
-		if(d>D)return {};
-		Polynomial poly(d+1);
-		for(int dd=0;dd<d;dd++){
-			poly[dd]=n&(1<<dd)?-1:1;
+	class Derbyshire : public Generator{
+	public:
+		Polynomial get(i64 n){
+			n+=2;
+			const int d=ilog(n);
+			if(d>D)return {};
+			Polynomial poly(d+1);
+			for(int dd=0;dd<d;dd++){
+				poly[dd]=n&(1<<dd)?-1:1;
+			}
+			poly[d]=1;
+			return poly;
 		}
-		poly[d]=1;
-		return poly;
-	}
+
+		i64 total(){
+			return (1<<(D+1))-2;
+		}
+	};
 }
 
 
@@ -352,8 +334,8 @@ void buildImage(unsigned char *image,const int *tally,int width,int height){
 		// image[3*i+0]=interp((double)min(tally[i],300)/300,0.4)*255;
 		image[3*i+0]=interp((double)min(tally[i],40)/40,0.4)*255;
 		// image[3*i+0]=interp((double)tally[i],0.4)*255;
-		image[3*i+1]=0;
-		image[3*i+2]=0;
+		image[3*i+1]=image[3*i+0];
+		image[3*i+2]=image[3*i+0];
 	}
 }
 
@@ -369,52 +351,70 @@ void writePlot(int *tally,int width,int height){
 }
 
 
-int main(int argc,char **argv){
+struct Settings{
+	double minx,maxx;
+	double miny,maxy;
+	int width,height;
+};
+
 
+void tallyThreadFunc(Gens::Generator &gen,i64 start,i64 end,int *tally,const Settings &S){
+	for(i64 index=start;index<end;index++){
+		const Polynomial &poly=gen.get(index);
+		vector<Com> r=FR::allroots(poly);
+		for(const Com &c : r){
+			int x=(c.r-S.minx)/(S.maxx-S.minx)*S.width;
+			int y=(c.i-S.miny)/(S.maxy-S.miny)*S.height;
+			if(x<0||x>=S.width||y<0||y>=S.height)continue;
+			tally[S.width*y+x]++;
+		}
+	}
+}
+
+
+int main(int argc,char **argv){
 	int *tally;
-	int width,height;
+
+	Settings S={
+		.minx=0,.maxx=1.5,
+		.miny=0,.maxy=1.5,
+		.width=2000,.height=0
+	};
+	S.height=S.width/(S.maxx-S.minx)*(S.maxy-S.miny);
+	// const double minx=-3,maxx=3,miny=-3,maxy=3;
 
 	if(argc==1){
-		// const double minx=-3,maxx=3,miny=-3,maxy=3;
-		const double minx=0,maxx=1.5,miny=0,maxy=1.5;
-		width=2000;
-		height=width/(maxx-minx)*(maxy-miny);
-
-		tally=new int[width*height]();
-
-		int prevdegree=0;
-		clock_t clk=clock();
-
-		// for(const Polynomial &poly : Generator<Gens::christensen<4,5>>()){
-		for(const Polynomial &poly : Generator<Gens::derbyshire<24>>()){
-			// cerr<<"poly = "<<poly<<endl;
-			int degree=poly.degree();
-			if(degree!=prevdegree){
-				if(prevdegree!=0){
-					clock_t c2=clock();
-					cerr<<'('<<(double)(c2-clk)/CLOCKS_PER_SEC<<"sec)"<<endl;
-					clk=c2;
-				}
-				prevdegree=degree;
-				cerr<<"d="<<degree<<' ';
-			}
-			vector<Com> r=FR::allroots(poly);
-			for(const Com &c : r){
-				int x=(c.r-minx)/(maxx-minx)*width;
-				int y=(c.i-miny)/(maxy-miny)*height;
-				if(x<0||x>=width||y<0||y>=height)continue;
-				tally[width*y+x]++;
+		// Gens::Christensen<4,5> gen;
+		Gens::Derbyshire<15> gen;
+		i64 ntotal=gen.total();
+		int ncores=thread::hardware_concurrency();
+		cerr<<"Factoring "<<ntotal<<" polynomials..."<<endl;
+		cerr<<"Using "<<ncores<<" threads"<<endl;
+		thread ths[ncores];
+		int **tallies=new int*[ncores];
+		for(int i=0;i<ncores;i++){
+			tallies[i]=new int[S.width*S.height]();
+			ths[i]=thread(tallyThreadFunc,gen,ntotal*i/ncores,ntotal*(i+1)/ncores,tallies[i],S);
+		}
+		tally=new int[S.width*S.height];
+		for(int i=0;i<ncores;i++){
+			ths[i].join();
+			for(int j=0;j<S.width*S.height;j++){
+				tally[j]+=tallies[i][j];
 			}
+			delete[] tallies[i];
 		}
 
-		writeData("out.txt",tally,width,height);
+		cerr<<"Writing data..."<<endl;
+		writeData("out.txt",tally,S.width,S.height);
 	} else if(argc==2){
-		tie(tally,width,height)=readData(argv[1]);
+		tie(tally,S.width,S.height)=readData(argv[1]);
 	} else {
 		cerr<<"Usage: "<<argv[0]<<" [out.txt]"<<endl;
 		return 1;
 	}
 
-	writePlot(tally,width,height);
+	cerr<<"Writing image..."<<endl;
+	writePlot(tally,S.width,S.height);
 	delete[] tally;
 }