/* Important note: This program is written in magma, so you need this software on your computer. It also makes use of the geng program in Brendan McKay's graph toolkit and nauty. You can download these two programs from http://cs.anu.edu.au/people/bdm/nauty/ The following command sets a path variable to point to the authors copy of the program called geng. You need to change it to point to your copy. */ gengpath := "/u/keevash/gtools10/geng"; /* Written by Peter Keevash - January 2002. 1. Purpose of program: find smallest value of a maximal fractional packing of monochromatic triangles in a 2-coloured complete graph. 2. Algorithm (description taken from paper): Initialisation: Choose a starting value for $n$, the number of vertices, and $d$, the `search depth'. Check all graphs on $n$ vertices. Find the $d$ smallest values of $t^*$ that occur. Build a list $L_n$ of all graphs for which $t^*$ is one of the $d-1$ smallest values. (Our lists will always be `reduced': containing no isomorphic or complementary pairs.) Iteration: Check all one vertex extensions of the graphs in $L_n$. By averaging, any other graph on $n+1$ vertices must have $t^*$ at least $\frac{n+1}{n-1}$ times the $d$th smallest value for graphs on $n$ vertices. We refer to this lower bound as the {\em level} at $n+1$. The smallest values below the level that we find by checking these extensions are indeed the smallest values for all graphs on $n+1$ vertices. We find as many of the smallest values as possible up to a maximum of $d$. If there are $d$ values below the level, then we let $L_{n+1}$ be the graphs with the $d-1$ smallest values and procede with the iteration. However, there may not be $d$ values below the level, so then we have to reduce the depth of the search: we let $L_{n+1}$ be all the graphs we have found below the level, and calculate the level for $n+2$ to be $\frac{n+2}{n}$ times the level for $n+1$. Then we procede with the iteration. Termination: There may come a point at which no values are found below the level, at which point nothing further can be said. To get results for larger $n$ it is then necessary to start again with a larger search depth. 3. Operation: the user is prompted to input 3 starting variables - the starting value of N, the search depth d, and the value of N to stop at. 4. Remarks on implementation: The function Fract(G) computes value of maximum fractional packing of triangles in G and its complement. This is a linear programming problem, for which magma has a built-in solver. Unfortunately, it will only return (approximate) real solutions rather than (exact) rational solutions, so it is necessary to use the procedure Rat, which reconstructs the rational solution under the assumption that is has denominator dividing 360360, which suffices for most purposes. */ Rat := function(x) // Reconstructs rational from real approximation assuming denominator divides 360360 err := false; y := 360360*x; temp := Round(y); if (Abs(y-temp) gt 0.1) then err := true; end if; return temp/360360,err; end function; Red := procedure(~glist) // eliminates isomorphs and complementary pairs from graph list length := #glist; x := 1; while (x lt length) do H1 := glist[x]; H2 := Complement(H1); y := x+1; while (y le length) do H3 := glist[y]; if ( IsIsomorphic(H1,H3 : IgnoreLabels := true) or IsIsomorphic(H2,H3 : IgnoreLabels := true) ) then Remove(~glist,y); length -:= 1; else y +:= 1; end if; end while; x +:= 1; end while; end procedure; Fract := function(G); // finds value of maximum fractional packing of triangles in G and its complement V := VertexSet(G); E := EdgeSet(G); N := #V; EDG := N*(N-1) div 2; // number of edges NB. Slightly wasteful to compute each time! TR := N*(N-1)*(N-2) div 6; // number of triples trips := [ x: x in Subsets( {1..N} , 3) ]; // potential triangles (empty or complete) edges := [ x: x in Subsets( {1..N} , 2) ]; /* We run through all triples. When we find a complete or empty triangle we: 1. Add its `trips' index to a sequence `triangle' 2. For each of its edges we: i) Set the bit of `edge_is_real' corresponding to its index in `edges' to 1 ii) Add the number of the triangle to a sequence belonging to the edge. Further explanation of 2.ii: `triangle_on_edge' is a sequence of sequences: The first parameter is the index of an edge in `edges' The second parameter is the number of the triangle */ ntriang := 0; triangle := []; // list of triangles edge_is_real := [ 0: x in [1..EDG]]; // indicates edge in triangle triangle_on_edge := [ []: x in [1..EDG]]; // sequence of triangles on edge for x in { 1..TR } do current_trip := trips[x]; H := sub< G | {V!i : i in current_trip} >; if IsComplete(H) or IsEmpty(H) then ntriang +:= 1; Append(~triangle, x); // add triangle to list current_edges := Subsets(current_trip,2); for y in current_edges do current_index := Index(edges,y); edge_is_real[current_index] := 1; // mark y as a real edge Append( ~(triangle_on_edge[current_index]), ntriang ); // add triangle to y list end for; end if; end for; if (ntriang eq 0) then // only happens for N<=5 return 0; end if; // Now we count the edges that belong to triangles nreal_edges := 0; real_edge := []; for x in { 1..EDG } do if (edge_is_real[x] eq 1) then nreal_edges +:= 1; Append(~real_edge, x); end if; end for; // Now build the edge-triangle incidence matrix R := RealField(); lhs := Matrix(R, nreal_edges, ntriang, [ 0 : i in [1..nreal_edges*ntriang] ] ); for x in { 1..nreal_edges } do edge_index := real_edge[x]; for y in triangle_on_edge[edge_index] do ; lhs[ x, y ] := 1; end for; end for; // Solve the linear program rhs := Matrix(R, nreal_edges, 1, [ 1 : i in [1..nreal_edges] ] ); rel := Matrix(R, nreal_edges, 1, [ -1 : i in [1..nreal_edges] ] ); // negative values - leq relation obj := Matrix(R, 1, ntriang, [ 1 : i in [1..ntriang] ] ); vector := MaximalSolution( lhs, rel, rhs, obj )[1]; value := 0; for x in {1..ntriang} do value +:= vector[x]; end for; return value; end function; exhaust := function(N,k) /* Checks all graphs on N vertices looking for best k values. Output is the number of values given (up to max k), the values, and a list of graphs for up to k-1 values. */ best := []; graphs := []; nvalues := 0; F := OpenGraphFile("cmd " cat gengpath cat " " cat IntegerToString(N),0,0); while true do more, G := NextGraph(F); if (not more) then break; // reached end of F end if; realvalue := Fract(G); value, err := Rat(realvalue); // convert to rational if err then print "Rounding error: ",realvalue," to ",value; print G; end if; // Find best[index] <= value < best[index+1], with equality flag index := nvalues; equality := false; for j in {1..nvalues} do test := best[j]; if (value le test) then if (value eq test) then index := j; equality := true; break; else index := j-1; equality := false; break; end if; end if; end for; // update best and graphs if equality then Append(~graphs[index],G); else if (index lt k) then if (nvalues lt k) then nvalues +:= 1; end if; Insert(~best,index+1,value); Insert(~graphs,index+1,[G]); end if; end if; end while; returnlist := []; nlists := Minimum( {nvalues,k-1} ); for j in {1..nlists} do Red(~graphs[j]); returnlist cat:= graphs[j]; end for; return nvalues,best,returnlist; end function; fullextend := function(N,list,k,level) /* Checks one vertex extensions of list (a list of graphs on N-1 vertices) looking for the best values <= level, up to a maximum of k. Output is the number of values given, the values, and a list of graphs for all but the worst value. */ best := []; graphs := []; nvalues := 0; length := #list; for i in { 1..length } do G := list[i]; /* One vertex extensions of G are represented by N-1 bits, indicating adjacency to the Nth vertex. We check in binary order. */ bit := [0 : j in {1..N}]; // Nth bit indicates overflow AddVertex(~G,N); while (bit[N] eq 0) do realvalue := Fract(G); value, err := Rat(realvalue); // convert to rational if err then print "Rounding error: ",realvalue," to ",value; print G; end if; if (value le level) then // Find best[index] <= value < best[index+1], with equality flag index := nvalues; equality := false; for j in {1..nvalues} do test := best[j]; if (value le test) then if (value eq test) then index := j; equality := true; break; else index := j-1; equality := false; break; end if; end if; end for; // update best and graphs if equality then Append(~graphs[index],G); else if (index lt k) then if (nvalues lt k) then nvalues +:= 1; end if; Insert(~best,index+1,value); Insert(~graphs,index+1,[G]); end if; end if; end if; // move to next graph carry := true; x := 0; while carry do x +:= 1; bit[x] +:= 1; if (bit[x] eq 1) then carry := false; AddEdge(~G,x,N); else bit[x] := 0; RemoveEdge(~G,x,N); end if; end while; end while; end for; returnlist := []; nlists := Minimum( {nvalues,k-1} ); for j in {1..nlists} do Red(~graphs[j]); returnlist cat:= graphs[j]; end for; return nvalues,best,returnlist; end function; // start of program readi d,"Search depth: "; readi N,"Value to exhaust from: "; readi Nmax,"Stop at: "; t := Cputime(); nvals,best,graphs:= exhaust(N,d); print N," vertices: ",[best[j] : j in [1..nvals]]," ",Cputime(t)," seconds."; level := best[nvals]; N +:= 1; while (N le Nmax) do if (nvals eq d) then level := best[nvals]*N/(N-2); else level := level*N/(N-2); end if; print "level: ",level; // minimum value if not extension of member of graphs t := Cputime(); nvals,best,graphs := fullextend(N,graphs,d,level); if (nvals eq 0) then print "No more values!"; break; end if; print Cputime(t)," seconds."; print N," vertices: ",[best[j] : j in [1..nvals]]; N +:= 1; end while;