#include <errno.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

struct STATE
{
  int8_t code[42];
  int8_t board[81 * 42];
  double score;
};

/*
WELL512 algorithm
http://www.iro.umontreal.ca/~panneton/WELLRNG.html
implemented by Chris Lomont
http://lomont.org/papers/2008/Lomont_PRNG_2008.pdf
*/
uint32_t buf[16];
uint32_t pos;
uint32_t wellrng512()
{
  uint32_t a, b, c, d;
  a = buf[pos];
  c = buf[(pos + 13) & 15];
  b = a ^ c ^ (a << 16) ^ (c << 15);
  c = buf[(pos + 9) & 15];
  c ^= (c >> 11);
  a = buf[pos] = b ^ c;
  d = a ^ ((a << 5) & 0xDA442D24UL);
  pos = (pos + 15) & 15;
  a = buf[pos];
  buf[pos] = a ^ b ^ d ^ (a << 2) ^ (b << 18) ^ (c << 28);
  return buf[pos];
}

double score(struct STATE *state)
{
  double s;
  int8_t i, j, k, l, dx, dy, lp, lm, o, x, y, xn, yn;
  int8_t dp, dm, fp, fm, gp, gm, maxp, maxm, minp, minm;
  int8_t ma[81], mp[81], mm[81], np[81], nm[81];

  s = 0.0;
  fp = 0;
  fm = 0;
  memset(ma, 0, 81);
  memset(mp, 0, 81);
  memset(mm, 0, 81);

  for(j = 0; j < 42; ++j)
  {
    i = state->code[j];
    l = abs(i) - 1;
    k = (i > 0) - (i < 0);

    gp = fp;
    gm = fm;

    memcpy(np, mp, 81);
    memcpy(nm, mm, 81);

    memcpy(state->board + j * 81, ma, 81);

    for(i = 0; i < 81; ++i)
    {
      o = ma[i];
      if(i == l && k != 0)
      {
        if(ma[i] == 0) ma[i] = 5 * k;
      }
      else
      {
        ma[i] += (np[i] == 3) * ((ma[i] == 0) * 5 - (ma[i] < 0) * ma[i]) - (nm[i] == 3) * ((ma[i] == 0) * 5 + (ma[i] > 0) * ma[i]);
      }
      if(ma[i] == o && o != 0)
      {
        lp = np[i] == 3 || np[i] == 4 ? abs(ma[i]) < 5 : -1;
        lm = nm[i] == 3 || nm[i] == 2 ? abs(ma[i]) < 5 : -1;
        ma[i] += (ma[i] > 0) * lp - (ma[i] < 0) * lm;
      }
      if(ma[i] != o && (ma[i] == 0 || o == 0))
      {
        dp = (ma[i] > 0) - (o > 0);
        dm = (ma[i] < 0) - (o < 0);
        fp += dp;
        fm += dm;
        x = i % 9;
        y = i / 9;
        for(dx = -1; dx < 2; ++dx)
        {
          for(dy = -1; dy < 2; ++dy)
          {
            if(dx == 0 && dy == 0) continue;
            yn = y + dx;
            yn = yn < 0 ? 8 : yn > 8 ? 0 : yn;
            xn = x + dy;
            xn = xn < 0 ? 8 : xn > 8 ? 0 : xn;
            mp[yn * 9 + xn] += dp;
            mm[yn * 9 + xn] += dm;
          }
        }
      }
    }
    maxp = fp > gp ? fp : gp;
    maxm = fm > gm ? fm : gm;
    minp = fp < gp ? fp : gp;
    minm = fm < gm ? fm : gm;
    if(maxp != 0 && maxm != 0)
    {
      s = s / (1 + abs(k) * log10(sqrt(j + 1))) + 1.0 * fp * fm * minp * minm / maxp / maxm;
    }
  }
  return s;
}

void init(struct STATE *state)
{
  memset(state->code, 0, 42);

  state->code[0] = 1;

  memset(state->board, 0, 81 * 42);

  state->score = score(state);
}

void random_walk(struct STATE *state)
{
  int8_t i, imax, j, l;
  double s, smax;

  for(l = 0; l < 2; ++l)
  {
    while(1)
    {
      i = wellrng512() % 163 - 81;
      j = wellrng512() % 10 + 1;

      if(i == 0 || state->board[j * 81 + abs(i) - 1] == 0)
      {
        state->code[j] = i;
        break;
      }
    }
  }

  for(l = 0; l < 1; ++l)
  {
    for(j = 1; j < 12; ++j)
    {
      imax = 0;
      smax = 0.0;
      for(i = -81; i < 82; ++i)
      {
        state->code[j] = i;
        s = score(state);
        if(smax < s)
        {
          imax = i;
          smax = s;
        }
      }
      state->code[j] = imax;
      state->score = smax;
    }
  }
}

void print(struct STATE *state)
{
  int8_t i, j, k, l;

  printf("%.9f ", state->score);
  for(j = 0; j < 42; ++j)
  {
    i = state->code[j];
    l = abs(i) - 1;
    k = (i > 0) - (i < 0);
    printf("(%d, %d), ", k * (l % 9 + 1), k * (l / 9 + 1));
  }
  printf("\n");
}

int main(int argc, char *argv[])
{
  struct STATE current, next;
  double delta, result, t;
  unsigned seed;
  int8_t i;

  result = 0.0;

  while(1)
  {
    seed = time(NULL);
    srand(seed);
    printf("seed: %d\n", seed);
    pos = 0;
    for(i = 0; i < 16; ++i)
    {
      buf[i] = rand() ^ (rand() << 16) ^ (rand() << 31);
    }

    init(&current);

    t = 1.0;

    while(t > 0.001)
    {
      t *= 0.999;

      next = current;

      random_walk(&next);

      if(result < next.score)
      {
        result = next.score;
        print(&next);
      }

      delta = next.score - current.score;
      if(delta >= 0 || wellrng512() < RAND_MAX * exp(delta / t))
      {
        current = next;
      }
    }
  }

  return EXIT_SUCCESS;
}