-- qt.lua,15 (one edge vector)
-- Julia sets via interval cell-mapping (quadtree version)

--require"julia" local f=f

local io=io
local root,exterior
local cx,cy
local Rxmin,Rxmax,Rymin,Rymax=-2.0,2.0,-2.0,2.0
local white=1.0
local black=0.0
local gray=0.5
local N=0
local nE=0
local E={}
local write=io.write

local function output(a1,a2,a3,a4,a5,a6)
	--write(
	--a1 or ""," ",
	--a2 or ""," ",
	--a3 or ""," ",
	--a4 or ""," ",
	--a5 or ""," ",
	--a6 or ""," \n")
end

local function imul(xmin,xmax,ymin,ymax)
	local mm=xmin*ymin
	local mM=xmin*ymax
	local Mm=xmax*ymin
	local MM=xmax*ymax
	local m,M=mm,mm
	if m>mM then m=mM elseif M<mM then M=mM end
	if m>Mm then m=Mm elseif M<Mm then M=Mm end
	if m>MM then m=MM elseif M<MM then M=MM end
	return m,M
end

local function isqr(xmin,xmax)
	local u=xmin*xmin
	local v=xmax*xmax
	if xmin<=0.0 and 0.0<=xmax then
		if u<v then return 0.0,v else return 0.0,u end
	else
		if u<v then return u,v else return v,u end
	end
end

local function f(xmin,xmax,ymin,ymax)
	local x2min,x2max=isqr(xmin,xmax)
	local y2min,y2max=isqr(ymin,ymax)
	local xymin,xymax=imul(xmin,xmax,ymin,ymax)
	return x2min-y2max+cx,x2max-y2min+cx,2.0*xymin+cy,2.0*xymax+cy
end

local function outside(xmin,xmax,ymin,ymax)
	local x,y
	if 0.0<xmin then x=xmin elseif 0.0<xmax then x=0.0 else x=xmax end
	if 0.0<ymin then y=ymin elseif 0.0<ymax then y=0.0 else y=ymax end
	return x^2+y^2>4.0
end

local function inside(xmin,xmax,ymin,ymax)
	return	xmin^2+ymin^2<=4.0 and xmin^2+ymax^2<=4.0 and
		xmax^2+ymin^2<=4.0 and xmax^2+ymax^2<=4.0
end

local function newcell()
	return {nil,nil,nil,nil,color=gray}
end

local function addedge(a,b)
	nE=nE+1
	E[nE]=b
end

local function refine(q)
	if q.color==gray then
		if q[1]==nil then
			q[1]=newcell()
			q[2]=newcell()
			q[3]=newcell()
			q[4]=newcell()
		else
			refine(q[1])
			refine(q[2])
			refine(q[3])
			refine(q[4])
		end
	end
end

local function clip(q,xmin,xmax,ymin,ymax,o,oxmin,oxmax,oymin,oymax)
	local ixmin,ixmax,iymin,iymax
	if xmin>oxmin then ixmin=xmin else ixmin=oxmin end
	if xmax<oxmax then ixmax=xmax else ixmax=oxmax end
	if ixmin>=ixmax then return end
	if ymin>oymin then iymin=ymin else iymin=oymin end
	if ymax<oymax then iymax=ymax else iymax=oymax end
	--if ixmin<=ixmax and iymin<=iymax then
	if iymin<iymax then
		if q[1]==nil then
			addedge(o,q)
		else
			local xmid=(xmin+xmax)/2.0
			local ymid=(ymin+ymax)/2.0
			clip(q[1],xmin,xmid,ymid,ymax,o,oxmin,oxmax,oymin,oymax)
			clip(q[2],xmid,xmax,ymid,ymax,o,oxmin,oxmax,oymin,oymax)
			clip(q[3],xmin,xmid,ymin,ymid,o,oxmin,oxmax,oymin,oymax)
			clip(q[4],xmid,xmax,ymin,ymid,o,oxmin,oxmax,oymin,oymax)
		end
	end
end

local function map(q,xmin,xmax,ymin,ymax)
	--xmin,xmax,ymin,ymax=f(xmin,xmax,ymin,ymax,cx,cy)
	xmin,xmax,ymin,ymax=f(xmin,xmax,ymin,ymax)
	if outside(xmin,xmax,ymin,ymax) then
		q.color=white
	else
		if not inside(xmin,xmax,ymin,ymax) then addedge(q,exterior) end
		clip(root,Rxmin,Rxmax,Rymin,Rymax,q,xmin,xmax,ymin,ymax)
	end
end

