Add EnumeratorOfPartitionsSet enumerator

[reiniscirpons]

This PR adds an enumerator for the PartitionsSet object, which returns all set partitions of a finite set, with the possibility of restricting the number of parts too. This PR resolves #6362.

Interestingly, the enumerator is faster than the iterator for a fixed k (when k is not fixed, they seem to run comparably fast):

gap> A := EnumeratorOfPartitionsSet([1..20], 10);;
gap> A[1000]; time;
[ [ 1, 2, 5, 6, 12 ], [ 3, 4, 7, 8, 9, 10, 11 ], [ 13 ], [ 14 ], [ 15 ], [ 16 ], [ 17 ], [ 18 ], [ 19 ], [ 20 ] ]
1
gap> B := IteratorOfPartitionsSet([1..20], 10);;
gap> NextIterator(B); time;
[ [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ], [ 12 ], [ 13 ], [ 14 ], [ 15 ], [ 16 ], [ 17 ], [ 18 ], [ 19 ], [ 20 ] ]
584

I also used the Rocq theorem prover at the initial stages of writing this function to help get the encoding/decoding logic correct (see Further details section for more details).

Text for release notes

see title

Further details

The enumerator is implemented using EnumeratorByFunctions. To do this we needed to implement two functions: an encoder, which takes a partition and returns its number with respect to some fixed order (NumberElement), and a decoder which does the opposite (ElementNumber). It was not initially clear to me how to do this, I knew I would make use of the identity Stirling2(n, k) = Stirling2(n-1, k-1) + k * Stirling2(n-1, k), but how exactly to do this would take a lot of iteration on my part if I were to do it directly in GAP. Another option would be to sketch it out on paper, but I decided to instead model it in Rocq. I am including below an annotated proof script I came up with as a result of this, for those interested. It should work with Rocq 9 and mathcomp >= 2.5:

Require Import ssreflect ssrbool ssrfun ssrbool.
From mathcomp Require Import ssrnat div seq zify.

Section PartitionSet.

(* A function for recursively computing Stirling numbers
   of the second kind. *)
Fixpoint Stirling2 (n k: nat): nat :=
  match n, k with
  | 0, 0 => 1
  | n'.+1, k'.+1 => Stirling2 n' k' + k * Stirling2 n' k
  | _, _ => 0
  end.

(* We model partitions as an inductively defined type.
   Instead of using lists of lists of elements as is used
   in GAP I decided to model them as a sequence of operations
   for that incrementally build up the partition starting with
   the empty partition as follows. *)
(* partition n k is the type of partitions of a set of n elements
   into k parts. *)
Inductive partition: nat -> nat -> Type :=
  (* The empty partition serves as our base case - it is a partition
     with 0 elements and 0 parts. *)
| EmptyPartition: partition 0 0
  (* We can construct a partition with n+1 elements and k+1 parts by
     taking an existing partition p with n elements and k parts
     and adding a new singleton part to p using the NewPart constructor. *)
| NewPart (n k: nat) (p: partition n k):
    partition (n.+1) (k.+1)
  (* We can construct a partition with n+1 elements and k parts by
     taking an existing partition p with n elements and k parts
     and adding a new element to the i-th part with the AddToPart
     constructor. *)
| AddToPart (n k i: nat) (p: partition n k):
    partition (n.+1) k.

(* The following function checks if a partition is well defined, namely
   we check that we only ever add an integer to an existing part. *)
Fixpoint partition_well_defined {n k: nat} (p: partition n k): bool :=
  match p with
  | EmptyPartition => true
  | NewPart _ _ p => partition_well_defined p
  | AddToPart _ k' i p => (i < k') && partition_well_defined p
  end.

(* The pickle function encodes a partition as an integer. *)
Fixpoint partition_pickle (n k: nat) (p: partition n k): nat :=
  match p with
  | EmptyPartition => 0
  | NewPart n' k' p' => partition_pickle n' k' p'
  | AddToPart n' k' i p' =>
    partition_pickle n' k' p' + i * Stirling2 n' k' + Stirling2 n' k'.-1 
  end.

(* We can formulate certain properties we want to hold for the encoder
   function and stipulate them as lemmas. For example, the lemma
   partition_pickle_ltn stipulates that when encoding a partition
   with n elements and k parts does not exceed Stirling2 n k.
   
   My initial implementation of the function was wrong, and I was
   unable to prove this lemma, which prompted me to tweak the definition,
   until eventually it was right. *)
Lemma partition_pickle_ltn: forall n k (p: partition n k),
  partition_well_defined p -> partition_pickle n k p < Stirling2 n k.
