/* Proposes hyperbolic components of the Mandelbrot set for verify.py to prove.

   This search is not part of the proof and code/run doesn't run it: it only decides which
   components data/components.txt.gz lists, in double precision, and a component it got
   wrong would only fail verify.py's check. It wrote that file with
       cc -O2 -o search code/search.c -lm
       ./search 1e-11 0.0005 128 > components.txt && gzip -9 -n components.txt
   Floating-point details differ between compilers and machines, so a rerun elsewhere can
   propose a slightly different list; any list gives a valid lower bound once proven.

   It grows each component's tree of satellites: the p/q satellite of a component W of
   period n has period n q, is attached where W's multiplier is e^(2 pi i p/q), and has
   radius about |phi'(e^(2 pi i p/q))| / q^2, where phi is the inverse of W's multiplier
   map. Its center is found by Newton's method from that prediction. Satellites whose
   predicted area is below a threshold are skipped. The trees start from the main
   cardioid and from every center Newton's method finds when started on a grid at each
   period where |f_c^k(0)| reaches a new minimum (the period of an atom domain), up to a
   maximum period. No component above a second maximum period is kept: the proof's work
   for a component grows about as the square of its period, and the many components of
   high period hold little area. Components are kept in the closed upper half-plane, and
   verify.py counts the mirror image of each one off the real axis. Each center is written
   to a thousandth of the component's predicted radius, which Newton's method in verify.py
   refines. */
#include <complex.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define KMAX 12
typedef double complex cplx;

static void smul(const cplx *a, const cplx *b, cplx *out, int K) {
  cplx t[KMAX];
  for (int k = 0; k < K; k++) { cplx s = 0; for (int i = 0; i <= k; i++) s += a[i] * b[k - i]; t[k] = s; }
  memcpy(out, t, sizeof(cplx) * K);
}

/* g with f(g(t)) = t mod t^K, for f[0] = 0 and f[1] != 0 */
static void revert(const cplx *f, cplx *g, int K) {
  memset(g, 0, sizeof(cplx) * K);
  g[1] = 1.0 / f[1];
  for (int k = 2; k < K; k++) {
    cplx comp[KMAX] = {0}, p[KMAX] = {0};
    p[0] = 1;
    for (int j = 1; j <= k; j++) { smul(p, g, p, K); for (int i = 0; i < K; i++) comp[i] += f[j] * p[i]; }
    g[k] -= comp[k] / f[1];
  }
}

/* a[1..K-1]: Taylor coefficients of the inverse multiplier map at the center c0 of period n */
static int inverse_multiplier(cplx c0, int n, int K, cplx *a) {
  cplx z[KMAX] = {0}, w[KMAX], cs[KMAX] = {0};
  cs[0] = c0;
  if (K > 1) cs[1] = 1;
  for (int it = 0; it < K; it++) {
    memcpy(w, z, sizeof w);
    for (int j = 0; j < n; j++) { smul(w, w, w, K); for (int i = 0; i < K; i++) w[i] += cs[i]; }
    memcpy(z, w, sizeof z);
  }
  cplx lam[KMAX] = {0};
  lam[0] = 1;
  memcpy(w, z, sizeof w);
  for (int j = 0; j < n; j++) {
    cplx two_w[KMAX];
    for (int i = 0; i < K; i++) two_w[i] = 2 * w[i];
    smul(lam, two_w, lam, K);
    smul(w, w, w, K);
    for (int i = 0; i < K; i++) w[i] += cs[i];
  }
  if (cabs(lam[0]) > 1e-6) return 0;
  lam[0] = 0;
  if (cabs(lam[1]) == 0 || !isfinite(creal(lam[1]))) return 0;
  revert(lam, a, K);
  return 1;
}

static double series_area(const cplx *a, int K) {
  double s = 0;
  for (int k = 1; k < K; k++) s += k * creal(a[k] * conj(a[k]));
  return M_PI * s;
}
static cplx sval(const cplx *a, int K, cplx l) { cplx s = 0, p = 1; for (int k = 0; k < K; k++) { s += a[k] * p; p *= l; } return s; }
static cplx sder(const cplx *a, int K, cplx l) { cplx s = 0, p = 1; for (int k = 1; k < K; k++) { s += k * a[k] * p; p *= l; } return s; }

static int nucleus(cplx *c, int n) {
  cplx cc = *c;
  for (int it = 0; it < 80; it++) {
    cplx z = 0, dz = 0;
    for (int j = 0; j < n; j++) { dz = 2 * z * dz + 1; z = z * z + cc; if (cabs(z) > 1e10) return 0; }
    if (dz == 0) return 0;
    cplx step = z / dz;
    cc -= step;
    if (!isfinite(creal(cc)) || !isfinite(cimag(cc))) return 0;
    if (cabs(step) <= 1e-15 * fmax(1.0, cabs(cc))) { *c = cc; return 1; }
  }
  return 0;
}

static int exact_period(cplx c, int n) {
  cplx z = 0;
  for (int j = 1; j <= n; j++) { z = z * z + c; if (j < n && n % j == 0 && cabs(z) < 1e-9) return 0; }
  return cabs(z) < 1e-9;
}

typedef struct { int n; cplx c; double area; } comp_t;
static comp_t *comps;
static size_t ncomps, capcomps;
static uint64_t *htab;
static size_t hcap;

