Jump to content

File:Algebraicszoom.png

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikipedia, the free encyclopedia

Original file (1,920 × 1,080 pixels, file size: 2.01 MB, MIME type: image/png)

Summary

Description
English: Visualisation of the (countable) field of algebraic numbers in the complex plane. Colours indicate degree of the polynomial the number is a root of (red = linear, i.e. the rationals, green = quadratic, blue = cubic, yellow = quartic...). Points becomes smaller as the integer polynomial coefficients become larger. View shows integers 0,1 and 2 at bottom right, +i near top.
Date
Source I (Stephen J. Brooks (talk)) created this work entirely by myself.
Author Stephen J. Brooks (talk) Source code in C with OpenGL.
Other versions leadingcoeff.png

C source code

Here's the source code. OpenGL graphics stuff is mixed in with maths stuff. The mathematical routines are findroots_inner (arguments given in findroots) and precalc (returns a set of algebraic numbers in the Point structure, x+iy is the value, o is the order of the polynomial that produced them and h is the complexity measure of the polynomial). LSet is just a container object (like Vector<Complex> or Vector<Point> in C++). I is the complex number i. frnd(x) produces a random double-precision number on the interval [0,x). Blocks with FILE *out=fopen(...) are logfiles, can be removed if necessary.

#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();
	}
}

Licensing

I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons
attribution
This file is licensed under the Creative Commons Attribution 3.0 Unported license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

Captions

This image shows algebraic numbers on the complex plane, colored by degree. Red = linear, green = quadratic, blue = cubic, yellow = quartic.

Items portrayed in this file

depicts

23 March 2010

image/png

c64a47c1f6b36bee273e95949394d0309ec88265

2,103,680 byte

1,080 pixel

1,920 pixel

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current21:48, 27 March 2010Thumbnail for version as of 21:48, 27 March 20101,920 × 1,080 (2.01 MB)Stephen J. Brooks{{Information |Description = Visualisation of the (countable) field of algebraic numbers in the complex plane. Colours indicate degree of the polynomial the number is a root of (red = linear, i.e. the rationals, green = quadratic, blue = cubic, yello

Global file usage

The following other wikis use this file: