#include <iostream>
#include <fstream>
#include <vector>
#include <tuple>
#include <functional>
#include <thread>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <cstdint>
#include <ctime>
#include <cassert>
#include "lodepng.h"

using namespace std;

using i64 = int64_t;

const double EPSILON=1e-8;


i64 ipow(i64 b,i64 e){
	assert(e>=0);
	if(e==0)return 1;
	if(e==1)return b;
	i64 x=1;
	while(e!=0){
		if(e&1)x*=b;
		b*=b;
		e>>=1;
	}
	return x;
}

int ilog(i64 n){
	assert(n!=0);
	int l=0;
	while(n>>=1)l++;
	return l;
}


class Polynomial;
ostream& operator<<(ostream &os,const Polynomial &poly);

class Polynomial : public vector<double> {
public:
	Polynomial() = default;
	Polynomial(size_t sz):vector<double>(sz){}
	Polynomial(initializer_list<double> init):vector<double>(init){}

	Polynomial d() const {
		Polynomial res(size()-1);
		for(size_t i=1;i<size();i++){
			res[i-1]=operator[](i)*i;
		}
		return res;
	}

	double ev(double x) const {
		if(size()==0)return 0;
		if(size()==1)return operator[](0);
		assert(operator[](size()-1)!=0);
		double res=operator[](0);
		double xn=x;
		for(size_t i=1;i<size();i++){
			res+=operator[](i)*xn;
			xn*=x;
		}
		return res;
	}

	bool odd() const {
		return size()%2==0;
	}
	bool even() const {
		return size()%2==1;
	}

	int degree() const {
		assert(size()>0);
		return size()-1;
	}

	// Returns quotient, *this becomes remainder
	Polynomial div(const Polynomial &p){
		// cerr<<*this<<" / "<<p<<" = ";
		assert(p.size()>0);
		Polynomial quo(size()-p.size()+1);
		for(ssize_t i=size()-p.size();i>=0;i--){
			double q=operator[](i+p.size()-1)/p.back();
			quo[i]=q;
			for(size_t j=0;j<p.size();j++){
				operator[](i+j)-=q*p[j];
			}
		}
		while(size()>1&&abs(back())<EPSILON)pop_back();
		// cerr<<quo<<"  r. "<<*this<<endl;
		return quo;
	}
};

ostream& operator<<(ostream &os,const Polynomial &poly){
	if(poly.degree()==0)return os<<poly[0];
	bool first=true;
	for(size_t i=0;i<poly.size();i++){
		if(poly[i]==0)continue;
		if(first)first=false;
		else os<<" + ";
		if(i==0)os<<poly[i];
		else {
			if(poly[i]==1)os<<"X";
			else if(poly[i]==-1)os<<"-X";
			else os<<poly[i]<<" X";
			if(i!=1)os<<'^'<<i;
		}
	}
	return os;
}

struct Com{
	double r,i;
};

ostream& operator<<(ostream &os,const Com &com){
	return os<<'('<<com.r<<','<<com.i<<')';
}

namespace FR {
	double bisect(const Polynomial &poly,int niters){
		assert(poly.odd());
		double L=-10,R=10;
		if(poly.ev(L)==0)return L;
		if(poly.ev(R)==0)return R;
		while((poly.ev(L)>0)==(poly.ev(R)>0)){
			L-=20;
			R+=20;
		}
		double M=(L+R)/2;
		for(int i=0;i<niters;i++){
			M=(L+R)/2;
			if(poly.ev(M)==0)break;
			if((poly.ev(M)>0)==(poly.ev(L)>0))L=M;
			else R=M;
		}
		return M;
	}

	double newton(const Polynomial &poly){
		assert(poly.odd());
		Polynomial deriv=poly.d();
		double x=0,xn;
		do {
			double d=deriv.ev(x);
			while(d==0){
				x+=0.01;
				d=deriv.ev(x);
			}
			xn=x-poly.ev(x)/d;
		} while(xn-x>EPSILON);
		return xn;
	}

	pair<double,double> bairstow(const Polynomial &poly){
		assert(poly.even()&&poly.degree()>=2);
		const size_t n=poly.degree();
		double u=poly[n-1]/poly[n],v=poly[n-2]/poly[n];
		if(abs(u)<EPSILON)u=1;
		if(abs(v)<EPSILON)v=1;
		for(int iter=0;iter<30;iter++){
			double barr[n+1],farr[n+1];
			barr[n]=barr[n-1]=0;
			farr[n]=farr[n-1]=0;
			for(ssize_t i=n-2;i>=0;i--){
				barr[i]=poly[i+2]-u*barr[i+1]-v*barr[i+2];
				farr[i]=barr[i+2]-u*farr[i+1]-v*farr[i+2];
			}
			double c=poly[1]-u*barr[0]-v*barr[1];
			double d=poly[0]-v*barr[0];
			double g=barr[1]-u*farr[0]-v*farr[1];
			double h=barr[0]-v*farr[1];

			double discr=1/(v*g*g+h*(h-u*g));
			// cerr<<"u="<<u<<" v="<<v<<" g="<<g<<" h="<<h<<" c="<<c<<" d="<<d<<" discr="<<discr<<endl;
			if(abs(c)<EPSILON&&abs(d)<EPSILON){
				break;
			}

			v-=discr*(-g*v*c+(g*u-h)*d);
			u-=discr*(-h*c+g*d);
		}
		// cerr<<"bairstow("<<poly<<") = {"<<u<<','<<v<<'}'<<endl;
		return {u,v};
	}

