// Every n-card layout up to the deck's symmetries (AGL(4,3), 1,965,150,720
// maps), by McKay's canonical augmentation with nauty. Independent of the
// counting argument: it reads no maxima, it just lists one layout per class
// and counts its sets.
//
// A layout is a colouring of the deck's point-line incidence graph (81 points,
// 1080 lines): chosen points, other points, lines. That graph's automorphisms
// are exactly AGL(4,3), so nauty's canonical labelling decides equivalence.
// One level: read the classes of (k-1)-card layouts, extend each by one card
// from each orbit of its stabilizer, and keep a child only if the added card is
// in the orbit of the child's canonical card to delete. That keeps exactly one
// layout per class (McKay 1998). The card to delete is chosen among the cards
// with the largest invariant (sets through the card, then the sum over the 40
// hyperplanes through it of the square of the layout's cards in each), first
// in nauty's canonical order. A child whose added card doesn't have the
// largest invariant is rejected before nauty runs.
//
// nauty isn't in the repo: build it in a temp dir, then
//   cc -O2 -I$NAUTY -o orbits orbits.c $NAUTY/nauty.a
//   ./orbits <in-file> <out-file> <part> <parts> [max-sets]
// With max-sets, only layouts holding at most that many sets are kept (a
// property every sub-layout shares, so the augmentation stays complete):
// max-sets 0 lists the caps, the layouts with no set at all.
// in-file holds one layout per line as 81 '0'/'1' characters ("-" = start from
// nothing). Output: the accepted layouts, and on stderr a summary line:
//   k classes maxsets classes-at-max  then the histogram "sets:classes ..."
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "nausparse.h"

#define P 81
#define NL 1080
#define NV (P + NL)

static int third[P][P];
static int line_pts[NL][3];
static int pt_lines[P][40];
static int pt_hyp[P][40];  // the 40 hyperplanes through each card
static char hyp[120][P];
static sparsegraph g;

static void build(void) {
  for (int a = 0; a < P; a++)
    for (int b = 0; b < P; b++) {
      int c = 0, pa = a, pb = b, w = 1;
      for (int i = 0; i < 4; i++, pa /= 3, pb /= 3, w *= 3) c += ((6 - pa % 3 - pb % 3) % 3) * w;
      third[a][b] = c;
    }
  int nl = 0, deg[P] = {0};
  for (int a = 0; a < P; a++)
    for (int b = a + 1; b < P; b++) {
      int c = third[a][b];
      if (c > b) {
        line_pts[nl][0] = a; line_pts[nl][1] = b; line_pts[nl][2] = c;
        pt_lines[a][deg[a]++] = nl; pt_lines[b][deg[b]++] = nl; pt_lines[c][deg[c]++] = nl;
        nl++;
      }
    }
  if (nl != NL) { fprintf(stderr, "%d lines\n", nl); exit(1); }
  int nh = 0, hd[P] = {0};
  for (int w = 1; w < P; w++) {
    int d[4] = {w % 3, w / 3 % 3, w / 9 % 3, w / 27}, lead = 0;
    while (!d[lead]) lead++;
    if (d[lead] != 1) continue;
    for (int t = 0; t < 3; t++, nh++)
      for (int x = 0; x < P; x++)
        if ((hyp[nh][x] = (d[0] * (x % 3) + d[1] * (x / 3 % 3) + d[2] * (x / 9 % 3) + d[3] * (x / 27)) % 3 == t))
          pt_hyp[x][hd[x]++] = nh;
  }
  SG_INIT(g);
  SG_ALLOC(g, NV, 2 * 3 * NL, "build");
  g.nv = NV; g.nde = 2 * 3 * NL;
  size_t pos = 0;
  for (int v = 0; v < NV; v++) {
    g.v[v] = pos;
    if (v < P) { g.d[v] = 40; for (int i = 0; i < 40; i++) g.e[pos++] = P + pt_lines[v][i]; }
    else { g.d[v] = 3; for (int i = 0; i < 3; i++) g.e[pos++] = line_pts[v - P][i]; }
  }
}

