// clang Ge.c Gb.o -Wmost -O3
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
typedef signed char co;
int Lq(co *); int Gq(int);
int const W[12];
static int weight(int x) {int c=0; while(x){++c; x &= x-1;} return c;}
static co v[24], u[24]; int nn[8]; int qw=0;
co sph[196560][24];
void pr(int y){for(int q=0; q<24; ++q) printf("%2d ", sph[y][q]); printf("\n");}
void p(co * v){int lq = Lq(v);
  if(0) {int s=0; for(int e=0; e<24; ++e) s += v[e]*v[e]; 
   if((s!=32) || !lq) {
     for(int k=0; k<24; ++k) printf("%3d", v[k]);
     printf(" pl %d\n", lq);}}
   else {for(int j=0; j<24; ++j) sph[qw][j] = v[j]; ++qw;}}
static void vo(int j){for(int k=j+1; k<24; ++k){
  v[k] = 4; p(v); v[k] = -4; p(v); v[k] = 0;}}
static void G(int k, int P){
  if(k){int j=nn[8-k]; v[j] = 2; G(k-1, P); v[j] = -2; G(k-1, !P); v[j] = 0;}
  else if(P) p(v);}
static void g(int i, int b){int j=i-1;
  if(i) {g(j, b); g(j, b^W[j]);}
  else {assert(Gq(b)); if(weight(b)==8) {int y=8;
    for(int a=0; a<24; ++a) if((b>>a)&1) nn[--y] = 23-a; G(8, 1);}
    for(int j=0; j<24; ++j) u[j] = (b>>(23-j))&1?1:-1;
    for(int j=0; j<24; ++j) {co q=u[j]; u[j] = -3*u[j]; p(u); u[j]=q;}}}
int ip(int i, int k){int t=0;
   for(int j=0; j<24; ++j) t += sph[i][j]*sph[k][j]; return t;}

int main(){for(int j=0; j<24; ++j) v[j]=0;
  for(int j=0; j<24; ++j) {v[j] = 4; vo(j); v[j] = -4; vo(j); v[j] = 0;}
  g(12, 0); assert(qw==196560);
// We have computed coordinate of all central neighbors in sph.
  if(0){ // For each neighbor compute distribution
          // of inner products with all other neighbors.
          // this distribution should match H!
    static int H[9] = {1, 0, 4600, 47104, 93150, 47104, 4600, 0, 1};
    for(int k=0; k<196560; ++k) {int h[9]; for(int j=0; j<9; ++j) h[j]=0;
    for(int j=0; j<196560; ++j) {int p = ip(j, k);
       assert(abs(p) < 33 && p % 8 == 0 ); ++h[4+(p>>3)];}
    if(!(k%1000)) printf("%d\n", k/1000);
    if(0) {for(int j=0; j<9; ++j) printf("%d ", h[j]); printf("\n");}
    else for(int j=0; j<9; ++j) assert(H[j]==h[j]);}}
  else {int hs[25]; for(int j=0; j<25; ++j) hs[j]=0;
    for (int z=0; z<1000; ++z) {
   // lmp will hold the indexes (for sph) of the growing clump.
   // Origin ball inplicitly included.
      int lmp[24]; lmp[0] = random()%196560;
   // We now have a 2-clump.
      int lx = 1; // how many balls in lmp
      int nn[4600], nnx=0; // the shrinking set of candidates
      for(int n=0; n<196560; ++n) if(ip(lmp[0], n) == 16) nn[nnx++]=n;
      assert(nnx==4600);
   // Now nn hods initial nnx candidates.
   // Each iteration of the while loop below shrinks the candidate list.
      while(nnx){ // Loop Invariant: 
     // for j=0 to j=nnx, nn[j] holds the index into sph of
     // ball j among those nnx balls touching all of those named
     // in lmp[0] thru lmp[lx-1].     
      lmp[lx++] = nn[random()%nnx];
   // now we have a (lx+1) clump: the origin, lmp[0], ... lmp[lx-1].
      {int nt=0; for(int k=0; k<nnx; ++k)
        if(ip(nn[k], lmp[lx-1]) == 16) nn[nt++] = nn[k];
        if(z<100) printf("%d ", nnx); nnx = nt;}}
    if(z<100) printf("\n%3d lx=%d\n", z, lx);
    assert(lx<25); ++hs[lx];
    if(z<20) for(int k=0; k<lx; ++k) pr(lmp[k]);}
  for(int j=0; j<25; ++j) if(hs[j]) printf("%3d %5d\n", j, hs[j]);
    printf("\n");}
  return 0;}