/* C6: paint off, one pixel, decay y_t = a y_{t-1} + (1-a) b_t. Exhaustive search over every periodic 0/1
   sequence of period L <= LMAX (DFS, exact steady state), for 32 levels f = (j+.5)/32: the least long-run
   mean of (y - f)^2 found is an attained upper bound on the optimum. Prints RMS over levels and per level.
   build/run: cc -O2 -o c6 c6_periodic.c -lm && ./c6 <tau_ms> <LMAX> */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NF 32
static double a, best[NF], fs[NF], apow[64];
static int L;
static void leaf(double c, double P, double C1, double C2, double S0, double Q) {
  double y0 = c / (1 - apow[L]);
  for (int j = 0; j < NF; j++) {
    double f = fs[j];
    double cost = (S0 * y0 * y0 + 2 * y0 * (P - f * Q) + C2 - 2 * f * C1 + L * f * f) / L;
    if (cost < best[j]) best[j] = cost;
  }
}
/* t = number of bits placed; c = c_t (state from y0 = 0); sums over t' = 1..t of a^t' c_t', c_t', c_t'^2 */
static void dfs(int t, double c, double P, double C1, double C2, double S0, double Q) {
  if (t == L) { leaf(c, P, C1, C2, S0, Q); return; }
  for (int b = 0; b < 2; b++) {
    if (t == 0 && b == 0 && L > 1) {} /* allow all; rotations are duplicates but harmless */
    double cn = a * c + (1 - a) * b; double at = apow[t + 1];
    dfs(t + 1, cn, P + at * cn, C1 + cn, C2 + cn * cn, S0 + at * at, Q + at);
  }
}
int main(int argc, char **argv) {
  double ms = atof(argv[1]); int LMAX = atoi(argv[2]);
  a = exp(-1000.0 / 60.0 / ms);
  apow[0] = 1; for (int i = 1; i < 64; i++) apow[i] = apow[i - 1] * a;
  for (int j = 0; j < NF; j++) { fs[j] = (j + .5) / NF; best[j] = 1e9; }
  for (L = 1; L <= LMAX; L++) dfs(0, 0, 0, 0, 0, 0, 0);
  double m = 0; for (int j = 0; j < NF; j++) m += best[j]; m /= NF;
  printf("decay %g ms (a=%.5f), periods <= %d: RMS %.5f\n", ms, a, LMAX, sqrt(m));
  /* flip every frame, exactly: */
  double fl = 0; for (int j = 0; j < NF; j++) { double y0 = (1 - a) / (1 - a * a), y1 = a * y0; double f = fs[j];
      /* flip: q alternates 1,0 when f in (0,1) per rounding with d alternating ... */ (void)y0; (void)y1; (void)f; }
  for (int j = 0; j < NF; j++) printf(" %.6f", best[j]); printf("\n");
  return 0;
}