// canonical labelling of the layout `in`; fills lab and orbits
static void canon(const char* in, int k, int* lab, int* orbits) {
  static DEFAULTOPTIONS_SPARSEGRAPH(opt);
  statsblk st;
  int ptn[NV];
  opt.getcanon = TRUE;
  opt.defaultptn = FALSE;
  int i = 0;
  for (int v = 0; v < P; v++) if (in[v]) lab[i++] = v;
  for (int v = 0; v < P; v++) if (!in[v]) lab[i++] = v;
  for (int v = P; v < NV; v++) lab[i++] = v;
  for (int j = 0; j < NV; j++) ptn[j] = 1;
  if (k > 0) ptn[k - 1] = 0;
  ptn[P - 1] = 0; ptn[NV - 1] = 0;
  SG_DECL(cg);
  sparsenauty(&g, lab, ptn, orbits, &opt, &st, &cg);
  SG_FREE(cg);
}

// the invariant of card v in layout `in` (see the top)
static long invariant(const char* in, const int* cnt, int v) {
  int deg = 0;
  long sq = 0;
  for (int i = 0; i < 40; i++) {
    const int* l = line_pts[pt_lines[v][i]];
    deg += in[l[0]] && in[l[1]] && in[l[2]];
    sq += (long)cnt[pt_hyp[v][i]] * cnt[pt_hyp[v][i]];
  }
  return deg * 1000000L + sq;
}

static int count_sets(const char* in) {
  int c = 0;
  for (int l = 0; l < NL; l++) c += in[line_pts[l][0]] && in[line_pts[l][1]] && in[line_pts[l][2]];
  return c;
}

int main(int argc, char** argv) {
  if (argc < 5) { fprintf(stderr, "usage: orbits in out part parts\n"); return 1; }
  int part = atoi(argv[3]), parts = atoi(argv[4]);
  int maxsets = argc > 5 ? atoi(argv[5]) : NL;
  build();
  FILE* fi = strcmp(argv[1], "-") ? fopen(argv[1], "r") : NULL;
  FILE* fo = fopen(argv[2], "w");
  char buf[256], in[P];
  int lab[NV], orbits[NV], clab[NV], corb[NV];
  long idx = 0, classes = 0, hist[NL + 1] = {0};
  int k = -1;
  for (;;) {
    if (fi) { if (!fgets(buf, sizeof buf, fi)) break; }
    else if (idx > 0) break;
    if (idx++ % parts != part) continue;
    int kk = 0;
    for (int v = 0; v < P; v++) { in[v] = fi ? buf[v] == '1' : 0; kk += in[v]; }
    k = kk + 1;
    canon(in, kk, lab, orbits);  // orbits of the parent's stabilizer
    for (int x = 0; x < P; x++) {
      if (in[x] || orbits[x] != x) continue;  // one card per orbit outside the layout
      in[x] = 1;
      if (count_sets(in) > maxsets) { in[x] = 0; continue; }
      int cnt[120];
      for (int h = 0; h < 120; h++) { cnt[h] = 0; for (int v = 0; v < P; v++) cnt[h] += in[v] && hyp[h][v]; }
      long inv[P], best = -1;
      for (int v = 0; v < P; v++) if (in[v]) { inv[v] = invariant(in, cnt, v); if (inv[v] > best) best = inv[v]; }
      if (inv[x] < best) { in[x] = 0; continue; }  // x can't be the card to delete
      canon(in, k, clab, corb);
      int cd = -1;
      for (int i = 0; i < k && cd < 0; i++) if (inv[clab[i]] == best) cd = clab[i];
      if (corb[x] == corb[cd]) {  // x is the canonical card to delete: keep
        classes++;
        hist[count_sets(in)]++;
        for (int v = 0; v < P; v++) fputc('0' + in[v], fo);
        fputc('\n', fo);
      }
      in[x] = 0;
    }
  }
  fclose(fo);
  int mx = 0;
  for (int s = 0; s <= NL; s++) if (hist[s]) mx = s;
  fprintf(stderr, "%d %ld %d %ld ", k, classes, mx, hist[mx]);
  for (int s = 0; s <= NL; s++) if (hist[s]) fprintf(stderr, " %d:%ld", s, hist[s]);
  fprintf(stderr, "\n");
  return 0;
}