local function update(q,xmin,xmax,ymin,ymax)
	if q.color==gray then
		if q[1]==nil then
			local b=nE
			q[2]=nE+1
			map(q,xmin,xmax,ymin,ymax)
			q[3]=nE
		else
			local xmid=(xmin+xmax)/2.0
			local ymid=(ymin+ymax)/2.0
			update(q[1],xmin,xmid,ymid,ymax)
			update(q[2],xmid,xmax,ymid,ymax)
			update(q[3],xmin,xmid,ymin,ymid)
			update(q[4],xmid,xmax,ymin,ymid)
		end
	end
end

local function color(q)
	if q.color==gray then
		if q[1]==nil then
			for i=q[2],q[3] do
				if E[i].color~=white then return end
			end
			q.color=white N=N+1
		else
			color(q[1])
			color(q[2])
			color(q[3])
			color(q[4])
		end
	end
end

local function prewhite(q)
	if q.color==gray then
		if q[1]==nil then
			for i=q[2],q[3] do
				local c=E[i].color
				if c==white or c==-gray then
					q.color=-gray
					N=N+1
					return
				end
			end
		else
			prewhite(q[1])
			prewhite(q[2])
			prewhite(q[3])
			prewhite(q[4])
		end
	end
end

local function recolor(q)
	if q.color==-gray then
		q.color=gray
	elseif q.color==gray then
		if q[1]==nil then
			q.color=black
		else
			recolor(q[1])
			recolor(q[2])
			recolor(q[3])
			recolor(q[4])
		end
	end
end

local function area(q)
	if q[1]==nil then
		if q.color==white then return 0.0,0.0
		elseif q.color==black then return 0.0,1.0
		else return 1.0,0.0 end
	else
		local g1,b1=area(q[1])
		local g2,b2=area(q[2])
		local g3,b3=area(q[3])
		local g4,b4=area(q[4])
		return (g1+g2+g3+g4)/4.0, (b1+b2+b3+b4)/4.0
	end
end

local function colorup(q)
	if q[1]~=nil and q.color==gray then
		local c1=colorup(q[1])
		local c2=colorup(q[2])
		local c3=colorup(q[3])
		local c4=colorup(q[4])
		if c1==c2 and c1==c3 and c1==c4 then
if c1~=gray then
			q[1]=nil; --q[2]=nil; q[3]=nil; q[4]=nil
N=N+1 end
			q.color=c1
		end
	end
	return q.color
end

local function save(q,xmin,ymin,N)
	if q[1]==nil or N==1 then
		output(xmin,ymin,N,q.color)
	else
		N=N/2
		local xmid=xmin+N
		local ymid=ymin+N
		save(q[1],xmin,ymin,N)
		save(q[2],xmid,ymin,N)
		save(q[3],xmin,ymid,N)
		save(q[4],xmid,ymid,N)
	end
end
local function show(p)
	local N=2^10
	-- io.output(p..".box")
	output(N)
	save(root,0,0,N)
	-- io.close()
end

local t0=0
local function memory(s)
	-- local t=os.clock()
	--local dt=string.format("%f",t-t0)
	-- local dt=t-t0
	-- io.stdout:write(s,"\t",dt," sec\t",t," sec\t",math.floor(collectgarbage"count"/1024),"M\n")
	-- t0=t
end

local function do_(f,s)
	local a,b=f(root,Rxmin,Rxmax,Rymin,Rymax)
	memory(s)
	return a,b
end

local function julia(l,a,b)
memory"begin"
	cx=a	cy=b
	root=newcell()
	exterior=newcell()	exterior.color=white
	show(0)
	for i=1,l do print("\nstep",i)
		nE=0
		do_(refine,"refine")
		do_(update,"update")
		repeat
			N=0 color(root,Rxmin,Rxmax,Rymin,Rymax) print("color",N)
		until N==0 memory"color"
		repeat
			N=0 prewhite(root,Rxmin,Rxmax,Rymin,Rymax) print("prewhite",N)
		until N==0 memory"prewhite"
		do_(recolor,"recolor")
		do_(colorup,"colorup")		print("colorup",N)
		local g,b=do_(area,"area")	print("area",g,b,g+b)
		show(i) memory"output"
	print("edges",nE)
	end
end

--julia(14,0.25,0.35)
--julia(14, -.12, .74 )
--julia(14,0,0)
--julia(12,0.25,0)
--julia(9,0.24,0)
--julia(9,0.26,0)
--julia(13,0,1)
--julia(12,-1.1,0)
--julia(9,-0.12,0.5)

-- figures for paper
--julia(14,0,1)
--julia(14,-1,0)
--julia(12,-0.12, 0.64)
--julia(14,-0.12, 0.60)
--julia(14,-0.12, 0.30)

-- julia (level, a, b) -- julia set de c= a + b i
-- julia(10,-0.25, 0.74)
julia(7,-0.25, 0.74)
