// Classify caps (line-free sets) in AG(4,3) up to AGL(4,3), level by level.
// Canonical form: lexicographically largest 81-bit membership string over all ordered affine
// frames (p0,...,pd) of affinely independent points of S (d = affine dimension of S), mapping
// p0->0, p1->e4(idx 1), p2->e3(idx 3), p3->e2(idx 9), p4->e1(idx 27).
// Usage: ./caps maxlevel outdir
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <vector>
#include <string>
#include <unordered_set>
#include <algorithm>
typedef unsigned __int128 u128;
static int ADD[81][81], NEG[81], THIRD[81][81];
static inline int enc(const int* c) { return c[0] * 27 + c[1] * 9 + c[2] * 3 + c[3]; }
static inline void dec(int i, int* c) { c[0] = i / 27; c[1] = (i / 9) % 3; c[2] = (i / 3) % 3; c[3] = i % 3; }
static void init() {
    for (int x = 0; x < 81; x++) {
        int xc[4]; dec(x, xc);
        int nc[4]; for (int i = 0; i < 4; i++) nc[i] = (3 - xc[i]) % 3; NEG[x] = enc(nc);
        for (int y = 0; y < 81; y++) {
            int yc[4]; dec(y, yc); int s[4], t[4];
            for (int i = 0; i < 4; i++) { s[i] = (xc[i] + yc[i]) % 3; t[i] = (6 - xc[i] - yc[i]) % 3; }
            ADD[x][y] = enc(s); THIRD[x][y] = enc(t);
        }
    }
}
static inline int SUB(int x, int y) { return ADD[x][NEG[y]]; }
static inline bool bit(u128 m, int i) { return (m >> i) & 1; }

struct Canon {
    u128 S; int pts[81]; int n;
    u128 best; long nbest;
    // image[j] for j in 0..80: point p0 + c1 u1 + c2 u2 + c3 u3 + c4 u4 where j = c1 + 3c2 + 9c3 + 27c4
    int img[81];
    int dim;  // affine dimension found
    // prefix comparison: compare bits [lo, hi) of candidate (built in cand) against best
    u128 cand;
    // returns -1 if cand prefix < best prefix, 0 equal, 1 greater, over bit range [0, hi)
    int cmpPrefix(int hi) {
        // bits are compared with index 0 most significant: construct masks
        u128 m = (hi >= 128) ? ~(u128)0 : (((u128)1 << hi) - 1);
        u128 a = cand & m, b = best & m;
        if (a == b) return 0;
        // find the lowest differing bit index (most significant in our order)
        u128 d = a ^ b;
        int low;
        uint64_t lo64 = (uint64_t)d;
        if (lo64) low = __builtin_ctzll(lo64); else low = 64 + __builtin_ctzll((uint64_t)(d >> 64));
        return bit(a, low) ? 1 : -1;
    }
    void finish(int hi) {
        // all 81 bits set in cand (bits >= 3^(dim+1)... are zero). Compare whole.
        int c = cmpPrefix(81);
        if (c > 0) { best = cand; nbest = 1; }
        else if (c == 0) nbest++;
    }
    void stage4(int p0, int u1, int u2, int u3, u128 hyp) {
        // choose p4 in S not in hyperplane
        for (int i = 0; i < n; i++) {
            int p4 = pts[i]; if (bit(hyp, p4)) continue;
            int u4 = SUB(p4, p0);
            u128 c = cand;
            for (int j = 27; j < 81; j++) { img[j] = ADD[img[j - 27]][u4]; if (bit(S, img[j])) c |= (u128)1 << j; }
            u128 save = cand; cand = c; finish(81); cand = save;
        }
    }
    void stage3(int p0, int u1, int u2, u128 plane) {
        for (int i = 0; i < n; i++) {
            int p3 = pts[i]; if (bit(plane, p3)) continue;
            int u3 = SUB(p3, p0);
            u128 c = cand; u128 hyp = plane;
            for (int j = 9; j < 27; j++) { img[j] = ADD[img[j - 9]][u3]; hyp |= (u128)1 << img[j]; if (bit(S, img[j])) c |= (u128)1 << j; }
            u128 save = cand; cand = c;
            if (dim == 3) finish(27);
            else if (cmpPrefix(27) >= 0) stage4(p0, u1, u2, u3, hyp);
            cand = save;
        }
    }
    void stage2(int p0, int u1, u128 line) {
        for (int i = 0; i < n; i++) {
            int p2 = pts[i]; if (bit(line, p2)) continue;
            int u2 = SUB(p2, p0);
            u128 c = cand; u128 plane = line;
            for (int j = 3; j < 9; j++) { img[j] = ADD[img[j - 3]][u2]; plane |= (u128)1 << img[j]; if (bit(S, img[j])) c |= (u128)1 << j; }
            u128 save = cand; cand = c;
            if (dim == 2) finish(9);
            else if (cmpPrefix(9) >= 0) stage3(p0, u1, u2, plane);
            cand = save;
        }
    }
    int affdim() {
        // compute affine dimension of S
        if (n == 0) return -1;
        if (n == 1) return 0;
        int p0 = pts[0];
        // gaussian elimination over F3 on difference vectors
        std::vector<std::vector<int>> rows;
        for (int i = 1; i < n; i++) {
            int d = SUB(pts[i], p0); int c[4]; dec(d, c);
            std::vector<int> v(c, c + 4);
            for (auto& r : rows) {
                int piv = -1; for (int k = 0; k < 4; k++) if (r[k]) { piv = k; break; }
                if (v[piv]) { int f = (v[piv] * r[piv]) % 3; /* r[piv] in {1,2}; want v - (v[piv]/r[piv]) r */
                    int inv = r[piv]; /* inverse of 1 is 1, of 2 is 2 */ int mult = (v[piv] * inv) % 3;
                    for (int k = 0; k < 4; k++) v[k] = ((v[k] - mult * r[k]) % 3 + 3) % 3; (void)f; }
            }
            if (std::any_of(v.begin(), v.end(), [](int x) { return x != 0; })) {
                // reduce existing rows? Not needed for rank; but need echelon: ensure new row's pivot differs. Keep it simple: full reduce rows against v.
                int piv = -1; for (int k = 0; k < 4; k++) if (v[k]) { piv = k; break; }
                for (auto& r : rows) if (r[piv]) { int mult = (r[piv] * v[piv]) % 3; for (int k = 0; k < 4; k++) r[k] = ((r[k] - mult * v[k]) % 3 + 3) % 3; }
                rows.push_back(v);
                if (rows.size() == 4) return 4;
            }
        }
        return (int)rows.size();
    }
    u128 run(u128 s, long* aut) {
        S = s; n = 0; for (int i = 0; i < 81; i++) if (bit(s, i)) pts[n++] = i;
        best = 0; nbest = 0; dim = affdim();
        if (dim <= 0) { best = (n == 1) ? (u128)1 : 0; if (aut) *aut = (n == 1) ? 24261120L : 1965150720L; return best; }
        for (int i = 0; i < n; i++) {
            int p0 = pts[i];
            for (int k = 0; k < n; k++) {
                if (k == i) continue;
                int p1 = pts[k]; int u1 = SUB(p1, p0);
                img[0] = p0; img[1] = p1; img[2] = ADD[p1][u1];
                u128 line = ((u128)1 << p0) | ((u128)1 << p1) | ((u128)1 << img[2]);
                cand = (u128)3; if (bit(S, img[2])) cand |= 4;
                if (dim == 1) finish(3);
                else if (cmpPrefix(3) >= 0) stage2(p0, u1, line);
            }
        }
        // |Aut(S)| = (#optimal frames) * |pointwise stabilizer of a dim-d flat in AGL(4,3)|
        static const long fix[5] = {24261120L, 303264L, 3888L, 54L, 1L};
        if (aut) *aut = nbest * fix[dim];
        return best;
    }
};