static int seen(int n, cplx c) {
  int64_t x = llround(creal(c) * 1e11), y = llround(fabs(cimag(c)) * 1e11);
  uint64_t k = (uint64_t)n * 0x9E3779B97F4A7C15ULL ^ (uint64_t)x * 0xC2B2AE3D27D4EB4FULL ^ (uint64_t)y * 0x165667B19E3779F9ULL;
  if (!k) k = 1;
  size_t i = k & (hcap - 1);
  while (htab[i]) { if (htab[i] == k) return 1; i = (i + 1) & (hcap - 1); }
  htab[i] = k;
  return 0;
}

static double EPS;
static int MAXPERIOD;
static const int KSER = 10;
static int kfor(double area) { return area > 1e-7 ? KSER : area > 1e-10 ? 6 : 4; }

static size_t *stack;
static size_t nstack, capstack;
static void push(size_t i) {
  if (nstack == capstack) { capstack = capstack ? 2 * capstack : 1 << 16; stack = realloc(stack, capstack * sizeof(size_t)); }
  stack[nstack++] = i;
}

static void add(int n, cplx c, int K) {
  if (n > MAXPERIOD) return;
  if (cimag(c) < 0) c = conj(c);
  if (seen(n, c)) return;
  cplx a[KMAX];
  if (!inverse_multiplier(c, n, K, a)) return;
  if (ncomps == capcomps) { capcomps = capcomps ? 2 * capcomps : 1 << 16; comps = realloc(comps, capcomps * sizeof(comp_t)); }
  comps[ncomps] = (comp_t){n, c, series_area(a, K)};
  push(ncomps++);
}

static int gcd(int a, int b) { while (b) { int t = a % b; a = b; b = t; } return a; }

static void grow(void) {
  while (nstack) {
    comp_t W = comps[stack[--nstack]];
    int K = kfor(W.area);
    cplx a[KMAX];
    if (!inverse_multiplier(W.c, W.n, K, a)) continue;
    a[0] = W.c;
    double dmax = 0;
    for (int t = 0; t < 128; t++) { double v = cabs(sder(a, K, cexp(2 * M_PI * I * t / 128.0))); if (v > dmax) dmax = v; }
    int real_parent = fabs(cimag(W.c)) < 1e-13;
    for (int q = 2; M_PI * pow(dmax / ((double)q * q), 2) >= EPS && (long)W.n * q <= MAXPERIOD; q++) {
      for (int p = 1; p < q; p++) {
        if (gcd(p, q) != 1) continue;
        cplx lam = cexp(2 * M_PI * I * (double)p / q);
        cplx root = sval(a, K, lam), d = sder(a, K, lam);
        double rho = cabs(d) / ((double)q * q);
        if (M_PI * rho * rho < EPS) continue;
        cplx guess = root + rho * (d * lam) / cabs(d * lam);
        if (real_parent && cimag(guess) < -1e-13) continue;  /* its mirror image is the same satellite */
        cplx c = guess;
        if (nucleus(&c, W.n * q) && exact_period(c, W.n * q) && cabs(c - guess) < rho) add(W.n * q, c, kfor(M_PI * rho * rho));
      }
    }
  }
}

static void grid_seeds(double h, int maxp) {
  const double x0 = -2.05, x1 = 0.55, y0 = 0.0, y1 = 1.25;
  long nx = (long)((x1 - x0) / h), ny = (long)((y1 - y0) / h);
  for (long iy = 0; iy < ny; iy++)
    for (long ix = 0; ix < nx; ix++) {
      cplx c0 = (x0 + (ix + 0.5) * h) + I * (y0 + (iy + 0.5) * h), z = 0;
      double best = INFINITY;
      for (int k = 1; k <= maxp; k++) {
        z = z * z + c0;
        double az = cabs(z);
        if (az > 4) break;
        if (az < best) {
          best = az;
          cplx c = c0;
          if (nucleus(&c, k) && exact_period(c, k)) add(k, c, KSER);
        }
      }
    }
}

int main(int argc, char **argv) {
  if (argc != 5) {
    fprintf(stderr, "usage: search <area threshold> <grid spacing> <maximum grid period> <maximum period>\n");
    return 2;
  }
  EPS = atof(argv[1]);
  double h = atof(argv[2]);
  int maxp = atoi(argv[3]);
  MAXPERIOD = atoi(argv[4]);
  hcap = (size_t)1 << 26;
  htab = calloc(hcap, sizeof(uint64_t));
  add(1, 0, KSER);
  grow();
  grid_seeds(h, maxp);
  grow();
  double total = 0;
  for (size_t i = 0; i < ncomps; i++) total += (fabs(cimag(comps[i].c)) < 1e-13 ? 1 : 2) * comps[i].area;
  fprintf(stderr, "%zu components, double-precision area sum %.12f\n", ncomps, total);
  printf("# period, then the real and imaginary parts of an approximate center\n");
  for (size_t i = 0; i < ncomps; i++) {
    double radius = sqrt(comps[i].area / M_PI);
    int digits = (int)ceil(-log10(1e-3 * radius));
    if (digits < 3) digits = 3;
    if (digits > 17) digits = 17;
    double re = creal(comps[i].c), im = cimag(comps[i].c);
    if (fabs(im) < 1e-13) im = 0;      /* centers on the real axis are written as real */
    printf("%d %.*f %.*f\n", comps[i].n, digits, re, digits, im);
  }
  return 0;
}