// The best threshold array with a given period lattice, for every period lattice of index n in a range:
// a periodic dither is an ordering of the n cosets of its period lattice Lam (one level per coset, L = n).
// Its exact error over flat levels is
//   E = 1/(12 n^2) + (1/n^2) sum_{p,q} min(i_p, i_q) (K_Lam(q - p) - 1/n)
// (i_p in 0..n-1 the level of coset p; K = h*h the filter's autocorrelation folded onto the cosets),
// since sum_{k=1}^{n-1} [i_p >= k][i_q >= k] = min(i_p, i_q). We anneal the ordering for each lattice
// Lam = <(a, 0), (c, d)> (Hermite normal form, ad = n, 0 <= c < a), starting from the best rank-1
// ordering when Lam's quotient is cyclic, and report the best.
//   cc -O3 -o cosets cosets.c -lm && ./cosets nmin nmax moves seed < k5.txt
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

static int T; static int tdx[512], tdy[512]; static double tw[512];
static int Kn; static int kdx[2048], kdy[2048]; static double kw[2048];
static unsigned long long rs = 0x9E3779B97F4A7C15ULL;
static double rnd(void) { rs ^= rs << 13; rs ^= rs >> 7; rs ^= rs << 17; return (rs >> 11) * (1.0 / 9007199254740992.0); }

static int A_, C_, D_, n_;
static int coset(int x, int y) { // reduce (x, y) mod Lam = <(A,0),(C,D)> to an index in [0, n)
  int k = (int)floor((double)y / D_); y -= k * D_; x -= k * C_; x %= A_; if (x < 0) x += A_; return y * A_ + x;
}

static double energy(const int *lev, const double *KL, int n) { // E as above
  double s = 0;
  for (int p = 0; p < n; p++) for (int q = 0; q < n; q++) { int m = lev[p] < lev[q] ? lev[p] : lev[q]; s += m * (KL[p * n + q] - 1.0 / n); }
  return 1.0 / (12.0 * n * n) + s / ((double)n * n);
}

int main(int argc, char **argv) {
  int nmin = atoi(argv[1]), nmax = atoi(argv[2]); long moves = atol(argv[3]); rs += atoi(argv[4]);
  while (scanf("%d %d %lf", &tdx[T], &tdy[T], &tw[T]) == 3) T++;
  // K = autocorrelation of the taps
  for (int a = 0; a < T; a++) for (int b = 0; b < T; b++) {
    int dx = tdx[b] - tdx[a], dy = tdy[b] - tdy[a], f = -1;
    for (int k = 0; k < Kn; k++) if (kdx[k] == dx && kdy[k] == dy) { f = k; break; }
    if (f < 0) { f = Kn++; kdx[f] = dx; kdy[f] = dy; kw[f] = 0; }
    kw[f] += tw[a] * tw[b];
  }
  double gbest = 1e9;
  for (int n = nmin; n <= nmax; n++) {
    double *KL = malloc(sizeof(double) * n * n), *KC = malloc(sizeof(double) * n);
    int *lev = malloc(sizeof(int) * n), *best = malloc(sizeof(int) * n), *at = malloc(sizeof(int) * n);
    int *px = malloc(sizeof(int) * n), *py = malloc(sizeof(int) * n);
    double nbest = 1e9; int ba = 0, bc = 0, bd = 0;
    for (int a = 1; a <= n; a++) {
      if (n % a) continue;
      int d = n / a;
      for (int c = 0; c < a; c++) {
        A_ = a; C_ = c; D_ = d; n_ = n;
        for (int i = 0; i < n; i++) { px[i] = i % a; py[i] = i / a; }
        for (int i = 0; i < n; i++) KC[i] = 0;
        for (int k = 0; k < Kn; k++) KC[coset(kdx[k], kdy[k])] += kw[k];
        for (int p = 0; p < n; p++) for (int q = 0; q < n; q++) KL[p * n + q] = KC[coset(px[q] - px[p], py[q] - py[p])];
        // start: the best rank-1 ordering i = (alpha . p) mod n if the quotient is cyclic, else random
        double e0 = 1e9; int got = 0;
        for (int A = 0; A < n; A++) for (int B = 0; B < n; B++) {
          // a linear form constant on cosets: A a = 0 and A c + B d = 0 (mod n)
          if ((A * a) % n || (A * c + B * d) % n) continue;
          // and a bijection onto Z_n
          char *seen = calloc(n, 1); int ok = 1;
          for (int i = 0; i < n && ok; i++) { int v = ((A * px[i] + B * py[i]) % n + n) % n; if (seen[v]) ok = 0; seen[v] = 1; lev[i] = v; }
          free(seen);
          if (!ok) continue;
          double e = energy(lev, KL, n);
          if (e < e0) { e0 = e; memcpy(best, lev, sizeof(int) * n); got = 1; }
        }
        if (got) memcpy(lev, best, sizeof(int) * n);
        else { for (int i = 0; i < n; i++) lev[i] = i; for (int i = n - 1; i > 0; i--) { int j = (int)(rnd() * (i + 1)); int t = lev[i]; lev[i] = lev[j]; lev[j] = t; } }
        double e = energy(lev, KL, n), eb = e; memcpy(best, lev, sizeof(int) * n);
        for (int i = 0; i < n; i++) at[lev[i]] = i;
        // anneal: swap the levels of two cosets; delta in O(n)
        double T0 = 2e-6, T1 = 1e-9;
        for (long s = 0; s < moves; s++) {
          double temp = T0 * pow(T1 / T0, (double)s / moves);
          int u = (int)(rnd() * n), v; int r = lev[u] + (int)(rnd() * 9) - 4; if (r < 0 || r >= n || r == lev[u]) continue; v = at[r];
          int lu = lev[u], lv = lev[v]; double ds = 0;
          for (int q = 0; q < n; q++) {
            if (q == u || q == v) continue;
            int lq = lev[q];
            double du = (double)((lv < lq ? lv : lq) - (lu < lq ? lu : lq)); // u's level becomes lv
            double dv = (double)((lu < lq ? lu : lq) - (lv < lq ? lv : lq));
            ds += du * ((KL[u * n + q] - 1.0 / n) + (KL[q * n + u] - 1.0 / n)) + dv * ((KL[v * n + q] - 1.0 / n) + (KL[q * n + v] - 1.0 / n));
          }
          // diagonal terms (min(l,l) = l) change too; the u-v pair term is symmetric in the swap
          ds += (lv - lu) * (KL[u * n + u] - 1.0 / n) + (lu - lv) * (KL[v * n + v] - 1.0 / n);
          double de = ds / ((double)n * n);
          if (de <= 0 || rnd() < exp(-de / temp)) {
            lev[u] = lv; lev[v] = lu; at[lv] = u; at[lu] = v; e += de;
            if (e < eb - 1e-15) { eb = e; memcpy(best, lev, sizeof(int) * n); }
          }
        }
        double check = energy(best, KL, n);
        if (fabs(check - eb) > 1e-10) fprintf(stderr, "drift %d %d %d: %g vs %g\n", a, c, d, check, eb);
        if (check < nbest) { nbest = check; ba = a; bc = c; bd = d; }
      }
    }
    printf("n %d best rms %.6f  period lattice (%d,0),(%d,%d)\n", n, sqrt(nbest), ba, bc, bd); fflush(stdout);
    if (nbest < gbest) gbest = nbest;
    free(KL); free(KC); free(lev); free(best); free(at); free(px); free(py);
  }
  printf("overall best rms %.6f\n", sqrt(gbest));
  return 0;
}