File:Algebraicszoom.png
Original file (1,920 × 1,080 pixels, file size: 2.01 MB, MIME type: image/png)
Captions
Summary
[edit]DescriptionAlgebraicszoom.png |
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
[edit]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
[edit]- 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.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 21:48, 27 March 2010 | 1,920 × 1,080 (2.01 MB) | Stephen J. Brooks (talk | contribs) | {{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 |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
File usage on other wikis
The following other wikis use this file:
- Usage on en.wikipedia.org
- Usage on id.wikipedia.org
- Usage on ja.wikipedia.org
- Usage on pt.wikibooks.org
- Usage on ro.wikipedia.org
- Usage on tr.wikipedia.org
- Usage on yo.wikipedia.org