// Best periodic threshold tile (m x m, thresholds a permutation of (r+.5)/m^2) for a filter:
// minimise mean over tile pixels of E_p = integral_0^1 (S_p(f) - f)^2 df, S_p(f) = sum_q h_q [d(p+q) >= 1-f].
// Simulated annealing over swaps; exact E_p by sorting the window's thresholds. The filter is read from
// k5.txt (python3 kernel.py > k5.txt: Starlock's pass per tap). Usage: ./anneal m steps T0 seed [start-tile]
// (steps 0 just scores the start tile). Prints the best tile's RMS and the tile (ranks 0..m*m-1).
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
int T, m, *dx, *dy; double *h; double *d; double *E;
typedef struct { double t, w; } tw;
static int cmp(const void *a, const void *b) { double x = ((tw*)a)->t, y = ((tw*)b)->t; return x < y ? -1 : x > y; }
static tw buf[1024];
double pixerr(int px, int py) {
  for (int i = 0; i < T; i++) { int x = ((px + dx[i]) % m + m) % m, y = ((py + dy[i]) % m + m) % m; buf[i].t = 1 - d[y * m + x]; buf[i].w = h[i]; }
  qsort(buf, T, sizeof(tw), cmp);
  double W = 0, prev = 0, s = 0;
  for (int i = 0; i <= T; i++) {
    double nx = i < T ? buf[i].t : 1.0;
    // integral of (W - f)^2 over [prev, nx]
    double a = W - prev, b = W - nx; s += (a * a * a - b * b * b) / 3.0;
    if (i < T) { W += buf[i].w; prev = nx; }
  }
  return s;
}
unsigned long long rs = 88172645463325252ULL;
double rnd() { rs ^= rs << 13; rs ^= rs >> 7; rs ^= rs << 17; return (rs >> 11) * (1.0 / 9007199254740992.0); }
int main(int argc, char **argv) {
  m = atoi(argv[1]); long steps = atol(argv[2]); double T0 = atof(argv[3]); rs += atoi(argv[4]);
  FILE *f = fopen("k5.txt", "r"); int cap = 1024; dx = malloc(cap * 4); dy = malloc(cap * 4); h = malloc(cap * 8); T = 0;
  while (fscanf(f, "%d %d %lf", &dx[T], &dy[T], &h[T]) == 3) T++;
  int n = m * m; d = malloc(n * 8); E = malloc(n * 8);
  // start: optional init file (thresholds), else random permutation
  int *perm = malloc(n * 4); for (int i = 0; i < n; i++) perm[i] = i;
  for (int i = n - 1; i > 0; i--) { int j = (int)(rnd() * (i + 1)); int t = perm[i]; perm[i] = perm[j]; perm[j] = t; }
  if (argc > 5) { FILE *g = fopen(argv[5], "r"); for (int i = 0; i < n; i++) fscanf(g, "%d", &perm[i]); fclose(g); }
  for (int i = 0; i < n; i++) d[i] = (perm[i] + .5) / n;
  double tot = 0; for (int i = 0; i < n; i++) { E[i] = pixerr(i % m, i / m); tot += E[i]; }
  // which pixels see location L: p = L - q
  int *aff = malloc(2 * T * 4); double *newE = malloc(2 * T * 8); char *mark = calloc(n, 1);
  double best = tot; int *bestp = malloc(n * 4); memcpy(bestp, perm, n * 4);
  for (long s = 0; s < steps; s++) {
    double temp = T0 * (1 - (double)s / steps);
    int A = (int)(rnd() * n), B = (int)(rnd() * n); if (A == B) continue;
    double t = d[A]; d[A] = d[B]; d[B] = t;
    int na = 0;
    for (int k = 0; k < 2; k++) { int L = k ? B : A, lx = L % m, ly = L / m;
      for (int i = 0; i < T; i++) { int x = ((lx - dx[i]) % m + m) % m, y = ((ly - dy[i]) % m + m) % m, p = y * m + x; if (!mark[p]) { mark[p] = 1; aff[na++] = p; } } }
    double delta = 0;
    for (int i = 0; i < na; i++) { int p = aff[i]; newE[i] = pixerr(p % m, p / m); delta += newE[i] - E[p]; mark[p] = 0; }
    if (delta <= 0 || (temp > 0 && rnd() < exp(-delta / (temp * n)))) {
      for (int i = 0; i < na; i++) E[aff[i]] = newE[i];
      tot += delta; int tp = perm[A]; perm[A] = perm[B]; perm[B] = tp;
      if (tot < best) { best = tot; memcpy(bestp, perm, n * 4); }
    } else { t = d[A]; d[A] = d[B]; d[B] = t; }
    if (steps >= 10 && s % (steps / 10) == 0) fprintf(stderr, "step %ld rms %.5f best %.5f\n", s, sqrt(tot / n), sqrt(best / n));
  }
  printf("m %d best rms %.6f\n", m, sqrt(best / n));
  for (int i = 0; i < n; i++) printf("%d%c", bestp[i], (i % m == m - 1) ? '\n' : ' ');
  return 0;
}