// --------------------------------------------------------------
//
//       kamm  (German for comb)
//
//       Evaluation of perfect rulers of length L with M markers
//       Compile:   cc -o kamm kamm.c -O
//       Parameters M and L can be entered in a dialogue  (#define DIALOG)
//       or read from the File "kamm.param"             (#define noDIALOG)
//       Format of File "kamm.param":
//
//       M1 L1
//       M2 L2
//
//       Mn Ln
//       0
//
//       M=0 terminates the program.
//       kamm checks whether the markers < 10 and > L-10 have been
//       calculated before, possibly in different order.
//
//       To terminate after one solution is found:  #define FIRST_ONLY
//       To get all solutions:                      #define notFIRST_ONLY
//
//       If all solutions are wanted, it might be necessary to
//       increase MAXPROT and MAXMARK, otherwise the protocol overflows.
//
//       In the noFIRST_ONLY mode two counts are given:
//       First number:   All solutions counted.
//       In parentheses: Symmetric solutions counted only once.
//
//       #define RESULTS_ONLY: Only the counts will be printed.
//       Otherwise all solutions are printed, symmetric ones are
//       marked "<<< sym".
//
//       Klaus Nagel  2003/02/10
//       revised      2005/01/19
//       revised      2009/03/23
//
// --------------------------------------------------------------

#define DIALOG
#define notFIRST_ONLY
#define notRESULTS_ONLY
#define MAXPROT   1000000
#define MAXMARK   22

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <string.h>

int L;
int M;
int Protocoll[MAXPROT][MAXMARK-2] ;  // checking for equivalent solutions
int Prot_pointer ;
int Sol_counter ;     // counter for solutions (equivalent Solution once)
int Sol_counter1 ;    // counter for solutions (non-symmetric twice)
int Limit ;           // L - B(M,2)
int Loops,Loops1;     // counts recursive calls (Units, Billions =1.e9)
int perm[2000];
int Stop;
int mark[0x100000];
double sanz[100];
int schonda;


//    ------------------ make Luschny set ------------------


void mach_lus(int k,int a[], int lus[])
{
   int i,j ;
   for(i=0; i <= L ; i++) lus[i] = 0;
   for(i=1; i < k ; i++)
   for(j=0; j < i ; j++)
   {
      lus[abs(a[i]-a[j])]++ ;
   }
}

// ------------------ sym ------------------

int sym(int a[])
{
    int i,j,index[1100];

   memset(index,0,sizeof index);
   for (i=0; i < M ; i++) index[a[i]] = 1;  // For sorting
   j = L ;
    for (i=0; i < j ; i++)
    {
        if(index[i] != index[j--]) return 2;
    }
    return 1;
}
// ------------------ rep ------------------
// rep: Determines an unique representative from a class of
// equivalent solutions and writes it to the protocol at pointer

void rep(int a[],int pointer)
{
   int i,j,k,b[100],index[1100];


   memset(index,0,sizeof index);
   for(i=0 ; i < M ; i++)
   index[a[i]]++ ;

   j = 0;
   k = 0;
   for (i=1; i <= L ; i++)
   {
      if ( index[i])
      {
         b[M-2-k] = Protocoll[pointer][k] = i-j;
         j = i;
         k++;
      }
   }
   for (i=0; i < M-2 ; i++)
   {
      if(b[i] < Protocoll[pointer][i] ) goto b_is_smaller;
      if(Protocoll[pointer][i] < b[i] ) return ;
   }
   return;
b_is_smaller:
   for (i=0; i < M-2 ; i++) Protocoll[pointer][i] = b[i] ;
   return;
}

// ------------------ Key ------------------

int Key(int k,int a[])
{
   int i,x,m;
   x = 0;
   m = 1 << 19;
   for (i=0; i < k ; i++)
   {
      if(a[i] > 9 && a[i] < L - 9) return(-1) ;
      if(a[i] < 10)
      {
         x |= 1 << a[i] ;
      }
      if(a[i] > L-10)
      {
         x |= m >> (L - a[i]) ;
      }
   }
   return(x);
}

// ------------------ Invert ------------------

int invert(int x)
{
   int y,m ;
   y = 0;
   m = 1 << 19;
   while(x)
   {
      if(x & 1) y |= m;
      m >>= 1;
      x >>= 1 ;
   }
   return(y);
}


// ------------------ Set  Marks ------------------

void set_marks(int k,int a[])
{
   int x;
   x = Key(k,a);
   if(x != -1 )
   {
      mark[x] |= (1<<k);
      mark[invert(x)] |= (1<<k);
   }
}


// ---------- The recursive program Pruef ------------------

// The marks a[0] through a[k-1] are set
// a[k] is to be calculated

