(* The algorithm for the Hosoya (Z-) index is completely analogous to that for the Merrifield-Simmons index *)

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

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

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=Sort[l],n=Sort[l]},Do[m=DelDomin[n,j];n=m,{j,1,Length[l]}];n];

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}];

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]];

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

optpair[n_]:=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]]

<<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]]]]}}]];

graphs2=Map[ToAdjacencyLists,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,33}]]];

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

(* Here we check that the optimal trees for the Hosoya index are the same as for the Merrifield-Simmons index if the number of vertices is ≤100. *)

graphs1=<<trees100.m;

graphs1==graphs2

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