struct H128 { size_t operator()(const u128& x) const { uint64_t a = (uint64_t)x, b = (uint64_t)(x >> 64); return a * 0x9E3779B97F4A7C15ull ^ (b + 0x632BE59BD9B4E019ull) * 0xC2B2AE3D27D4EB4Full; } };

int main(int argc, char** argv) {
    init();
    int maxlevel = atoi(argv[1]); std::string outdir = argv[2];
    const double AGL = 1965150720.0;
    std::vector<u128> level; level.push_back((u128)1);  // the 1-point cap at origin (canonical)
    Canon C;
    for (int k = 1; k < maxlevel; k++) {
        std::unordered_set<u128, H128> next;
        for (u128 s : level) {
            // blocked = points on a secant of s
            u128 blocked = s;
            int pts[81], n = 0; for (int i = 0; i < 81; i++) if (bit(s, i)) pts[n++] = i;
            for (int i = 0; i < n; i++) for (int j = i + 1; j < n; j++) blocked |= (u128)1 << THIRD[pts[i]][pts[j]];
            for (int x = 0; x < 81; x++) if (!bit(blocked, x)) next.insert(C.run(s | ((u128)1 << x), nullptr));
        }
        level.assign(next.begin(), next.end());
        std::sort(level.begin(), level.end(), [](const u128& a, const u128& b) { return a < b; });
        // orbit-stabilizer check: sum |AGL|/|Aut| = number of labeled caps of size k+1
        double labeled = 0; long autMin = 1L << 40, autMax = 0;
        for (u128 s : level) { long a; C.run(s, &a); labeled += AGL / a; autMin = std::min(autMin, a); autMax = std::max(autMax, a); }
        printf("size %2d: %zu orbits, labeled caps = %.0f, |Aut| in [%ld, %ld]\n", k + 1, level.size(), labeled, autMin, autMax);
        fflush(stdout);
        if (k + 1 >= 9) {
            FILE* f = fopen((outdir + "/caps" + std::to_string(k + 1) + ".txt").c_str(), "w");
            for (u128 s : level) { for (int i = 0; i < 81; i++) if (bit(s, i)) fprintf(f, "%d ", i); fprintf(f, "\n"); }
            fclose(f);
        }
        if (level.empty()) break;
    }
}