aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore3
-rw-r--r--roots.cpp208
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;
}