	pair<Com,Com> quadratic(double u,double v){
		double D=u*u-4*v;
		if(D>=0){
			double sqrtD2=sqrt(D)/2;
			return {{-u/2+sqrtD2,0},{-u/2-sqrtD2,0}};
		} else {
			double sqrtD2=sqrt(-D)/2;
			return {{-u/2,sqrtD2},{-u/2,-sqrtD2}};
		}
	}

	vector<Com> allroots(const Polynomial &poly){
		assert(poly.degree()>=1);
		if(poly.degree()==1){
			return {{-(double)poly[0]/poly[1],0}};
		}
		if(poly.odd()){
			double rt=bisect(poly,30);
			vector<Com> r=allroots(Polynomial(poly).div({-rt,1}));
			r.push_back({rt,0});
			return r;
		}
		pair<double,double> quadr=bairstow(poly);
		vector<Com> r=poly.degree()==2
			? vector<Com>()
			: allroots(Polynomial(poly).div({quadr.second,quadr.first,1}));
		pair<Com,Com> rts=quadratic(quadr.first,quadr.second);
		r.push_back(rts.first);
		r.push_back(rts.second);
		return r;
	}
}

void check(const Com &c){
	assert(!::isnan(c.r)&&!::isnan(c.i));
}
void check(const vector<Com> &v){
	for(const Com &c : v)check(c);
}

namespace Gens {
	class Generator{
	public:
		virtual Polynomial get(i64 n) const = 0;
		virtual i64 total() const = 0;
	};

	template <int N,int D>
	class Christensen : public Generator{
	public:
		Polynomial get(i64 n) const {
			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;
		}

		i64 total() const {
			i64 n=0;
			for(int d=1;d<=D;d++){
				n+=ipow(2*N+1,d);
			}
			return n;
		}
	};

	template <int D>
	class Derbyshire : public Generator{
	public:
		Polynomial get(i64 n) const {
			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;
		}

		i64 total() const {
			return (1<<(D+1))-2;
		}
	};
}


void writeData(const char *fname,int *tally,int width,int height){
	ofstream f(fname);
	assert(f);
	f<<width<<' '<<height<<'\n';
	for(int i=0;i<width*height;i++){
		f<<tally[i]<<'\n';
	}
}

tuple<int*,int,int> readData(const char *fname){
	ifstream f(fname);
	assert(f);
	int width,height;
	f>>width>>height;
	assert(width>0&&width<10000&&height>0&&height<10000);
	int *tally=new int[width*height];
	for(int i=0;i<width*height;i++){
		f>>tally[i];
	}
	return make_tuple(tally,width,height);
}


inline double interpCurve(double x){
	// return atan(30*x);
	// return pow(x,1.0/2);
	// return sqrt(log(18*pow(x,0.7)+1));
	// return x;
	return log(x+1);
}

inline double interp(double x,double clip){
	return min(interpCurve(x)/interpCurve(clip),1.0);
}

void buildImage(unsigned char *image,const int *tally,int width,int height){
	int maxval=0;
	for(int i=0;i<width*height;i++){
		maxval=max(maxval,tally[i]);
	}

	int *cumul=new int[maxval+1];
	memset(cumul,0,(maxval+1)*sizeof(int));
	for(int i=0;i<width*height;i++){
		cumul[tally[i]]++;
	}
	for(int i=1;i<maxval+1;i++){
		cumul[i]+=cumul[i-1];
	}
	assert(cumul[maxval]==width*height);

	for(int i=0;i<width*height;i++){
		// image[3*i+0]=interp((double)min(tally[i],capval)/capval,0.4)*255;
		image[3*i]=pow((double)cumul[tally[i]]/cumul[maxval],3)*255;
		image[3*i+1]=image[3*i+0];
		image[3*i+2]=image[3*i+0];
	}

	delete[] cumul;
}

void writePlot(int *tally,int width,int height){
	unsigned char *image=new unsigned char[width*height*3];
	buildImage(image,tally,width,height);
	int ret=lodepng::encode("out.png",image,width,height,LCT_RGB);
	if(ret){
		cerr<<lodepng_error_text(ret)<<endl;
		exit(1);
	}
	delete[] image;
}


struct Settings{
	double minx,maxx;
	double miny,maxy;
	int width,height;
};


void tallyThreadFunc(const 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;

	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){
		// Gens::Christensen<4,5> gen;
		Gens::Derbyshire<20> gen;
		i64 ntotal=gen.total();
		int ncores=thread::hardware_concurrency();
		cerr<<"Factoring "<<ntotal<<" polynomials..."<<endl;
		cerr<<"Using "<<ncores<<" thread"<<(ncores==1?"":"s")<<endl;
		thread ths[ncores];
		int **tallies=new int*[ncores];
		cerr<<"Spawning threads: ";
		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);
			if(i!=0)cerr<<", ";
			cerr<<i;
		}
		cerr<<endl<<"Joining threads: ";
		tally=new int[S.width*S.height];
		for(int i=0;i<ncores;i++){
			ths[i].join();
			if(i!=0)cerr<<", ";
			cerr<<i;
			for(int j=0;j<S.width*S.height;j++){
				tally[j]+=tallies[i][j];
			}
			delete[] tallies[i];
		}
		cerr<<endl;

		cerr<<"Writing data..."<<endl;
		writeData("out.txt",tally,S.width,S.height);
	} else if(argc==2){
		tie(tally,S.width,S.height)=readData(argv[1]);
	} else {
		cerr<<"Usage: "<<argv[0]<<" [out.txt]"<<endl;
		return 1;
	}

	cerr<<"Writing image..."<<endl;
	writePlot(tally,S.width,S.height);
	delete[] tally;
}