// Simulated annealing: max lines in an n-subset of AG(4,3). Usage: ./sa n seed iters restarts
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cmath>
#include <random>
#include <vector>
#include <algorithm>
static int third[81][81];
static inline int enc(int a, int b, int c, int d) { return a * 27 + b * 9 + c * 3 + d; }
int main(int argc, char** argv) {
    int n = atoi(argv[1]); uint64_t seed = atoll(argv[2]); long iters = atol(argv[3]); int restarts = atoi(argv[4]);
    for (int x = 0; x < 81; x++) for (int y = 0; y < 81; y++) {
        int xc[4] = {x / 27, (x / 9) % 3, (x / 3) % 3, x % 3}, yc[4] = {y / 27, (y / 9) % 3, (y / 3) % 3, y % 3}, zc[4];
        for (int i = 0; i < 4; i++) zc[i] = (6 - xc[i] - yc[i]) % 3;
        third[x][y] = enc(zc[0], zc[1], zc[2], zc[3]);
    }
    std::mt19937_64 rng(seed);
    int globalBest = -1; std::vector<int> bestSet;
    for (int r = 0; r < restarts; r++) {
        std::vector<int> perm(81); for (int i = 0; i < 81; i++) perm[i] = i;
        std::shuffle(perm.begin(), perm.end(), rng);
        std::vector<char> in(81, 0); std::vector<int> ins, outs;
        for (int i = 0; i < 81; i++) { if (i < n) { in[perm[i]] = 1; ins.push_back(perm[i]); } else outs.push_back(perm[i]); }
        // degree: deg[p] = number of lines through p with the other two points in the set (p need not be in set)
        auto degOf = [&](int p) { int c = 0; for (int y : ins) if (y != p && in[third[p][y]]) c++; return c / 2; };
        int L = 0; for (int p : ins) L += degOf(p); L /= 3;
        int best = L;
        double T0 = 2.0, T1 = 0.05;
        for (long it = 0; it < iters; it++) {
            double T = T0 * pow(T1 / T0, (double)it / iters);
            int i = rng() % n, j = rng() % (81 - n);
            int p = ins[i], q = outs[j];
            int dp = degOf(p);
            // remove p
            in[p] = 0; ins[i] = q;  // temporarily
            int dq = 0; for (int y : ins) if (y != q && in[third[q][y]]) dq++; dq /= 2;
            int delta = dq - dp;
            if (delta >= 0 || std::exp(delta / T) * 4294967296.0 > (double)(uint32_t)rng()) {
                in[q] = 1; outs[j] = p; L += delta;
                if (L > best) best = L;
                if (L > globalBest) { globalBest = L; bestSet.clear(); for (int y : ins) bestSet.push_back(y); }
            } else { ins[i] = p; in[p] = 1; }
        }
    }
    printf("n=%d best=%d set=", n, globalBest);
    std::sort(bestSet.begin(), bestSet.end());
    for (int y : bestSet) printf("%d ", y);
    printf("\n");
}