read "funcarith.txt":

ourVertices:=proc(lvl)
	if lvl>0 then
		RETURN(vertices(lvl));
	else
		RETURN([[0],[1],[2]]);
	end if;
end;

clearFunc:=proc(T,level,func)
	local vm, v;

	vm:=ourVertices(level);
	for v in vm do
		T[v,func]:=0;
	end do;

end;

makeIntVers:=proc(lvl)
        local ver;
        ver:=ourVertices(lvl);
        ver:=subsop(1=NULL,2=NULL,3=NULL,ver);
        RETURN(ver);
end;

makeRing5:=proc(T,lvl,func,triangle)
	local success, nb, direction, vert, nextvert, sign;
	
	clearFunc(T,lvl,func);

	sign:=1;
	success:=0;
	vert:=[op(triangle),0,2];
	nextvert:=[op(triangle),0,1];
	direction:=4;
	nb:=3;	

	vert:=T[vert,2];
	
 	if (not (vert = nextvert)) then  
		T[T[vert,nb-1],func]:=sign;
		sign:=-1*sign;
		success:=1;
	end if;

	while (not (vert = nextvert)) do
		T[T[vert,nb],func]:=sign;
		sign:=-1*sign;
		vert:=T[vert,direction];		
	end do;
	
	
	nextvert:=[op(triangle),1,2];	
	vert:=T[vert,2];

	if (not (vert = nextvert)) then 
                T[T[vert,nb-1],func]:=sign;
                sign:=-1*sign;
        end if;

	while (not (vert = nextvert)) do
                T[T[vert,nb],func]:=sign;
                sign:=-1*sign;
                vert:=T[vert,direction];
        end do;

	nextvert:=[op(triangle),0,2];
        vert:=T[vert,2];
                
        if (not (vert = nextvert)) then 
                T[T[vert,nb-1],func]:=sign;
                sign:=-1*sign;
        end if;

        while (not (vert = nextvert)) do
                T[T[vert,nb],func]:=sign;
                sign:=-1*sign;
                vert:=T[vert,direction];
        end do;	

	return success;

end;


makeEdge5:=proc(T,lvl,func,edgein)	
	local success, edge, vert, sign, nextvert, nb, direction;
	
	clearFunc(T,lvl,func);

	success:=0;	
	edge:=edgein;
	direction:=1;
	nb:=2;

	sign:=-1;

	if (edge = [0,2]) then 
		edge:=[2,0];
	end if;

	vert:=[edge[1]];
	nextvert:=[edge[2]];
	vert:=T[vert,1];	

	if (not (vert = nextvert)) then  
                T[T[vert,nb+1],func]:=sign;
                sign:=-1*sign;
        	success:=1;
	end if;

        while (not (vert = nextvert)) do
                T[T[vert,nb],func]:=sign;
                sign:=-1*sign;
                vert:=T[vert,direction];
        end do;

	return success;
end;

makeFan6:=proc(T,lvl,func, vertex)
	local i;
	clearFunc(T,lvl,func);	
	T[vertex,func]:=2;
	for i from 1 to 4 do
		T[T[vertex,i],func]:=-1;	
	end do;

	T[T[T[vertex,2],3],func]:=1;
	T[T[T[vertex,3],2],func]:=1;

end;

#ls is the list of local eigen functions, the new ones are put in startfunc+n
#finallevel=how far will we extend?
makeAllFan6:=proc(T,lvl,nls2,startfunc,finallevel,lambdas,epsilons)
	local nls,ls,ver,sver,vertex,i,j,possibleEps,eps,firstTime,rememberFunc;
	ls:=nls2;
	ver:=vertices(lvl);
	make_neighbors(T,lvl,ver);
	sver:=ourVertices(lvl-2);
	sver:=subsop(1=NULL,2=NULL,3=NULL,sver);
	i:=0; #we are on the first fan
	#numFansPerVertex:=max(1,2^(finallevel-lvl-1));
	if finallevel>lvl then
		if finallevel-lvl-1>0 then
			possibleEps:=rp([-1,1],finallevel-lvl-1);
		else
			possibleEps:=[[]];
		end if;
		for j from 1 to nops(possibleEps) do
			possibleEps[j]:=[1,op(possibleEps[j])];
		end do;
	else
		possibleEps:=[];
	end if;
	for vertex in sver do
		firstTime:=1;
		for eps in possibleEps do
			makeFan6(T,lvl,startfunc+i,vertex);
			oneGramSchmidt(T,[op(ls),startfunc+i]);
			if firstTime=1 then
				rememberFunc:=startfunc+i;
				firstTime:=0;
			end if;
			lambdas[lvl,startfunc+i]:=6;
			epsilons[startfunc+i]:=eps;
			i:=i+1;
		end do;
		if finallevel=lvl then
			makeFan6(T,lvl,startfunc+i,vertex);
			oneGramSchmidt(T,[op(ls),startfunc+i]);
			if firstTime=1 then
				rememberFunc:=startfunc+i;
				firstTime:=0;
			end if;
			lambdas[lvl,startfunc+i]:=6;
			epsilons[startfunc+i]:=[];
			i:=i+1;
		end if;
		ls:=[op(ls),rememberFunc];
	end do;
	RETURN(i);	
end;

