diff options
| -rw-r--r-- | .gitignore | 3 | ||||
| -rw-r--r-- | roots.cpp | 208 |
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 @@ -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; } |