void pruef(int k, int a[], int lfeld[])
{
   int lus[1100],i,j,lang,index[1100],p,sum,ii,x,syma ;

   sanz[k]++;
   if ( ++Loops == 1000000000)
   {
      Loops1++;
      Loops = 0;
   }
   if(k > M || Stop != 0) return;
   x = Key(k,a);
   if(k < 10 && x != -1 && mark[x] )
   {
      return;
   }

// Determine all lengths yet measurable
// Excess  (A length can be measured excess times too often)

   sum = 0;
   for ( i=1; i <= L ; i++)
   {
      if(lfeld[i] > 1) sum += lfeld[i] - 1;
   }
   if(sum > Limit) return ;  // Not all lengths can be reached any more

// Look for greatest missing length
   for ( lang = L-2 ; lang > 1 ; lang--)  // The Marks 0,1,L are fixed
   {
      if(lfeld[lang] == 0) goto weiter;
   }

// Solution, because no missing length was found

   rep(a,Prot_pointer) ;  // enter solution in protocol
   for (i=0; i < Prot_pointer; i++)
   {
      for ( j = 0 ; j < M-2 ; j++)
      {
         if (Protocoll[Prot_pointer][j]   != Protocoll[i][j] ) break;
      }
      if(j == M-2)  // No deviation between entries i und Prot_pointer.
                    // Solution is not new.
       {
          schonda++;
          return;
       }
   }

// No equivalent solution found in protocol

   Prot_pointer++;
   if (Prot_pointer >= MAXPROT)
    {
         perror("Protocoll overflow\n");
         exit(0);
    }
   memset(index,0,sizeof index);
   for (i=0; i < M ; i++) index[a[i]] = 1;  // For sorting
   ++Sol_counter;
   syma = sym(a);
   Sol_counter1 += syma;

#ifndef RESULTS_ONLY
   printf(">");
   for (i=0; i <= L ; i++) if(index[i]) printf("%4d",i);
   printf(" %8d",Sol_counter);
   if(syma == 1) printf(" <<< sym");
   printf("\n");
#endif

#ifdef FIRST_ONLY
   Stop = 1;
#endif
   return ;

 weiter:

      memset (index, 0, sizeof index);
      for (i =0; i < k ; i++)
      {
         p = a[i] + lang;
         if ( p < L )
         {
            index[p]++;
         }
         p = a[i] - lang;
         if ( p > 0 )
         {
            index[p]++;
         }
      }
// For all marks satisfieing "lang"
      for(ii=3 ; ii <= L ; ii++ )
      {
         i = perm[ii] ;
         if(index[i])
         {
            a[k] = i;
            memcpy(lus,lfeld,sizeof lus);
            for(j=0 ; j < k; j++)
            {
               lus[abs(i-a[j])]++ ;
            }
            pruef(k+1,a,lus);   // This is the recursive call !!!!!!!!
            if ( k < 9)
             set_marks(k+1,a);
         }
      }
}

// ------------------ Main Program ------------------

int main()
{
    struct timeval tp;
    struct timezone tzp;
   int a[100],lus[1100],i,sec,usec,t;
   double zeit;

#ifndef DIALOG
   FILE* fd;
#endif
   printf("kamm, compiled %s   %s\n",__DATE__,__TIME__);
#ifndef DIALOG
   fd = fopen("kamm.param","r");
   if (!fd) perror("File kamm.param is missing\n");
#endif
    schonda = 0;
anfang:

#ifdef DIALOG
   printf("M und L? ");
   scanf("%d",&M);
   if(M==0) exit(0);
   scanf("%d",&L);
#endif

#ifndef DIALOG
   fscanf(fd,"%d",&M);
   if(M==0) exit(0);
   fscanf(fd,"%d",&L);
#endif
   memset(mark,0,sizeof mark);
   Limit = (M-1)*M/2 - L;
   t = 0;
   a[t++] = 0 ;  // Marks 0 und L necessary for length L
   a[t++] = L ;
   a[t++] = 1 ;  // without loss of generality for length L-1

   mach_lus(3,a,lus);

// Generate permutations starting on the outside working
// toward the centre
   for(i=0; i <= L ; i++)
   {
      if(i&1) perm[i] = L - (i>>1);
      else perm[i] = i>>1 ;
   }

   Sol_counter = 0;
   Sol_counter1 = 0;
   Prot_pointer = 0;
   memset(Protocoll,0,sizeof Protocoll);
   memset(sanz,0,sizeof sanz);
   Loops = Loops1 = 0;
   Stop = 0;

   gettimeofday(&tp, &tzp);
   sec = tp.tv_sec;
   usec = tp.tv_usec;
   pruef(t,a,lus);
   gettimeofday(&tp, &tzp);
   zeit = (tp.tv_sec-sec)+(tp.tv_usec-usec)/1000000. ;

   printf("\n%6d %09d recursive calls %12.3lf seconds\n",Loops1,Loops,zeit);
#ifdef FIRST_ONLY

   if(Sol_counter)
        printf("\n M = %3d   L = %4d has at least one solution\n\n",M,L);
   else
        printf("\n M = %3d   L = %4d has no solution\n\n",M,L);
#endif
#ifndef FIRST_ONLY
   printf("\n M = %3d   L = %4d has %8d (%8d) solutions\n\n",M,L,2*Sol_counter,Sol_counter1);
#endif
    printf("%8d times existing solutions encountered.\n",schonda);
#ifdef TEST
   for (i=3 ; i <= M ; i++)
   {
      printf("%3d %12.0lf\n",i,sanz[i]);
   }
#endif
goto anfang;
}