Constructing Steiner Triple Systems with OCaml

A steiner system, S(a,b,c), where a,b,c are all positive integers with a <= b <= c, is defined as a set of c elements, a collection of subsets of this set, each of size b, such that every a elements exist in exactly one subset.

e.g. S(2,3,7)

We need a set of 7 elements,

S = {0,1,2,3,4,5,6},

and we need to create subsets of S, each of size 3, where each pair of elements exist in exactly one subset:

b1 = {0,1,2}

b2 = {0,3,4}

b3 = {0,5,6}

b4 = {1,4,5}

b5 = {1,3,6}

b6 = {2,4,6}

b7 = {2,3,5}

From observation, it is clear that any two elements exist in exactly one subset b1,...,b7. We shall refer to each subset as a block.

Clearly, a steiner system does not necessarily exist for any given a,b,c. For example, there is no S(2,7,8) system. You cannot generate subsets of {0,1,2,3,4,5,6,7}, each of size 7, such that any two elements are in exactly one subset.

An S(a,b,c) is called a steiner triple system if a = 2, b = 3, and can be denoted STS(n). It can be easily proven that a S(2,3,n) exists for any n such that

$$n\cong1,3\mod6$$
But, this fact doesn't actually give us the triples of the STS(n), which is what we are interested in. The above example's 7 blocks were found by trial and error, and this isn't really a scalable method for constructing triple systems. It is not immediately clear how we can construct an STS(n) for a given n. This post creates the algorithm(s) necessary to create a STS(n) for an any applicable n.

The algorithm
We first split the problem into two parts. We will solve the problem for the case of $n\cong 1\mod 6$ and $n\cong3\mod6$ separately.

Case 1: $$n\cong3\mod6$$ . This is sometimes known as the Bose Construction. I won't go into details but will just present the OCaml source code, as it is not particularly long or complicated:

 let construct_for_residue_3 v =  
     let res = (v - 3) / 6 in  
     let arr = ref [] in  
     let sz = (2*res) + 1 in  
     let mult = (sz + 1) / 2 in  
     for i = 0 to sz - 1 do  
         arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
     done;  
     for j = 0 to sz - 2 do  
         for k = (j+1) to sz - 1 do  
             arr := (List.append (!arr) [[[j;0];[k;0];[(mult*(j + k)) mod sz; 1]];  
                             [[j;1];[k;1];[(mult*(j+k)) mod sz; 2]];  
                             [[j;2];[k;2];[(mult*(j+k)) mod sz; 0]]])  
         done  
     done;  
     !arr  

As an example,  let's try it for the case where n = 9.

 let l = construct_for_residue_3 9;;  

which gives

 [[[0; 0]; [0; 1]; [0; 2]]; [[1; 0]; [1; 1]; [1; 2]]; [[2; 0]; [2; 1]; [2; 2]];                         [[0; 0]; [1; 0]; [2; 1]]; [[0; 1]; [1; 1]; [2; 2]]; [[0; 2]; [1; 2]; [2; 0]];                         
   [[0; 0]; [2; 0]; [1; 1]]; [[0; 1]; [2; 1]; [1; 2]]; [[0; 2]; [2; 2]; [1; 0]];   
   [[1; 0]; [2; 0]; [0; 1]]; [[1; 1]; [2; 1]; [0; 2]]; [[1; 2]; [2; 2]; [0; 0]]]  

We have a List of 12 elements, where each element contains 3 sublists of ordered number pairs. e.g. the first element is

[[0; 0]; [0; 1]; [0; 2]];

Note that any two number pairs [i; j], [k;l] (where i,k is in range [0,2] and j,l is in range [0,2]),  exist in exactly one sublist.

Case 2:$n\cong1\mod6$
This case is slightly more complicated, but the approach is the same. The source code is given below:

 let construct_for_residue_1 v =  
     let res = (v - 1) / 6 in  
     let arr = ref [] in  
     for i = 0 to res - 1 do  
         arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
     done;  
     for j = 0 to 2 * res - 2 do  
         for k = j + 1 to 2 * res - 1 do  
             let z = qgf j k (2 * res) in  
             arr := (List.append (!arr) [[[j;0];[k;0];[z;1]];  
                             [[j;1];[k;1];[z;2]];  
                             [[j;2];[k;2];[z;0]]])  
         done;  
     done;  
     for l = 0 to res - 1 do  
         arr := (List.append (!arr) [[[(-1);(-1)];[l+res;0];[l;1]];  
                         [[(-1);(-1)];[l+res;1];[l;2]];  
                         [[(-1);(-1)];[l+res;2];[l;0]]])  
     done;  
     !arr  

