#include <lset.c>
#include <rnd/frnd.c>

char nonconv; int fq[5001];

void findroots_inner(Complex *c,const unsigned o,LSet *pr)
{
	Complex r;
	if (o==1)
	{
		r=-c[0]/c[1];
		LSet_add(pr,&r);
		return;
	}
	int n; Complex f,d,p,or;
	r=frnd(2)-1+I*(frnd(2)-1);
int i=0,j=0; // Complex h[1000];
	do
	{
if (j==500) {r=frnd(2)-1+I*(frnd(2)-1); j=0;} else j++;
if (i>=5000) {nonconv=1; break;}
/*{
	FILE *out=fopen("5000iters.log","at");
	fprintf(out,"-----\n");
	//for (i=0;i<1000;i++) fprintf(out,"h[%d]=%lg+%lgi\n",i,h[i].re,h[i].im);
	fclose(out);
	break;
}*/
//else h[i]=r;
i++;
		or=r; f=0; d=0; p=1;
		for (n=0;n<o;n++,p*=r)
		{
			f+=p*c[n];
			d+=p*c[n+1]*(n+1);
		}
		f+=p*c[o];
		r-=f/d;
	}
	while (modsquared(r-or)>1e-20);
fq[i]++;
	LSet_add(pr,&r);
	for (n=o;n>0;n--) c[n-1]+=r*c[n];
	for (n=0;n<o;n++) c[n]=c[n+1];
	findroots_inner(c,o-1,pr);
}

Complex *findroots(Complex *c,const unsigned o)
{ // c[0] to c[o] are coeffs of 1 to x^o; c is destroyed, return value is created
	LSet r=LSet(Complex);
	findroots_inner(c,o,&r);
	free(c);
	return r.a;
}

#include <graphics.c>
#include <rnd/eithertime.c>
#include <rnd/sq.c>
#include <rnd/Mini.c>

GLuint othertex(const unsigned sz)
{
	GLuint ret; glGenTextures(1,&ret);
	glBindTexture(GL_TEXTURE_2D,ret);
	glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_LINEAR_MIPMAP_LINEAR);
	glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);
	//aniso();
	int n,x,y,xs=sz,ys=sz;
	unsigned char *td=malloc(xs*ys*3); float f;
	for (y=ys-1;y>=0;y--) for (x=xs-1;x>=0;x--)
	{
		n=(y*xs+x)*3;
		f=sq((float)sz/2)/(1+sq((float)x-xs/2)+sq((float)y-ys/2));
		td[n]=td[n+1]=td[n+2]=Mini(0xFF,f);
	}
	gluBuild2DMipmaps(GL_TEXTURE_2D,3,xs,ys,GL_RGB,GL_UNSIGNED_BYTE,td);
	free(td);
	return ret;
}

void putblob(const float x,const float y,const float r)
{
	glTexCoord2f(1,1); glVertex2f(x+r*16,y+r*16);
	glTexCoord2f(1,0); glVertex2f(x+r*16,y-r*16);
	glTexCoord2f(0,0); glVertex2f(x-r*16,y-r*16);
	glTexCoord2f(0,1); glVertex2f(x-r*16,y+r*16);
}

typedef struct {double x,y; int h,o;} Point;

