diff --git a/historique/historique.c b/historique/historique.c
new file mode 100644
index 0000000..4a71dc1
--- /dev/null
+++ b/historique/historique.c
@@ -0,0 +1,235 @@
+#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;
+}
diff --git a/historique/historique_21960.c b/historique/historique_21960.c
new file mode 100644
index 0000000..8c33b7c
--- /dev/null
+++ b/historique/historique_21960.c
@@ -0,0 +1,232 @@
+#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;
+
+  seed = 1571773356;
+  srand(seed);
+  printf("seed: %d\n", seed);
+  pos = 0;
+  for(i = 0; i < 16; ++i)
+  {
+    buf[i] = rand() ^ (rand() << 16) ^ (rand() << 31);
+  }
+
+  init(&current);
+
+  result = 0.0;
+
+  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;
+}