Sorry your browser is not supported!

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#