// Exact per-split maximum M(a,b,c) test against a target t, for a >= 22:
//   any S with class split (a,b,c) and L >= t has L1 >= t - g(a) - g(c) - Tmax, L2 >= t - g(a) - g(b) - Tmax
//   (Tmax = min(ab,bc,ca)); enumerate such S1, S2 (mod translation), and optimize S0 = H0 \ R exactly:
//   L0 = 117 - 13|R| + C(|R|,2) - tri(R),  T = bc - sum_{q in R} r(-q),  r(w) = #{(y,z) in S1 x S2 : y+z = w}.
// Prints the maximum of L0+L1+L2+T found; if < t the split is impossible at target t.
// split3d (this file, from split3c): two sharper facts, for the kill decision only.
//   1. Pigeonhole: every point of the third has at least b + c - 27 representations
//      y + z, so the |R| points removed cost at least |R|(b + c - 27) crossing sets:
//      T <= bc - |R| max(0, b + c - 27). That raises the thresholds on L1 and L2.
//   2. Per pair, T <= bc - (sum of the |R| smallest r(-q)), and tri(R) >= 0, so a pair
//      whose best case is below t is skipped. It stops at the first pair reaching t.
// So it decides max L < t exactly, but when it says "below", the max it prints is
// only over the pairs it didn't skip. Parallel: ./split3d a b c t reps-file part parts
// split3e (from split3d): shape mode. A configuration reaching t only counts if the
// whole layout contains no full hyperplane (any of the 120); slice i is embedded as
// the points (y, i) of F_3^4, so a transversal line y0 + y1 + y2 = 0 is a line there.
// Usage: ./split3c a b c t
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <vector>
#include <algorithm>
static int ADD[27][27], NEG[27], THIRD[27][27];
static const int gmax[28] = {0,0,0,1,1,2,3,5,8,12,12,13,14,16,19,23,26,30,36,41,47,54,62,71,81,92,104,117};
static inline int enc(int a, int b, int c) { return a * 9 + b * 3 + c; }
uint32_t translate(uint32_t m, int v) { uint32_t r = 0; for (int x = 0; x < 27; x++) if (m >> x & 1) r |= 1u << ADD[x][v]; return r; }
bool transl_min(uint32_t m) { for (int v = 1; v < 27; v++) if (translate(m, v) < m) return false; return true; }
static uint64_t HYP[120][2];
static bool has_full_hyperplane(uint32_t s0, uint32_t s1, uint32_t s2) {
    uint64_t lo = 0, hi = 0;  // 81 bits: point (y, i) -> index y + 27 i
    for (int y = 0; y < 27; y++) {
        if (s0 >> y & 1) lo |= 1ull << y;
        if (s1 >> y & 1) lo |= 1ull << (27 + y);
        if (s2 >> y & 1) { int k = 54 + y; if (k < 64) lo |= 1ull << k; else hi |= 1ull << (k - 64); }
    }
    for (int h = 0; h < 120; h++) if ((HYP[h][0] & ~lo) == 0 && (HYP[h][1] & ~hi) == 0) return true;
    return false;
}
int main(int argc, char** argv) {
    int a = atoi(argv[1]), b = atoi(argv[2]), c = atoi(argv[3]), t = atoi(argv[4]); const char* s1file = argc > 5 && argv[5][0] != '-' ? argv[5] : nullptr;
    int part = argc > 7 ? atoi(argv[6]) : 0, parts = argc > 7 ? atoi(argv[7]) : 1;
    for (int x = 0; x < 27; x++) {
        int xc[3] = {x / 9, (x / 3) % 3, x % 3}; NEG[x] = enc((3 - xc[0]) % 3, (3 - xc[1]) % 3, (3 - xc[2]) % 3);
        for (int y = 0; y < 27; y++) { int yc[3] = {y / 9, (y / 3) % 3, y % 3};
            ADD[x][y] = enc((xc[0] + yc[0]) % 3, (xc[1] + yc[1]) % 3, (xc[2] + yc[2]) % 3);
            THIRD[x][y] = enc((6 - xc[0] - yc[0]) % 3, (6 - xc[1] - yc[1]) % 3, (6 - xc[2] - yc[2]) % 3); }
    }
    { int nh = 0;
      for (int w = 1; w < 81; w++) {
        int d[4] = {w / 27, (w / 9) % 3, (w / 3) % 3, w % 3}, lead = 0; while (!d[lead]) lead++;
        if (d[lead] != 1) continue;
        for (int tt = 0; tt < 3; tt++, nh++) { HYP[nh][0] = HYP[nh][1] = 0;
          for (int i = 0; i < 3; i++) for (int y = 0; y < 27; y++) {
            int yc[3] = {y / 9, (y / 3) % 3, y % 3};
            if ((d[0] * yc[0] + d[1] * yc[1] + d[2] * yc[2] + d[3] * i) % 3 == tt) {
              int k = y + 27 * i; if (k < 64) HYP[nh][0] |= 1ull << k; else HYP[nh][1] |= 1ull << (k - 64); } } } }
      if (nh != 120) { printf("hyperplanes: %d\n", nh); return 1; } }
    int rsz0 = 27 - a;
    int Tmax = std::min(a * b, std::min(c * a, b * c - rsz0 * std::max(0, b + c - 27)));
    int minL1 = t - gmax[a] - gmax[c] - Tmax, minL2 = t - gmax[a] - gmax[b] - Tmax;
    int rsz = 27 - a;
    // R subsets of size rsz with their collinear-triple counts
    std::vector<uint32_t> Rs; std::vector<int> Rtri;
    for (uint32_t m = 0; m < (1u << 27); m++) { /* too many; build by recursion instead */ break; }
    {
        std::vector<int> idx(rsz);
        // iterate combinations
        std::vector<int> comb(rsz); for (int i = 0; i < rsz; i++) comb[i] = i;
        while (true) {
            uint32_t m = 0; for (int i = 0; i < rsz; i++) m |= 1u << comb[i];
            int tri = 0; for (int i = 0; i < rsz; i++) for (int j = i + 1; j < rsz; j++) if (m >> THIRD[comb[i]][comb[j]] & 1) tri++;
            Rs.push_back(m); Rtri.push_back(tri / 3);
            int i = rsz - 1; while (i >= 0 && comb[i] == 27 - rsz + i) i--; if (i < 0) break;
            comb[i]++; for (int j = i + 1; j < rsz; j++) comb[j] = comb[j - 1] + 1;
            if (rsz == 0) break;
        }
    }
    int L0base = 117 - 13 * rsz + rsz * (rsz - 1) / 2;
    // candidates S1, S2
    std::vector<uint32_t> B, C; std::vector<int> LB, LC;
    {
        uint32_t mask = 0; int L = 0, n = 0;
        for (uint64_t i = 1; i < (1ull << 27); i++) {
            int p = __builtin_ctzll(i); int cnt = 0; uint32_t m = mask & ~(1u << p);
            while (m) { int y = __builtin_ctz(m); m &= m - 1; if (mask >> THIRD[p][y] & 1) cnt++; } cnt /= 2;
            if (mask >> p & 1) { mask &= ~(1u << p); L -= cnt; n--; } else { mask |= 1u << p; L += cnt; n++; }
            if (!s1file && n == b && L >= minL1 && transl_min(mask)) { B.push_back(mask); LB.push_back(L); }
            if (n == c && L >= minL2 && transl_min(mask)) { C.push_back(mask); LC.push_back(L); }
        }
    }
    if (s1file) { FILE* f = fopen(s1file, "r"); char line[1024];
        while (fgets(line, sizeof line, f)) { uint32_t m = 0; char* p = strtok(line, " \n"); while (p) { m |= 1u << atoi(p); p = strtok(nullptr, " \n"); }
            if (!m) continue; int L = 0; for (int x = 0; x < 27; x++) if (m >> x & 1) for (int y = x + 1; y < 27; y++) if (m >> y & 1) if (m >> THIRD[x][y] & 1) L++; L /= 3;
            if (L >= minL1) { B.push_back(m); LB.push_back(L); } } fclose(f); }
    printf("split (%d,%d,%d) target %d: U=%d; need L1>=%d (%zu cands mod transl), L2>=%d (%zu cands); |R|=%d (%zu subsets)\n",
           a, b, c, t, gmax[a] + gmax[b] + gmax[c] + Tmax, minL1, B.size(), minL2, C.size(), rsz, Rs.size());
    fflush(stdout);
    int best = -1;
    for (size_t i = part; i < B.size(); i += parts) {
        int ys[27], ny = 0; for (int x = 0; x < 27; x++) if (B[i] >> x & 1) ys[ny++] = x;
        for (size_t j = 0; j < C.size(); j++) {
            if (LB[i] + LC[j] + gmax[a] + Tmax < t) continue;
            int r[27] = {0}; for (int x = 0; x < 27; x++) if (C[j] >> x & 1) for (int k = 0; k < ny; k++) r[ADD[x][ys[k]]]++;
            int rn[27]; for (int w = 0; w < 27; w++) rn[w] = r[NEG[w]];   // rn[q] = r(-q)
            int srt[27]; for (int w = 0; w < 27; w++) srt[w] = rn[w];
            std::partial_sort(srt, srt + rsz, srt + 27);
            int minPen = 0; for (int k = 0; k < rsz; k++) minPen += srt[k];
            if (L0base + LB[i] + LC[j] + b * c - minPen < t) continue;   // can't reach t
            int bestPen = 1 << 20;
            for (size_t s = 0; s < Rs.size(); s++) {
                int pen = Rtri[s]; uint32_t m = Rs[s]; while (m) { int q = __builtin_ctz(m); m &= m - 1; pen += rn[q]; }
                if (L0base + LB[i] + LC[j] + b * c - pen >= t && has_full_hyperplane(((1u << 27) - 1) & ~Rs[s], B[i], C[j])) continue;
                if (pen < bestPen) bestPen = pen;
            }
            if (bestPen == (1 << 20)) continue;
            int tot = L0base + LB[i] + LC[j] + b * c - bestPen;
            if (tot > best) best = tot;
            if (best >= t) break;
        }
        if (best >= t) break;
    }
    printf("split (%d,%d,%d): exact max L = %d  -> %s target %d\n", a, b, c, best, best >= t ? "REACHES" : "below", t);
}