An algebraic number is a number that can be expressed as the root of a polynomial with integer coefficients. For example, 2 and -2 are roots of the polynomial x^2-4, so they are algebraic. The complex constant i is also algebraic, since it is a root of x^2
+1. So the algebraic numbers are a subset of the complex numbers. A number like sqrt(3)/3 is algebraic, while pi or e are not. This program creates a bunch of polynomials (that's what that ugly bunch of for statements is, creating polynomials like "-1+x^2-3x^4"), finds their roots (using newton's method), then returns them and plots them based on the polynomial's degree and the size of its coefficients.
type complex
a as float
b as float
endtype
type croot
a as float
b as float
csum as integer
order as integer
endtype
dim factors(0) as complex
dim plotlist(-1) as croot
generatePointList()
plotpoints()
sync
wait key
end
function plotpoints()
lock pixels
for n=0 to array count(plotlist())
mul#=-.1*plotlist(n).csum+1
if mul#<1.0/17 then mul#=1.0/17
rval=mul#*255*cos(90*plotlist(n).order)
gval=mul#*255*sin(80*plotlist(n).order)
bval=mul#*255*cos(45*plotlist(n).order+45)
ink rgb(rval,gval,bval),0
circle xunittopixel(plotlist(n).a),yunittopixel(plotlist(n).b),mul#*15
next n
unlock pixels
endfunction
function generatePointList()
maxh=11
counter=0
for h=2 to maxh
csum1=h
dim t(h)
for i=(1<<(h-1))-1 to 0 step -2
t(0)=0
k=0
for j=h-2 to 0 step -1
if (i>>j)&&1
t(k)=t(k)+1
else
inc k
t(k)=0
endif
next j
if k=0 then goto myspot1
order1=k
nz=0
for j=k to 0 step -1
if t(j)<>0 then inc nz
next j
for j=(1<<(nz-1))-1 to 0 step -1
Dim polynomial(k+1)
sp=1
for l=k to 0 step -1
if t(l)=0 or l=k
polynomial(l)=t(l)
else
if j&&sp
polynomial(l)=t(l)
else
polynomial(l)=-t(l)
sp=sp<<1
endif
endif
next l
polynomialfactor()
addfactorstolist(csum1,order1)
undim polynomial()
next j
myspot1:
next i
undim t()
next h
endfunction
function addfactorstolist(csum,order)
for n=1 to array count(factors())
array insert at bottom plotlist()
plotlist(array count(plotlist())).a=factors(n).a
plotlist(array count(plotlist())).b=factors(n).b
plotlist(array count(plotlist())).order=order
plotlist(array count(plotlist())).csum=csum
next n
endfunction
function polynomialfactorprint()
mystr$=""
arrcount=array count(factors())
for n=1 to arrcount
mystr$=mystr$+str$(factors(n).a)
if factors(n).b>=0
mystr$=mystr$+"+"
endif
mystr$=mystr$+str$(factors(n).b)+"i"
if arrcount>0 and n<arrcount
mystr$=mystr$+", "
endif
next n
print mystr$
endfunction
function polynomialprint()
mystr$=""
arrcount=array count(polynomial())
for n=0 to arrcount
if polynomial(n)<>0
mystr$=mystr$+str$(polynomial(n))+"*x^"+str$(n)
if arrcount>0 and n<arrcount-1
mystr$=mystr$+" + "
endif
endif
next n
print mystr$
endfunction
function polynomialfactor()
undim factors()
order=array count(polynomial())
if polynomial(order)=0
while polynomial(order)=0
dec order
endwhile
endif
if order=0 then exitfunction
dim cpoly(order) as complex
for n=0 to order
cpoly(n).a=polynomial(n)
cpoly(n).b=0
next n
dim factors(order)
polynomialfactorinside(order)
undim cpoly()
endfunction
//uses newton's method to find the complex roots of polynomial with coefficients cpoly.
function polynomialfactorinside(order)
//r=a+bi, "root"
//d=a+bi, "derivative"
//f=a+bi, "function"
if order=1
//factors(1)=-cpoly(0)/cpoly(1)
div#=cpoly(1).a*cpoly(1).a+cpoly(1).b*cpoly(1).b
factors(1).a=-(cpoly(0).a*cpoly(1).a+cpoly(0).b*cpoly(1).b)/div#
factors(1).b=-(cpoly(0).b*cpoly(1).a-cpoly(0).a*cpoly(1).b)/div#
exitfunction
endif
r as complex
r.a=rnd(1000)/500.0-1
r.b=rnd(1000)/500.0-1
d as complex
f as complex
p as complex
for i=0 to 500
//d=p(1)+p(2)*2*r+p(3)*2*r^2+ etc.
//f=p(0)+p(1)*z+p(2)*z^2+ etc.
//rnew=r-f/d (newton's method)
p.a=1
p.b=0
d.a=0
d.b=0
f.a=0
f.b=0
for n=0 to order-1 //construct d and f
//f=f+p*cpoly(n)
f.a=f.a+p.a*cpoly(n).a-cpoly(n).b*p.b
f.b=f.b+p.a*cpoly(n).b+cpoly(n).a*p.b
//d=d+p*cpoly(n+1)*(n+1)
d.a=d.a+(p.a*cpoly(n+1).a-cpoly(n+1).b*p.b)*(n+1)
d.b=d.b+(p.a*cpoly(n+1).b+cpoly(n+1).a*p.b)*(n+1)
temp1#=p.a //p*=r
temp2#=p.b
p.a=temp1#*r.a-r.b*temp2#
p.b=temp1#*r.b+r.a*temp2#
next n
f.a=f.a+p.a*cpoly(order).a-cpoly(order).b*p.b
f.b=f.b+p.a*cpoly(order).b+cpoly(order).a*p.b
//rnew=r-f/d
rnew as complex
div#=d.a*d.a+d.b*d.b
rnew.a=r.a-(f.a*d.a+f.b*d.b)/div#
rnew.b=r.b-(f.b*d.a-f.a*d.b)/div#
diff as complex
diff.a=rnew.a-r.a
diff.b=rnew.b-r.b
if diff.a*diff.a+diff.b*diff.b<.0000001
exit
endif
r.a=rnew.a
r.b=rnew.b
next i
factors(order).a=r.a
factors(order).b=r.b
for n=order to 1 step -1
//c[n-1]+=r*c[n]
cpoly(n-1).a=cpoly(n-1).a+r.a*cpoly(n).a-r.b*cpoly(n).b
cpoly(n-1).b=cpoly(n-1).b+r.a*cpoly(n).b+r.b*cpoly(n).a
next n
for n=0 to order-1
cpoly(n)=cpoly(n+1)
next n
polynomialfactorinside(order-1)
endfunction
function xunittopixel(x#)
x2#=x#*screen width()/6+screen width()/2
endfunction x2#
function yunittopixel(y#)
y2#=-y#*screen width()/6+screen height()/2
endfunction y2#
Inspired by this:
http://en.wikipedia.org/wiki/File:Algebraicszoom.png
[edit]
I'll also add: it plots about 16,000 points, each of those the result of a root finding operation, so it takes about 30 seconds to calculate it all. The picture I linked to has about 163,000 points, so that's one of the reasons my version looks worse