LSet precalc(const int maxh)
{
	LSet ret=LSet(Point); Point p;
	int h,i,j,k,nz,l,sp;
for (i=0;i<=5000;i++) fq[i]=0;
	int temps=0,eqns=0,roots=0;
	for (h=2;h<=maxh;h++) // Complexity measure sum(|c_n|+1)
	{
		p.h=h;
		int *t=malloc(h*sizeof(int));
		for (i=(1<<(h-1))-1;i>=0;i-=2) // 2 step stops t[k-1] being zero
		{
			t[0]=0;
			for (j=h-2,k=0;j>=0;j--)
				if ((i>>j)&1) t[k]++; else {k++; t[k]=0;}
			temps++;
			if (k==0) continue; // k is the order
			p.o=k;
			//p.o=t[k];
			nz=0;
			for (j=k;j>=0;j--) if (t[j]!=0) nz++;
			for (j=(1<<(nz-1))-1;j>=0;j--) // Signs loop
			{
				Complex *c=malloc((k+1)*sizeof(Complex));
				for (l=k,sp=1;l>=0;l--)
					if (t[l]==0 || l==k) c[l]=t[l];
					else {c[l]=(j&sp?t[l]:-t[l]); sp<<=1;}
				eqns++;
					nonconv=0;
Complex *cc=malloc((k+1)*sizeof(Complex)); memcpy(cc,c,(k+1)*sizeof(Complex));
				c=findroots(c,k);
					if (!nonconv)
				for (l=k-1;l>=0;l--)
				{
					roots++;
					p.x=c[l].re; p.y=c[l].im;
					LSet_add(&ret,&p);
				}
					else
				{
					FILE *out=fopen("nonconv.log","at");
					for (l=k;l>=0;l--) fprintf(out,"%+lg*z^%d%s",cc[l].re,l,(l?"":"\n"));
					fclose(out);
				}
				free(c);
free(cc);
			}
		}
		free(t);
	}
	FILE *out=fopen("stats.txt","at");
	fprintf(out,"temps=%d eqns=%d roots=%d\n",temps,eqns,roots);
	fclose(out);
out=fopen("histoiters.csv","wt");
for (i=0;i<=5000;i++) fprintf(out,"%d,%d\n",i,fq[i]);
fclose(out);
	return ret;
}

WINMAIN
{
	int n; gl_ortho=1;
	GRAPHICS(0,0,"Algebraic numbers [Stephen Brooks 2010]");
	GLuint tex=othertex(256),list=0;
	double ox=0,oy=0,zoom=yres/5,k1=0.125,k2=0.5;
	SetCursorPos(xres/2,yres/2);
	double ot=eithertime();
	LSet ps=precalc(15);
	LOOP
	{
		double dt=eithertime()-ot; ot=eithertime();
		ox+=(mx-xres/2)/zoom; oy+=(my-yres/2)/zoom;
		if (KEY(VK_O)) ox=oy=0;
		SetCursorPos(xres/2,yres/2);
		if (mb&1) zoom*=exp(dt*3); if (mb&2) zoom*=exp(-dt*3);
		if (KHIT(VK_Z)) {k1*=1.3; glDeleteLists(list,1); list=0;}
		if (KHIT(VK_X)) {k1/=1.3; glDeleteLists(list,1); list=0;}
		if (KHIT(VK_C)) {k2+=0.05; glDeleteLists(list,1); list=0;}
		if (KHIT(VK_V)) {k2-=0.05; glDeleteLists(list,1); list=0;}
		glMatrixMode(GL_MODELVIEW);
		glPushMatrix();
		glScaled(zoom,zoom,zoom);
		glTranslated((xres/2/zoom)-ox,(yres/2/zoom)-oy,0);
		if (!list)
		{
			list=glGenLists(1); glNewList(list,GL_COMPILE_AND_EXECUTE);
		glEnable(GL_BLEND);
		glBlendFunc(GL_ONE,GL_ONE);
		glDisable(GL_DEPTH_TEST);
		glEnable(GL_TEXTURE_2D);
		glBindTexture(GL_TEXTURE_2D,tex);
		glBegin(GL_QUADS);
		Point *p=ps.a;
		for (n=ps.m-1;n>=0;n--)
		{
			switch (p[n].o)
			{
				case 1: glColor3f(1,0,0); break;
				case 2: glColor3f(0,1,0); break;
				case 3: glColor3f(0,0,1); break;
				case 4: glColor3f(0.7,0.7,0); break;
				case 5: glColor3f(1,0.6,0); break;
				case 6: glColor3f(0,1,1); break;
				case 7: glColor3f(1,0,1); break;
				case 8: glColor3f(0.6,0.6,0.6); break;
				default: glColor3f(1,1,1); break;
			}
			putblob(p[n].x,p[n].y,k1*pow(k2,p[n].h-3));
		}
		glEnd();
			ot=eithertime();
			glEndList();
		}
		else if (list) glCallList(list);
		if (KEY(VK_L)) {glDeleteLists(list,1); list=0;}
		if (KEY(VK_CONTROL) && KHIT(VK_S)) screenshotauto();
		glMatrixMode(GL_MODELVIEW);
		glPopMatrix();
		ccl();
	}
}