Here, the [i,j] pairs range over different integer values but the triples are the same, i.e. a triple contains three [i,j] pairs. And every two pairs exist in exactly one triple.

For example lets try it with n = 7.

 let l = construct_for_residue_1 7;;  

This gives

  [[[0; 0]; [0; 1]; [0; 2]]; [[0; 0]; [1; 0]; [1; 1]]; [[0; 1]; [1; 1]; [1; 2]];

  [[0; 2]; [1; 2]; [1; 0]]; [[-1; -1]; [1; 0]; [0; 1]]; [[-1; -1]; [1; 1]; [0; 2]];                       
   [[-1; -1]; [1; 2]; [0; 0]]]  

Here, we can see 7 triples and it is easy to check any [[i;j];[k;l]] pair exists in one triple. We can map pairs to the b's in the example at the top of the page. For example

[0;0] ->   b1
[0;1] ->   b2
[0;2] ->   b3
[1;0] ->   b4
[1;1] ->   b5
[1;2] ->   b6
[-1;-1] -> b7

A good question at this point might be, "What's the point in this?", to which the answer is that steiner systems and combinatorial designs in general, have a wide range of connections to different areas of mathematics -  for example:

Projective and Affine Geometry (in fact the above STS(7) steiner system can be considered a projective plane);

Coding theory (e.g. Golay codes);

Group theory (Not going to go into it here, but the above STS(7) has 168 symmetries, and its group of
symmetries is isomorphic to PGL(3,2));

Experimental designs.

It's also worth noting that there are many more interesting steiner systems, beyond the triple system examples here. Probably the most noteworthy is S(5,8,24). Just a quick recap on what that means: 24 points, grouped into subsets, each of size 8, such that any given 5 points exist in exactly 1 subset. That's the definition, but the interesting point is that this object is immensely symmetric, and has over 140 million symmetries, and what's more, any 5 points can be mapped on to any other five points and still preserve symmetry.

source:

1:  (* steiner triple systems *)  
2:  open Core.Std  
3:  exception Impossible_construction  
4:    
5:    
6:  (* construct STS for vertex size 3 modulo 6 *)  
7:  let construct_for_residue_3 v =  
8:      let res = (v - 3) / 6 in  
9:      let arr = ref [] in  
10:      for i = 1 to (2 * res + 1) do  
11:          arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
12:      done;  
13:        
14:      for j = 1 to (2 * res + 1) do  
15:          for k = (j+1) to (2 * res + 1) do  
16:              arr := (List.append (!arr) [[[j;1];[k;1];[(j + k) / 2; 1]];  
17:                              [[j;2];[k;2];[(j+k) / 2; 2]];  
18:                              [[j;0];[k;0];[(j+k) / 2;0]]])  
19:          done  
20:      done;  
21:      !arr  
22:    
23:  (* quasigroup function -- just a function used for Skolem's construction *)  
24:  let qgf left right max =  
25:      let z = (left + right) mod max in  
26:      if z mod 2 = 0 then  
27:          z / 2  
28:      else  
29:          (z + max - 1) / 2  
30:    
31:  (* Skolem's construction - for vertex size 1 modulo 6 *)  
32:  let construct_for_residue_1 v =  
33:      let res = (v - 1) / 6 in  
34:      let arr = ref [] in  
35:      for i = 0 to res - 1 do  
36:          arr := (List.append (!arr) [[[i;0];[i;1];[i;2]]])  
37:      done;  
38:    
39:      for j = 0 to 2 * res - 2 do  
40:          for k = j + 1 to 2 * res - 1 do  
41:              let z = qgf j k (2 * res) in  
42:              arr := (List.append (!arr) [[[j;0];[k;0];[z;1]];  
43:                              [[j;1];[k;1];[z;2]];  
44:                              [[j;2];[k;2];[z;0]]])  
45:          done;  
46:      done;  
47:    
48:      for l = 0 to res - 1 do  
49:          arr := (List.append (!arr) [[[(-1);(-1)];[l+res;0];[l;1]];  
50:                          [[(-1);(-1)];[l+res;1];[l;2]];  
51:                          [[(-1);(-1)];[l+res;2];[l;0]]])  
52:      done;  
53:    
54:      !arr  
55:    
56:  (* find all blocks intersecting the element *)  
57:  let filter_element triplesystem e =  
58:      List.filter triplesystem (fun i -> List.mem i e)