#ls is a list of funcnums to gram schimdt the last with the first
oneGramSchmidt:=proc(T,lvl,ls)
	local i,m,ma,ma2;

	discLaplac(T,lvl,ls[-1],6000);
        funcQuotient(T,lvl,6000,ls[-1],6000);
        ma:=findFuncMaxDifferenceWithoutZero(T,lvl,6000,.00000001):


	#lvl:=nops(T[[0],1])-1;
	#print(lvl);
	for i from 1 to nops(ls)-1 do
		#m:=evalf(discInnerProduct(T,lvl,ls[i],ls[nops(ls)]),20);
		m:=discInnerProduct(T,lvl,ls[i],ls[nops(ls)]);
		funcSumWithConstant(T,lvl,-m,ls[nops(ls)],ls[i],ls[nops(ls)]);

        discLaplac(T,lvl,ls[-1],6000);
        funcQuotient(T,lvl,6000,ls[-1],6000);
        ma:=findFuncMaxDifferenceWithoutZero(T,lvl,6000,.00000001):
	ma2:=discInnerProduct(T,lvl,ls[i],ls[i]);



	end do;
	normalizeFunc(T,lvl,ls[nops(ls)]);

	discLaplac(T,lvl,ls[-1],6000);
        funcQuotient(T,lvl,6000,ls[-1],6000);
        ma:=findFuncMaxDifferenceWithoutZero(T,lvl,6000,.00000001):

end;

makeTri6:=proc(T,lvl,func,ls,whichVect)
	local A,B,C,i;
	clearFunc(T,lvl,func);	
	if whichVect=1 then
		A:=sqrt(2)/12+2*sqrt(11)/33;
		B:=sqrt(2)/12-sqrt(11)/33;
		C:=sqrt(2)/12-sqrt(11)/33;
	elif whichVect=2 then
		A:=sqrt(2)/12-sqrt(11)/33;
		B:=sqrt(2)/12+2*sqrt(11)/33;
		C:=sqrt(2)/12-sqrt(11)/33;
	else
		A:=sqrt(2)/12-sqrt(11)/33;
		C:=sqrt(2)/12+2*sqrt(11)/33;
		B:=sqrt(2)/12-sqrt(11)/33;
	end if;

	T[[op(ls),0,1],func]:=A*0+B*2+C*0;
	T[[op(ls),0,2],func]:=A*0+B*0+C*2;
	T[[op(ls),1,2],func]:=A*2+B*0+C*0;
	T[[op(ls),0,0,1],func]:=A*0+B*(-1)+C*1;
	T[[op(ls),0,0,2],func]:=A*0+B*(1)+C*(-1);
	T[[op(ls),0,1,2],func]:=A*0+B*(-1)+C*(-1);
	T[[op(ls),1,0,1],func]:=A*(1)+B*(-1)+C*(0);
	T[[op(ls),1,0,2],func]:=A*(-1)+B*(-1)+C*(0);
	T[[op(ls),1,1,2],func]:=A*(-1)+B*(1)+C*(0);
	T[[op(ls),2,0,1],func]:=A*(-1)+B*(0)+C*(-1);
	T[[op(ls),2,0,2],func]:=A*(1)+B*(0)+C*(-1);
	T[[op(ls),2,1,2],func]:=A*(-1)+B*(0)+C*(1);

end;

makeTri6Wrapper:=proc(T,startfunc,ls,whichvector,lvl,finallevel,lambdas,epsilons)
	local possibleEps,i,j,eps;
	i:=0;
	if finallevel>lvl then
		if finallevel-lvl-1>0 then
			possibleEps:=p([-1,1],finallevel-lvl-1);
		else
			possibleEps:=[[]];
		end if;
		for j from 1 to nops(possibleEps) do
			possibleEps[j]:=[1,op(possibleEps[j])];
		end do;
	else
		possibleEps:=[];
	end if;
		for eps in possibleEps do
			makeTri6(T,lvl,startfunc+i,ls,whichvector);
			lambdas[lvl,startfunc+i]:=6;
			epsilons[startfunc+i]:=eps;
			i:=i+1;
		end do;
		if finallevel=lvl then
			makeTri6(T,lvl,startfunc+i,ls,whichvector);
			lambdas[lvl,startfunc+i]:=6;
			epsilons[startfunc+i]:=[];
			i:=i+1;
		end if;
	RETURN (i);
end;

makeAllLambda6:=proc(T,lvl,finallevel,startfuncnum,lambdas,epsilons)
	local tris,tri,j,currentFunc,orthoList;
	lambdas:=table();
	orthoList:=[];
	if lvl>2 then
		tris:=p([0,1,2],lvl-2);
	else
		tris:=[[]];
	end if;
	currentFunc:=startfuncnum;
	for tri in tris do
		orthoList:=[op(orthoList),currentFunc];
		for j from 0 to 2 do
			currentFunc:=currentFunc+makeTri6Wrapper(T,currentFunc,tri,j,lvl,finallevel,lambdas,epsilons);
		end do;
	end do;
	currentFunc:=currentFunc+makeAllFan6(T,lvl,orthoList,currentFunc,finallevel,lambdas,epsilons);
	return (currentFunc-startfuncnum);
end;

setFuncAuto:=proc(T,finallevel,startfunc)
	local level;

	if (finallevel < 1) then return 0; end if;

	level:=1;









end;