#include <stdio.h>
#include <gmp.h>
#define GODFREY 0
#define EXT 1
#if GODFREY
  #define CHECK_BIT(var,pos) (((var) & (pos)) ? 2LL : -2LL)
  #define CHECK_BIT_G(var,pos) (((var) & (pos)) ? 1LL : -1LL)
#else
  #define CHECK_BIT(var,pos) (((var) & (pos)) ? 1LL : 0LL)
#endif
#define n 9
#if EXT
  #define NPOS (2*n+1)
#else
  #define NPOS (2*n)
#endif

mpz_t bigsum;
mpz_t tmp;
unsigned int multiproc= 0; // are we running as one of several processes?
unsigned int num_proc= 0; // log_2 of number of processes
unsigned int proc_id= 0; // id of this proces. Range 0 .. 2^num_proc - 1
unsigned int j, k, p, odd, ones;
unsigned long long int swp,min, inner, outer, sum, len, start, abortmul;
unsigned long long int gray;
unsigned long long int bigN = 1LL<<NPOS;
long long int prod[n];


inline void calc(unsigned long long int gray)
{
  for (k=1; k <= n; k++)
  {
    sum = 0;
    swp = 1LL;
    for(j=1; j <= (NPOS-k); j++, swp<<=1) // n^2 calc of sum
      #if GODFREY
	sum += CHECK_BIT_G(gray,swp)*CHECK_BIT_G(gray,swp<<k);
      #else
	if(gray & swp) sum += CHECK_BIT(gray,swp<<k);
      #endif
    prod[k-1] = sum;
  }
  #if EXT
    ones = 0;
    for (k=0; k<NPOS; k++) if (gray & (1LL<<k)) ones++;
  #endif
}

inline void updatesum(unsigned int odd, int times)
{
  #if EXT
    mpz_set_si(tmp, -ones);
  #else
    mpz_set_si(tmp,1);
  #endif
  for(k=0; k < n; k++)
  {
    if(prod[k] == 0LL) 
    {
      return;
    }
    mpz_mul_si(tmp,tmp,prod[k]);
  }
  for (k=0; k < times; k++)
  {
    if(odd) mpz_sub(bigsum,bigsum,tmp);
    else mpz_add(bigsum,bigsum,tmp);
  }
}

int main(int argc, char *argv[])
{
  mpz_init(bigsum);
  mpz_init(tmp);
  odd = 0;
  int arg= 1;
  while(arg<argc && *argv[arg]=='-') 
  { 
    switch(*(argv[arg]+1))
    { 
      case 'p':
	++arg;
	multiproc= 1;
	if (sscanf(argv[arg], "%d/%d", &proc_id, &num_proc) != 2) return -1;
	break;
      default:
	return -1;
    }
    ++arg;
  }

  for (p=0; p < n+1; p++) // phases
  {
    printf("\nPhase %d\n", p);
    for (outer=0; outer < 1ULL<<p; outer++)
    {
      gray = 0ULL;
      if(p!=n) 
	gray |= 1LL<<p; // set bit at pos p fixed

      for (j=0; j < p; j++)
	if(outer & 1LL<<j) 
	  gray |= 1LL<<j | 1LL<<(NPOS-j-1);

      if (proc_id==0)
      {
	calc(gray);
	odd = !odd;
	if(p<n)
	{
	  if(p == n-1) odd = 1;
	  updatesum(odd,2);
	}
	else
	  updatesum(0,1);
      }
      
      if(p<n)
      {
	// find length to calculate for this process
	if (n-p>num_proc) len = 1LL<<(2*(n-p-1)-num_proc+EXT);
	else
	{
	  if (proc_id>0) continue;
	  len = 1LL<<(2*(n-p-1)+EXT);
	}

	start = proc_id * len; // starting position used in multiproc
	if (proc_id!=0)
	{
	  start--;
	  len++;
	}
	gray ^= (start^(start>>1LL))<<(p+1); // initialise gray code
	calc(gray);
	for (inner=start+1; inner < len+start; inner++)
	{
	  min = inner & (-inner); // bit changed in gray code
	  min<<= p+1; // shift according to phase
	  gray ^= min; 
	  if(min&gray) // 1
	  {
	    ones++;
	    swp=min>>1;
	    for (k=0; k < n && swp; k++, swp>>=1) // updating bits to the left of min
	    {
	      prod[k] += CHECK_BIT(gray,swp);
	    }

	    swp=min<<1;
	    for (k=0; k < n && swp < bigN; k++, swp<<=1) // updating bits to the right of min
	    {
	      prod[k] += CHECK_BIT(gray,swp);
	    }
	  }
	  else // 0
	  {
	    ones--;
	    swp=min>>1;
	    for (k=0; k < n && swp; k++, swp>>=1)
	    {
	      prod[k] -= CHECK_BIT(gray,swp);
	    }
	    
	    swp=min<<1;
	    for (k=0; k < n && swp < bigN; k++, swp<<=1)
	    {
	      prod[k] -= CHECK_BIT(gray,swp);
	    }
	  }
	  odd = !odd;
	  updatesum(odd,2);
	}
      }
      #if EXT
	if(p==n && proc_id==0)
	{
	  gray |= 1LL<<n;
	  calc(gray);
	  updatesum(1,1);
	}
      #endif
    }
  }
  #if GODFREY
    gmp_printf("bigSum %Zd\n", bigsum);
    mpz_fdiv_q_2exp(bigsum, bigsum, NPOS);
  #endif
    if(multiproc) printf("process (%d/%d) - ", proc_id, num_proc);
  #if EXT
    printf("enxtended ");
  #endif
    gmp_printf("skolem sequences of order %d: %Zd\n",n, bigsum);
  return 0;
}