Proof.
  move => n k.
  elim => [//|n' k' p IH /=|n' [//|k'] i p IH /=].
  - by lia.
  - by nia.
Qed.

(* The unpickle function does the decoding. We allow it to fail on invalid
   inputs by returning None. *)
Let helper_i (n k m: nat): nat := (m - Stirling2 n k) %/ (Stirling2 n k.+1).
Let helper_m (n k m: nat): nat := (m - Stirling2 n k) %% (Stirling2 n k.+1).
Fixpoint partition_unpickle (n k m: nat): option (partition n k) :=
  if m < Stirling2 n k then
    match n, k, m with
    | 0, 0, 0 => Some EmptyPartition
    | n'.+1, k'.+1, m => 
      if m < Stirling2 n' k' then
        if partition_unpickle n' k' m is Some p then
          Some (NewPart n' k' p)
        else
          None
      else
        if partition_unpickle n' k'.+1 (helper_m n' k' m) is Some p then
          Some (AddToPart n' k'.+1 (helper_i n' k' m) p)
        else
          None
    | _, _, _ => None
    end
  else
    None.

(* Here we require that the decoding function is defined if and only if
   the input number is less than Stirling2 n k. *)
Lemma partition_unpickleP: forall n k m,
  reflect (exists p, partition_unpickle n k m = Some p) (m < Stirling2 n k).
Proof.
  move => n k m; apply (iffP idP).
  elim: n k m => [|n IH] [|k m /= H] //;
    first by case => [_|//] /=; by exists EmptyPartition.
  - rewrite H /=; case: ifP => Hm.
  -- move: (IH _ _ Hm) => [p' ->]. by exists (NewPart n k p'). 
  -- have: (helper_m n k m < Stirling2 n k.+1) =>
       [|/IH [p' ->]]; first by rewrite /helper_m; lia.
     by exists (AddToPart n k.+1 (helper_i n k m) p').
  - case: n => [|n]; case: k => [|k] [p] //=;
      first by case: m.
    by case: ifP.
Qed.

(* We also prove that the result of the decoding function is a well defined
   partition. *)
Lemma partition_unpickle_well_defined: forall n k m p,
  partition_unpickle n k m = Some p -> partition_well_defined p.
Proof.
  elim => [|n IH]; case => [|k] //;
    first by case => [p /= [] <-|].
  move => m p /=.
  case: ifP => [Hm |//]; case: ifP => Hm';
  case H: (partition_unpickle _ _ _) => [p'|//]; case => <- /=.
  - by apply IH with m.
  - apply/andP; split;
      last by apply (IH _ _ _ H).
    by rewrite /helper_i; nia.
Qed.

(* Finally, we prove that the pickle and unpickle functions are mutually
   inverse when the partition is well defined. *)
Lemma partition_pickleK: forall n k p,
  partition_well_defined p ->
  partition_unpickle n k (partition_pickle n k p) = Some p.
Proof. 
  move => n k; elim => [//|n' k' p' /= H Hp|n' k' i p' /= H /andP [Hi Hp]];
  move: (H Hp) => {}H.
  - rewrite H;
    have: (partition_pickle n' k' p' < Stirling2 n' k') => [|H'];
      first by apply/partition_unpickleP; exists p'.
    rewrite H'; case: ifP => [//|]; by lia.
  - case: k' p' H Hp Hi => [//|k' p' H Hp Hi].
    have: (partition_pickle n' k'.+1 p' < Stirling2 n' k'.+1) => [|H'];
    first by apply/partition_unpickleP; exists p'.
    case: ifP => [_|] /=; last by nia.
    case: ifP => [|_]; first by lia.
    set m' := helper_m n' k' _.
    set i' := helper_i n' k' _.
    have: (m' = partition_pickle n' k'.+1 p') => [|->];
      first by rewrite /m' /helper_m -addnBA // subnn addn0 addnC modnMDl modn_small.
    have: (i' = i) => [|->];
      first by rewrite /i' /helper_i; nia.
    by rewrite H.
Qed.

End PartitionSet.

It was useful since I could detect when the encode/decode functions are wrong: if a function was wrong, I would get stuck when proving a lemma. But the way I got stuck in the proof informed me on what I needed to change in the function to make it more correct. It took about 2-3 iterations until I got the functions right and got all of the proofs correct.

This gave me a good idea for what the GAP functions should look like. I could not translate the functions 1-1 to GAP, since the Rocq script uses a functional style. If I had more Rocq-fu I could probably write it closer to an imperative style. So the main complexity was unrolling the recursion when implementing these functions in GAP and translating the inductive partition definition to the list definition. But this was not too bad, and the big advantage was that I did not have to second guess the mathematical foundations of my encode/decode function.

[codecov]

Codecov Report

❌ Patch coverage is 93.96552% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 78.70%. Comparing base (aa52ee2) to head (cc23354).
⚠️ Report is 8 commits behind head on master.

Files with missing lines	Patch %	Lines
lib/combinat.gi	93.91%	7 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff            @@
##           master    #6384    +/-   ##
========================================
  Coverage   78.69%   78.70%            
========================================
  Files         684      684            
  Lines      292892   293008   +116     
  Branches     8660     8686    +26     
========================================
+ Hits       230497   230607   +110     
- Misses      60581    60587     +6     
  Partials     1814     1814            

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[fingolfin]

Any reason not do just do

Suggested change

	## gap> i := 0;; for c in PartitionsSet(m) do i := i+1; od;
	## gap> i;
	## gap> Length(PartitionsSet(m));

[fingolfin]

How about this instead:

Suggested change

	## gap> Position(cm, [[1, 3], [2, 5], [4], [6, 7, 8, 9]]);
	## 11001
	## gap> cm[11001];
	## [ [ 1, 3 ], [ 2, 5 ], [ 4 ], [ 6, 7, 8, 9 ] ]
	## gap> Position(cm, last);
	## 1000

[fingolfin]

Very nice achievement, @reiniscirpons, thanks!

[fingolfin]

The documentation seems fine, the feature is useful. I didn't try to understand the code in depth, that'd require a lot of work. But I trust @reiniscirpons who put in a lot of work, the code is derived from a formalization and well-documented, and there are tests.

Maybe two very minor tweaks can be considered for the documentation, but they are not required.

[limakzi]

I am unsure with this comment.
The iterator accepts <positive integer> and <set> as first argument.

PartitionsSet accepts only <set>.
GAP proposes other function for <positive integer> and its called Partitions.
I would profer to have primary function, focusing entirely each of these.

Same goes here, I would prefer to accept positive integer here too.