// Determines the steiner systems via interpolation of quadrics through
// contact points.
function ComputeSteinerSetsBySyzygy(C, Tgt_list)

  P3 := AmbientSpace(C);
  L := LinearSystem(P3, 2);
  linearSystemTable := AssociativeArray();

  // Construct a table of linear systems
  for ell in Tgt_list do
	  // WARNING: code fails if there is a flex line.
	  contact_points := ReducedSubscheme(C meet Scheme(P3, ell));
	  assert Degree(contact_points) eq 3;
	  L2 := LinearSystem( L, contact_points);
	  linearSystemTable[ell] := L2;
  end for;

  function IsSyzygetic(pr1, pr2)
	  if not IsEmpty(pr1 meet pr2) then
	    return false;
	  end if;

	  L3 := &meet [linearSystemTable[ell] : ell in pr1 join pr2];
	  if Dimension(L3) eq 1 then
	    return true;
	  elif Dimension(L3) eq 0 then
	    return false;
	  else
	    error "Something bad happened", Dimension(L2);
	  end if;
  end function;

  pairs := Subsets( Seqset(Tgt_list), 2);
  SS := [* {} : i in [1 .. 255] *];

  // Partition the pairs into Steiner sets
  while not IsEmpty(pairs) do
	  print #pairs;
	  p1 := Random(pairs);
	  i := 1;

	  while i le 255 do
	    if IsEmpty(SS[i]) then
		    SS[i] := {p1};
		    Exclude(~pairs, p1);
		    break;
	    elif IsSyzygetic(p1, Random(SS[i])) then
		    Include(~SS[i], p1);
		    Exclude(~pairs, p1);
		    break;
	    end if;
	    i +:= 1;
	  end while;

	  error if i gt 255, "Steiner set not found for pair", [#ss : ss in SS];
  end while;

  assert &and [#ss eq 28 : ss in SS];
  return { ss : ss in SS};
end function;


// Determines the steiner systems via Riemann Roch spaces.
function ComputeSteinerSetsByLinearEquivalence(C, Tgt_list)

  P3 := AmbientSpace(C);
  L := LinearSystem(P3, 2);
  divisorTable := AssociativeArray();

  // Construct a table of linear systems
  for ell in Tgt_list do
	  D := Divisor(C, C meet Scheme(P3, ell)) div 2;
	  assert Degree(D) eq 3;
	  divisorTable[ell] := D;
  end for;

  function IsSyzygetic(pr1, pr2)
	  if not IsEmpty(pr1 meet pr2) then
	    return false;
	  end if;

	  cmp1:= [divisorTable[ell] : ell in Setseq(pr1)];
	  D1  := cmp1[1] - cmp1[2];
	  cmp2:= [divisorTable[ell] : ell in Setseq(pr2)];
	  D2  := cmp2[1] - cmp2[2];

	  return IsLinearlyEquivalent(D1,D2);
  end function;

  pairs := Subsets( Seqset(Tgt_list), 2);
  SS := [* {} : i in [1 .. 255] *];

  // Partition the pairs into Steiner sets
  while not IsEmpty(pairs) do
	  print #pairs;
	  p1 := Random(pairs);
	  i := 1;

	  while i le 255 do
	    if IsEmpty(SS[i]) then
		    SS[i] := {p1};
		    Exclude(~pairs, p1);
		    break;
	    elif IsSyzygetic(p1, Random(SS[i])) then
		    Include(~SS[i], p1);
		    Exclude(~pairs, p1);
		    break;
	    end if;
	    i +:= 1;
	  end while;

	  error if i gt 255, "Steiner set not found for pair", [#ss : ss in SS];
  end while;

  assert &and [#ss eq 28 : ss in SS];
  return { ss : ss in SS};
end function;
