(* coll[n] is the collection of all α-optimal trees of size n. First, the initial value is given. *)

coll[0]={{{2,1},{}}};

(* Combines three rooted trees to a single new tree. *)

myprod[{a_,b_},{c_,d_},{e_,f_}]:={a*c*e+b*d*f,a*c*e};

(* Trees that are strictly dominated by others are removed *)

domin[{{a_,b_},_},{{c_,d_},_}]:=(a>c)&&(a+b>c+d);

DelDomin[l_,j_]:=Module[{m={}},If[j<Length[l],Do[If[domin[l[[j]],l[[k]]],m=Append[m,{k}]],{k,j+1,Length[l]}]];Delete[l,m]];

NoDomin[l_]:=Module[{m=Reverse[Sort[l]],n=Reverse[Sort[l]]},Do[m=DelDomin[n,j];n=m,{j,1,Length[l]}];n];

(* This calculates the upper envelope *)

FindEnvelope[l_]:=Module[{m={l[[1]]},n=Delete[l,1],x={0},y},While[Length[n]>0,y=-(m[[-1,1,1]]-n[[1,1,1]])/(m[[-1,1,2]]-n[[1,1,2]]);
Which[y>1,n=Delete[n,1],y≥x[[-1]],x=Append[x,y];
m=Append[m,n[[1]]];n=Delete[n,1],True,m=Delete[m,-1];
x=Delete[x,-1]]];{m,x}];

(* Now, coll[n] can be constructed recursively from the previous lists *)

coll[n_]:=coll[n]=FindEnvelope[NoDomin[Union[Flatten[Table[{myprod[coll[i1][[j1,1]],coll[i2][[j2,1]],coll[n-i1-i2-1][[k,1]]],{{i1,j1},{i2,j2},{n-i1-i2-1,k}}},{i1,0,Floor[(n-1)/3]},{i2,i1,Floor[(n-1-i1)/2]},{j1,1,Length[coll[i1]]},{j2,If[i1==i2,j1,1],Length[coll[i2]]},{k,If[n-i1-i2-1==i2,j2,1],Length[coll[n-i1-i2-1]]}],4]]]][[1]];

(* Finds optimal pairs of rooted trees (optimal trees of size 3n for some n) *)

myprod2[{a_,b_},{c_,d_}]:=a*c+b*d;

optpair[n_]:=Reverse[Sort[Flatten[Table[{myprod2[coll[i1][[j1,1]],coll[n-i1][[j2,1]]],{{i1,j1},{n-i1,j2}}},{i1,0,Floor[n/2]},{j1,1,Length[coll[i1]]},{j2,If[i1==n/2,j1,1],Length[coll[n-i1]]}],2]]][[1]]

(* Here the optimal trees are constructed explicitly *)

<<Combinatorica`

assoctree[n_,j_]:=If[n==0,Star[1],AddEdges[GraphUnion[Star[1],assoctree[coll[n][[j,2,1,1]],coll[n][[j,2,1,2]]],assoctree[coll[n][[j,2,2,1]],coll[n][[j,2,2,2]]],assoctree[coll[n][[j,2,3,1]],coll[n][[j,2,3,2]]]],{{1,2},{1,2+V[assoctree[coll[n][[j,2,1,1]],coll[n][[j,2,1,2]]]]},{1,2+V[assoctree[coll[n][[j,2,1,1]],coll[n][[j,2,1,2]]]]+V[assoctree[coll[n][[j,2,2,1]],coll[n][[j,2,2,2]]]]}}]];

graphs1=Flatten[Table[{assoctree[n,1],AddEdges[GraphUnion[Star[1],assoctree[n,-1]],{{1,2}}],AddEdges[GraphUnion[Star[1],assoctree[optpair[n][[2,1,1]],optpair[n][[2,1,2]]],assoctree[optpair[n][[2,2,1]],optpair[n][[2,2,2]]]],{{1,2},{1,2+V[assoctree[optpair[n][[2,1,1]],optpair[n][[2,1,2]]]]}}]},{n,0,6}]];

(* Finally, the results are checked by comparing them to the trees found by a brute-force search *)

SetDirectory["C:\\Documents and Settings\\swagner\\My Documents\\MerrifieldSimmons"];

<<gopt.m;

graphs2=Map[FromAdjacencyLists,%,{2}];

Table[IsomorphicQ[graphs1[[i]],graphs2[[i,1]]],{i,1,20}]

(* A list of optimal values *)

Flatten[Table[{coll[n][[1,1,1]],coll[n][[-1,1,1]]+coll[n][[-1,1,2]],optpair[n][[1]]},{n,0,20}]]