// commentaar in Magma staat achter twee /'s

// Het is handig om voor een aantal eenvoudige berekingen kleine
// functies aan te maken.
// Hieronder zijn een aantal voorbeelden met commentaar verzameld.

// berekenen van n!
// toepassing: factorial(10) geeft het resultaat 10! = 3628800
// syntax voor een for-loop is:
// for var in lijst do
//   een of meerdere aanwijzingen (eindigend op ;)
// end for;
// het return-statement geeft de waarde aan die teruggeven wordt
factorial := function(n)
  fac := 1;
  for i in [2..n] do
    fac := i * fac;
  end for;
  return fac;
end function;

// recursie: binnen een functie mag de functie zelf gebruikt worden
// de naam wordt hierbij door $$ vervangen
// syntax van if-statement:
// if voorwaarde then
//   een of meerdere aanwijzingen
// elif aandere voorwaarde then
//   een of meerdere aanwijzingen
// else 
//   een of meerdere aanwijzingen
// end if; 
rec_factorial := function(n)
  if n eq 0 then
    return 1;
  else
    return n * $$(n-1);
  end if;
end function;

// berekenen van de binomiaalcoefficient volgens de definitie n!/(k! (n-k)!)
// merk op: a div b geeft het gehele deel van a/b terug,
// terwijl a/b een rationaal getal terug geeft (ook in het geval dat
// b een deler van a is)
binomial := function(n, k)
  if k lt 0 or k gt n then
    return 0;
  elif k eq 0 or k eq n then
    return 1;
  else
    return factorial(n) div (factorial(k) * factorial(n-k));
  end if;
end function;

// recursie voor de binomiaalcoefficienten volgens
// (n // k) = (n-1 // k-1) + (n-1 // k)
// let op bij het toepassen van deze functie, loopt snel vast
// met CTRL-C laat zich Magma stoppen
// met twee keer CTRL-C binnen een seconde maak je een emergency exit  
rec_binomial1 := function(n, k)
  if k lt 0 or k gt n then
    return 0;
  elif k eq 0 or k eq n then
    return 1;
  else
    return $$(n-1, k-1) + $$(n-1, k);
  end if;
end function;

// recursie voor de binomiaalcoefficienten volgens
// (n // k) = (n // k-1) * (n-k+1)/k
rec_binomial2 := function(n, k)
  if k lt 0 or k gt n then
    return 0;
  elif k eq 0 or k eq n then
    return 1;
  else
    return $$(n, k-1) * (n-k+1)/k;
  end if;
end function;

// functie om in een loop over alle deelverzamelingen van een verzameling
// te lopen (algoritme 3.12.2 in Cameron)
// deze functie vormt de op Y volgende deelverzameling
// merk op: in Magma bestaan ook verzamelingen, maar hier worden de
// verzamelingen als (rijen) sequences gemaakt, waarbij de elementen 
// in een rij echter altijd opstijgend zijn
// algemeen is een verzameling gegeven door {1, 9, 4}, waarbij
// {1, 9, 4} = {9, 4, 1} = {1, 1, 9, 4, 9}
// terwijl rijen geordend zijn, dus [1, 9, 4] en [1, 4, 9] zijn verschillend
// syntax voor rij: 
// [ x : iets waarvan x afhangt | voorwaarde waar x aan voldoet ]
// #rij is het aantal elementen in de rij, dus #[1, 9, 4] = 3
next_subset := function(n, Y)
  i := n;
  while i ge 1 and i in Y do
    i := i-1;
  end while;
  if i eq 0 then
    return [1..n+1];
  else
// Z wordt de rij van de elementen van de rij Y die kleiner zijn dan i
    Z := [ Y[j] : j in [1..#Y] | Y[j] lt i ];
    Z := Z cat [i]; // aanvoegen van i aan het eind van Z
  end if;
  return Z;
end function;

// procedure om in een loop over alle deelverzamelingen van een verzameling
// te lopen (algoritme 3.12.2 in Cameron)
// deze functie is de hoofdloop
// een procedure heeft (in tegenstelling tot een functie) geen return-statement
subsets := procedure(n)
  Y := [];
  count := 1;
  while #Y le n do
    print count, Y;
    Y := next_subset(n, Y);
    count := count+1;
  end while;
end procedure;

// eenvoudige manier om de deelverzamelingen van {1..n} aan te maken,
// gebaseerd op de bijectie tussen de deelverzamelingen en de
// 2-tallige ontwikkelingen van de getallen 0, 1, ..., 2^n-1
binary_subsets := procedure(n)
  for i in [0..2^n-1] do
// s := IntegerToSequence(n, m) maakt de m-tallige ontwikkeling van n (als rij)
// met coefficienten voor de opstijgende machten van m, dus
// n = s[1] m^0 + s[2]*m^1 + ... + s[k]*m^(k-1)
    s := IntegerToSequence(i, 2);
    Y := { j : j in [1..#s] | s[j] eq 1 };
    print i+1, Y;
  end for;
end procedure;

// functie om in een loop over alle k-elementen deelverzamelingen van een 
// verzameling met n elementen te lopen (algoritme 3.12.3 in Cameron)
// deze functie vormt de op Y volgende deelverzameling
next_k_subset := function(n, k, Y)
  i := 1;
  while i lt k and Y[i+1] eq Y[i]+1 do
    i := i+1;
  end while;
  if i eq k and Y[i] eq n then
    return [];
  else
    Z := [ j : j in [1..k] ];
    Z[i] := Y[i]+1;
    for j in [i+1..k] do
      Z[j] := Y[j];
    end for;
    return Z;
  end if;
end function;

// procedure om in een loop over alle k-elementen deelverzamelingen van een 
// verzameling met n elementen te lopen (algoritme 3.12.3 in Cameron)
// deze functie is de hoofdloop
k_subsets := procedure(n, k)
  Y := [ i : i in [1..k] ];
  count := 1;
  while #Y gt 0 do
    print count, Y;
    Y := next_k_subset(n, k, Y);
    count := count+1;
  end while;
end procedure;

// recursieve functie om de partities van {1..n} aan te maken
// (algoritme 3.13.3 in Cameron)
// Let op: probeer dit niet voor X met meer dan 10 elementen!
// Exclude haalt het element x uit de rij X weg (als het er in zit)
// Seqset maakt van een rij (sequence) een verzameling (set)
// Setseq maakt van een verzameling weer een rij
// Append(~part, p) hangt de nieuwe partitie p achter de rij part
// zonder een kopie van part te maken
partitions := function(X)
  if #X eq 0 then
    return [ {} ];
  else
    part := [];
    x := X[1];
    XX := Exclude(X, x);
    for i in [0..2^#XX-1] do
      s := IntegerToSequence(i, 2);
      Y := [ XX[j] : j in [1..#s] | s[j] eq 1 ];
      XXX := Setseq(Seqset(XX) diff Seqset(Y));
      partX := $$(XXX);
      for j in [1..#partX] do
        Append(~part, partX[j] join {[x] cat Y});
      end for;
    end for;
    return part;
  end if;
end function;

