// TABU SEARCH ALGORITHM FOR THE JOB-SHOP SCHEDULING PROBLEM
//
// Copyright: E.Nowicki and C.Smutnicki, Wroclaw University of Technology
//
// Source algorithm: E.Nowicki,C.Smutnicki, A FAST TABOO SEARCH ALGORITHM 
//          FOR THE JOB SHOP PROBLEM, Management Science 6, 1966, 797-813
//
// Coded in C by C.Smutnicki (smutnick@ict.pwr.wroc.pl) release 20.12.1992
//   modified 04.10.1994 (improved general outlook, 20-25% faster)
//
// INSTANCE DATA
// job-shop instance is defined by n,m,o,p[1..o],ni[1..o],nj[1..n] where,
// n         - number of jobs,
// m         - number of machines,
// o         - total number of operations,
// p[j]      - the processing time of operation j,
// ni[j]     - the machine dedicated for processing of operation j,
// nj[j]     - the number of operations in job j,
// operations must be indexed from 1 to o in the following order: 
//             by jobs and within a job in the technological order of 
//             processing (e.g. the first job 1,2,.., n[1], the second -
//             n[1]+1,n[1]+2,...,n[1]+n[2], etc.)
// snj[1..n] - array that defines the partition of p[] and ni[] into
//             individual jobs, i.e. operations with indices 
//             snj[j-1]+1,...,snj[j] belong to job j, snj[0]=0,
//             this table is NOT considered as the source data since follows
//             from nj[] so it will be skiped in parameters of functions

// PROCESSING ORDER 
// A solution is defined by the processing order per machines pi[1..o].
// The partition of pi[] into individual permutations (for machines) is
// given by array snk[0..m], where area (pi[j], j=snk[i-1]+1...snk[i]), 
// defines the sequence of operation processed on machine i, snk[0]=0.
// Array snk[] is NOT considered as source data since its value follows from
// m, ni and snj, so it is always skipped in parametrs of functions.

// GRAPH
// Graph G corresponds one-to-one to the processing order pi[].
// The date structure of the graph is adjusted to carry only job-shop, 
// so each node of the graph has at most two immediate successors (also at
// most two immediate predecessors).
// Graph is constructed once at the begin and modified during TS run.
// Graph is defined by three arrays a[1..n], b[1..n] and c[1..n], where
// n    - the number of nodes in the graph (in fact n equals to o),
// a[j] - index of immediate sequential successor of j, zero otherwise
// b[j] - index of immediate technological successor of j, zero otherwise
// c[j] - equals the number of immediate predecessors of j: 0,1 or 2,
// p[j] - weight of node j (array p[] is the same as for the instance data).

// BENCHMARK ft10 (notorious instance of Fisher and Thompson 10 jobs 10 machines)
// 10 10
// 0 29 1 78 2  9 3 36 4 49 5 11 6 62 7 56 8 44 9 21
// 0 43 2 90 4 75 9 11 3 69 1 28 6 46 5 46 7 72 8 30
// 1 91 0 85 3 39 2 74 8 90 5 10 7 12 6 89 9 45 4 33
// 1 81 2 95 0 71 4 99 6  9 8 52 7 85 3 98 9 22 5 43
// 2 14 0  6 1 22 5 61 3 26 4 69 8 21 7 49 9 72 6 53
// 2 84 1  2 5 52 3 95 8 48 9 72 0 47 6 65 4  6 7 25
// 1 46 0 37 3 61 2 13 6 32 5 21 9 32 8 89 7 30 4 55
// 2 31 0 86 1 46 5 74 4 32 6 88 8 19 9 48 7 36 3 79
// 0 76 1 69 3 76 5 51 2 85 9 11 6 40 7 89 4 26 8 74
// 1 85 0 13 2 61 6  7 8 64 9 76 5 47 3 52 4 90 7 45

// PROGRAM TEST
// set accordingly to the source algorithm of Nowicki & Smutnicki 
//   parameters of TS as follows: maxt=8, maxl=5, maxiter=2500, 
//   diter=400, maxd=100, maxc=2, and start from insa() method
// for the instance ft10 we have to obtain the following search trajectory: 
// starting makespan provided by insa() is 994 and then
// makespan - iteration - no of neighbours
//   986           1          7  
//   980           2          6  
//   964           4          5  
//   963           5          5  
//   962          40          2  
//   959         148          5  
//   958        1973          4  
//   957        2096          1  
//   952        2629          3  
//   938        2640          5  
// back jump to solution 938 iteration 5141
// back jump to solution 938 iteration 7241
//   930        7326 

// OTHER NOTES
// all "if" instructions treat zero int as false and non-zero as true
// program was run under borland c++ on a PC

#include <dos.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <alloc.h>
#include <string.h>
#include <values.h>
#include <math.h>

#define MAXUNSINT 65520

#define MS   sizeof(int)
#define LS   sizeof(long)

#define MSX  16              // no of bites on 2b int

int min(int x, int y) { return (x<y)?x:y; }

int max(int x, int y) { return (x>y)?x:y; }

// returns time in seconds from t to the current time moment
float showtime(struct time t)
{ struct time s;
  float ta,tb;

  gettime(&s);
  ta=3600.0*t.ti_hour+60.0*t.ti_min+t.ti_sec+0.01*t.ti_hund;
  tb=3600.0*s.ti_hour+60.0*s.ti_min+s.ti_sec+0.01*s.ti_hund;
  return tb-ta;
}

// own proposal of
// Taillard's uniform[low,high] auxiliary generator for random instances, 
// see http://mscmga.ms.ic.ac.uk/info.html files jobshop?.txt for details
// or Taillard's www page
int unif(long *seed, int low, int high)
{ long k;
  double value_0_1;
  long m=2147483647l, a=16807l, b=127773l, c=2836l;

  k=*seed/b; *seed=a*(*seed%b)-k*c;
  if (*seed < 0) *seed=*seed+m;
  value_0_1=*seed; value_0_1/=m;
  return low+(int)floor(value_0_1*(high-low+1));
}

// swaps two integers
void swap_(int *a, int *b) { int temp; temp=*a; *a=*b; *b=temp; }

// own proposal of
// Taillard's generator for random job-shop instances, 
// see http://mscmga.ms.ic.ac.uk/info.html files jobshop?.txt for details
// or Taillard's www page
void generator(long *t_seed, long *m_seed, int n, int m, int p[], int ni[])
{ int i,j,b=1,e=99;
  for (j=1;j<=n;j++)
    for (i=1;i<=m;i++)
      p[m*(j-1)+i]=unif(t_seed,b,e);
  for (j=1;j<=n;j++) for (i=1;i<=m;i++) ni[m*(j-1)+i]=i;
  for (j=1;j<=n;j++)
    for (i=1;i<=m;i++)
      swap_(ni+m*(j-1)+i,ni+m*(j-1)+unif(m_seed,i,m));
}

// adds an arc to the graph
void _arc(int i, int j, int ab[], int c[]) { ab[i]=j; c[j]++; }

// adds an arc to the extended graph
void _arc_(int i, int j, int ab[], int _ab[], int c[]) { ab[i]=j; _ab[j]=i; c[j]++; }

// removes an arc from the graph
void _darc(int i, int j, int ab[], int c[]) { ab[i]=0; c[j]--; }

// returns the length of the longest path passing through selected node w,
// returns -1 if the graph contains a cycle
int rpq(int n, int w, int a[], int b[], int c[], int p[])
{ int i,j,k,l,x,y,q=0,nn=n+1;
  int z=0,ls=0,lp=1;
  int *cp;
  int *stos,*d;
  
  cp=(int *)malloc(nn*MS); stos=(int *)malloc(nn*MS); d=(int *)calloc(nn,MS);

  for(i=1;i<=n;i++) { cp[i]=c[i]; d[i]=0; if (!cp[i]) stos[++ls]=i; }

  while (lp <= ls)
  { i=stos[lp++]; if (!(i-w)) k=lp;
    x=d[i]+p[i]; if (z < x) z=x;
    j=a[i]; if (j) { if (!--cp[j]) stos[++ls]=j; if (x > d[j]) d[j]=x; }
    j=b[i]; if (j) { if (!--cp[j]) stos[++ls]=j; if (x > d[j]) d[j]=x; }
  }
  if (ls < n) { free(cp); free(stos); free(d); return -1; }
  y=d[w];

  for (i=1;i<=n;d[i++]=0);

  q=x=d[w]+p[w];
  j=a[w]; if (j) if (x > d[j]) d[j]=x;
  j=b[w]; if (j) if (x > d[j]) d[j]=x;
  for (l=k;l<=n;l++)
  { i=stos[l];
    if (d[i])
    { x=d[i]+p[i]; if (x > q) q=x;
      j=a[i]; if (j) if (x > d[j]) d[j]=x;
      j=b[i]; if (j) if (x > d[j]) d[j]=x;
    }
  }
  free(cp); free(stos); free(d); return y+q;
}

// quicksort; returns permutation such that a[pi[i]]<=a[pi[i+1]] for a[n..m]
void sort(int n, int m, int a[], int pi[])
{  int i,j,p;
   int x;

   if (m <= n) return;
   i=n; j=m; x=a[pi[(i+j)/2]];
   do
   {  while (a[pi[i]] < x) i++;
      while (x < a[pi[j]]) j--;
      if (i <= j) { p=pi[j]; pi[j]=pi[i]; pi[i]=p; i++; j--; }
   }  while (i < j);
   sort(n,j,a,pi); sort(i,m,a,pi);
}

// quicksort; returns permutation such that a[pi[i]]>=a[pi[i+1]] for a[n..m]
void _sort(int n, int m, int a[], int pi[])
{  int i,j,p;                  
   int x;

   if (m <= n) return;
   i=n; j=m; x=a[pi[(i+j)/2]];
   do
   {  while (a[pi[i]] > x) i++;
      while (x > a[pi[j]]) j--;
      if (i <= j) { p=pi[j]; pi[j]=pi[i]; pi[i]=p; i++; j--; }
   }  while (i < j);
   _sort(n,j,a,pi); _sort(i,m,a,pi);
}

void mysort(int l, int p[], int pi[])
{
  int i,j,a,b;

  _sort(1,l,p,pi);

// lexicographic order of elements with equal p[] values
  for (i=2;i<=l;i++)
  {
    b=pi[i]; a=p[b]; j=i;
et: if (p[pi[j-1]] > a) goto re;
    if (pi[j-1] < b) goto re;
    pi[j]=pi[j-1]; j=j-1;
    if (j > 1) goto et;
re: pi[j]=b;
  }
}

// swaps elements j and j+1 in the processing order pi[]
// and modifies appropriately graph G
void __swap(int j, int pi[], int a[], int c[])
{ int _i,i,i_,k;

  _i=pi[j-1]; i=pi[j]; i_=pi[j+1];

  a[i]=0; c[i_]--;
  if (a[_i]) { c[i]--; a[_i]=i_; c[i_]++; } if (a[i_]) a[i]=a[i_];
  a[i_]=i; c[i]++;

  pi[j]=i_; pi[j+1]=i;
}

// swaps elements j and j+1 in the processing order pi[]
// and modifies appropriately graph G (invpi[] ???)
void swap(int j, int pi[], int a[], int c[])
{ int _i,i,i_,k;

  _i=pi[j-1]; i=pi[j]; i_=pi[j+1];

  a[i]=0; c[i_]--;
  if (a[_i]) { c[i]--; a[_i]=i_; c[i_]++; } if (a[i_]) a[i]=a[i_];
  a[i_]=i; c[i]++;

  pi[j]=i_; pi[j+1]=i;
}

// swaps elements j and j+1 in the processing order pi[]
// and modifies appropriately extended graph G
void _swap_(int j, int pi[], int a[], int _a[], int c[])
{ int _i,i,i_,k;

  _i=pi[j-1]; i=pi[j]; i_=pi[j+1];

  a[i]=0; c[i_]--;
  if (a[_i]) { c[i]--; a[_i]=i_; c[i_]++; } if (a[i_]) a[i]=a[i_];
  a[i_]=i; c[i]++;

  pi[j]=i_; pi[j+1]=i;
}

// provides ad hoc processing order pi[] for the job-shop instance
void simple(int n, int m, int o, int nj[], int p[], int ni[], int pi[])
{
  int i,j,k;
  int *nk,*card;

  nk=(int *)calloc(m+1,MS); card=(int *)calloc(m+1,MS);

  for (k=1;k<=m;k++) nk[k]=0;

  for (k=1;k<=o;k++) nk[ni[k]]++;
  card[1]=0; for(k=1;k<m;k++) card[k+1]=card[k]+nk[k];

  k=0;
  for (j=1;j<=n;j++)
  {
    for (i=k+1;i<=k+nj[j];i++) pi[++card[ni[i]]]=i;
    k+=nj[j];
  }

  free(nk); free(card);
}

// provides starting processing order pi[] using SPT rule for the job-shop instance
void spt(int n, int m, int o, int nj[], int p[], int ni[], int pi[])
{ int i,j,k,ls=0,pm,ix=0;
  unsigned rm;
  int *nk,*snk,*snj,*car,*stos,*a;
  unsigned *t,*r;

  nk=(int *)calloc(m+1,MS); snk=(int *)calloc(m+1,MS); car=(int *)calloc(m+1,MS);
  snj=(int *)calloc(n+1,MS);
  r=(unsigned *)calloc(o+1,MS); stos=(int *)calloc(o+1,MS); t=(unsigned *)calloc(m+1,MS);
  a=(int *)calloc(o+1,MS);

  for (i=1;i<=m;i++) nk[i]=0;
  for (i=1;i<=o;i++) a[i]=0;

  for (i=1;i<=o;i++) nk[ni[i]]++;
  j=1; for (i=1;i<=m;i++) { snk[i]=j; j+=nk[i]; }
  j=1; for (i=1;i<=n;i++) { snj[i]=j; j+=nj[i]; }
  j=0; for (i=1;i<=m;i++) { car[i]=j; j+=nk[i]; }
  for (i=1;i<=n;i++) for (j=snj[i];j<snj[i]+nj[i]-1;j++) a[j]=j+1;

  for (j=1;j<=n;j++) stos[++ls]=snj[j];
  while (ls >= 1)
  { ix++; if (100*(ix/100)==ix) printf("%8d",ix);   /* !!!!!!!!!!! */
    rm=MAXUNSINT; pm=MAXINT;
    for (i=1;i<=ls;i++) rm=min(rm,r[stos[i]]);
    for (i=1;i<=ls;i++)
      { j=stos[i]; if (!(r[j]-rm)) if (pm > p[j]) { pm=p[j]; k=i; }}
    j=stos[k]; stos[k]=stos[ls]; ls--;
    k=ni[j]; pi[++car[k]]=j; t[k]=r[j]+p[j];
    if (a[j]) { stos[++ls]=j+1; r[j+1]=max(r[j+1],t[k]); }
    for (i=1;i<=o;i++) if (!(ni[i]-k)) r[i]=max(r[i],t[k]);
  }
  free(nk); free(snk);  free(snj); free(car);
  free(r);  free(stos); free(t);   free(a);
}

// provides starting processing order pi[] with accordance to INSA method,
// see E.Nowicki,C.Smutnicki, Management Science 6, 1966, 797-813, for details
void insa(int n, int m, int o, int nj[], int p[], int ni[], int pi[])
{
  int i,j,k,l,s,t,v,nn,ix,op,ma,im,rem;
  int va,opt;
  int *nk,*snk,*snj,*sig,*pip,*car;
  int *a,*b;
  int *c;

  nk=(int *)calloc(m+1,MS); snk=(int *)calloc(m+1,MS); snj=(int *)calloc(n+1,MS);
  sig=(int *)calloc(o+1,MS); pip=(int *)calloc(o+1,MS); car=(int *)calloc(m+1,MS);
  a=(int *)calloc(o+1,MS); b=(int *)calloc(o+1,MS); c=(int *)calloc(o+1,MS);

  for (i=1;i<=m;i++) nk[i]=0;
  for (i=0;i<=o;i++) sig[i]=0;

  for (i=1;i<=o;i++) nk[ni[i]]++;
  j=1; for (i=1;i<=m;i++) { snk[i]=j; j+=nk[i]; }
  j=1; for (i=1;i<=n;i++) { snj[i]=j; j+=nj[i]; }
  j=0; for (i=1;i<=m;i++) { car[i]=j; j+=nk[i]; }

  t=0;
  for (j=1;j<=n;j++)
  {
    s=0; for (i=snj[j];i<=snj[j]+nj[j]-1;i++) s+=p[i];
    if (s > t) { t=s; v=j; }
  }
  l=0;
  for (j=1;j<=n;j++) if (j-v) for (i=snj[j];i<=snj[j]+nj[j]-1;i++) pip[++l]=i;

  nn=o-nj[v]; 
  mysort(nn,p,pip);          
  for (i=snj[v];i<=snj[v]+nj[v]-1;i++) sig[++car[ni[i]]]=i;

  for (i=0;i<=o;i++) { a[i]=0; b[i]=0; c[i]=0; }
  for (j=1;j<=n;j++) for (i=snj[j];i<snj[j]+nj[j]-1;i++) _arc(i,i+1,b,c);

  for (ix=1;ix<=nn;ix++)
  {
     if (10*(ix/10)==ix) printf("%4d\n",ix);          // trace !!!!!!!!!!! 

     op=pip[ix]; ma=ni[op];
     for (i=car[ma];i>=snk[ma];i--) sig[i+1]=sig[i];
     sig[snk[ma]]=op; car[ma]++; _arc(op,sig[snk[ma]+1],a,c);

     opt=MAXINT;
     for (im=snk[ma];im<=car[ma];im++)
     {
       va=rpq(o,op,a,b,c,p); if (va!=-1) if (va < opt) { opt=va; rem=im; }
       if (im!=car[ma]) swap(im,sig,a,c);
     }
     for (i=car[ma]-1;i>=rem;i--) swap(i,sig,a,c);
  }
  for (i=1;i<=o; i++) pi[i]=sig[i];
  free(nk); free(snk); free(snj); free(car); free(sig); free(pip);
  free(a); free(b); free(c);
}

// local function only for lb() 
// returns value t[?] for the operation of job j processed on machine k 
// returns MAXINT if job j has no operation on machine k
int get(int t[], int ni[], int nj[], int snj[], int k, int j)
{ int i;
  for (i=snj[j];i<=snj[j]+nj[j]-1;i++) if (!(ni[i]-k)) return t[i];
  return MAXINT;
}

// returns a lower bound for the job-shop instance in accordance to
// Taillard's proposal
unsigned  lb(int n, int m, int o, int nj[], int p[], int ni[])
{
  int i,j,k,rm,qm;
  unsigned s,ta=0,tb=0;
  int *snj,*r,*q;

  snj=(int *)calloc(n+1,MS); r=(int *)calloc(o+1,MS); q=(int *)calloc(o+1,MS);

  for (i=1;i<=o;i++) r[i]=q[i]=0;

  j=1; for (i=1;i<=n;i++) { snj[i]=j; j+=nj[i]; }
  for (j=1;j<=n;j++)
  {  for (i=snj[j]+1;i<=snj[j]+nj[j]-1;i++) r[i]=r[i-1]+p[i-1];
     for (i=snj[j]+nj[j]-1;i>=snj[j]+1;i--) q[i-1]=q[i]+p[i];
  }

  k=0;
  for (j=1;j<=n;j++)
  {
    s=0; for (i=k+1;i<=k+nj[j];s+=p[i++]); k+=nj[j]; if (s > ta) ta=s;
  }
  for (k=1;k<=m;k++)
  {
    rm=MAXINT; qm=MAXINT; s=0;
    for (i=1;i<=n;i++)
    { rm=min(rm,get(r,ni,nj,snj,k,i));
      if (get(p,ni,nj,snj,k,i)!=MAXINT) s+=get(p,ni,nj,snj,k,i);
      qm=min(qm,get(q,ni,nj,snj,k,i));
    }
    if (rm+s+qm > tb) tb=rm+s+qm;
  }
  free(snj); free(r); free(q);
  return max(ta,tb);
}

// outputs to CRT the current processing order
void outd(int n, int m, int o, int nj[], int p[], int ni[], int pi[])
{
  int i,j,k;
  int *nk,*snk,*pip;

  nk=(int *)calloc(m+1,MS); snk=(int *)calloc(m+1,MS);

  for (k=1;k<=m;k++) nk[k]=0;

  for (k=1;k<=o;k++) nk[ni[k]]++;
  snk[1]=1; for (k=1;k<m;k++) snk[k+1]=snk[k]+nk[k];

  for (k=1;k<=m;k++)
  {
    pip=pi+snk[k]-1; j=nk[k];
    for(i=1;i<=j;i++) printf("%3d",pip[i]);
    printf("\n");
  }

  free(nk); free(snk);
}

// returns the length of the critical path in the graph G
// finds nodes on the critical path sc[1..lsc]
// existence of the cycle in the graph is NOT checked (graph is surely acyclic)
unsigned path(int n, int a[], int b[], int c[], int p[], int sc[], int *lsc)
{ int i,j,l;
  int ls=0,lp=1,nn;
  unsigned z=0,x;
  int *cp;
  int *stos,*e;
  unsigned *d;

  cp=(int *)calloc(n+1,MS); stos=(int *)calloc(n+1,MS); e=(int *)calloc(n+1,MS); d=(unsigned *)calloc(n+1,MS);

  for (i=1;i<=n;i++) { cp[i]=c[i]; d[i]=0; e[i]=0; if (!cp[i]) stos[++ls]=i; }

  while (lp <= ls)
  { i=stos[lp++]; x=d[i]+p[i]; if (z < x) { z=x; *lsc=i; } j=a[i];
    if (j) { if (!--cp[j]) stos[++ls]=j; if (x > d[j]) { d[j]=x; e[j]=i; }}
    j=b[i];
    if (j) { if (!--cp[j]) stos[++ls]=j; if (x > d[j]) { d[j]=x; e[j]=i; }}
  }
  i=*lsc; j=1; sc[j]=*lsc;
  while (e[i]) { sc[++j]=e[i]; i=e[i]; } nn=j>>1;
  for (i=1;i<=nn;i++) { l=sc[i]; sc[i]=sc[j+1-i]; sc[j+1-i]=l; }
  *lsc=j;

  free(cp); free(stos); free(d); free(e); return z;
}

// returns the makespan for the processing order pi[]
unsigned cmax(int n, int m, int o, int nj[], int p[], int ni[], int pi[])
{ int i,j,k=0;
  unsigned d;
  int *a,*b,*sc;
  int *c;

  a=(int *)calloc(o+1,MS); b=(int *)calloc(o+1,MS); sc=(int *)calloc(o+1,MS); c=(int *)calloc(o+1,MS);

  for (i=1;i<=o;i++) { a[i]=0; b[i]=0; c[i]=0; }

  for (j=1;j<=n;j++) { for (i=k+1;i<k+nj[j];i++) _arc(i,i+1,b,c); k+=nj[j]; }
  for (i=1;i<o;i++) if(ni[pi[i]]==ni[pi[i+1]]) _arc(pi[i],pi[i+1],a,c);

  d=path(o,a,b,c,p,sc,&i); free(a); free(b); free(c); free(sc); return d;
}

// THE MAIN TABU SEARCH ALGORITHM
// Data for the algorithm:
//    n,m,o,nj[0..n],p[1..o],ni[1..o] see INSTANCE DATA,
//    sig[1..o] initial processing order (use simple(), spt() or insa()),
//    maxt    - length of the tabu list,
//    maxl    - length of the list L (for back jumps),
//    maxiter - limit of iterations without improving the makespan,
//    diter   - number of iterations subtracted from maxiter after each backtrack,
//    clen    - length of the cycle detector.
// Results:
//    sig[]   - the best found processing order (!),
//            - the best makespan returns as the value of function taboo().
// For seting TS control variable values see E.Nowicki,C.Smutnicki, 
//        Management Science 6, 1966, 797-813. 
// Other notes:
//    A move is seen basically as a pair of integers (x,y).
//    To simplify technique of operating on moves, the pair of integers (x,y) 
//        is coded as single long variable v. 
//    There are functions pack(), unpx(), unpy() to operate on (x,y) and v.
//    Tabu list tabp[0..maxt-1] has cyclic addressation and stores packed moves;
//        tab[lt] is the youngest element in tabu list (last introduced), so 
//        element next to lt determines the entry point to tab[].
//    List L (long-memory for back jumps) has cyclic addressation [0..maxl] 
//        with each element defined by struct list; element dl is the youngest.

// cyclic forward value of counter k in cycle d
// int pt(int k, int d) { ++k; return k%=d; }          // old
int pt(int k, int d) { return (k+1)%d; }               // modified

// cyclic backward value of counter k in cycle d
// int _pt(int k, int d) { --k; k+=d; return k%=d; }   // old
int _pt(int k, int d) { return (k+d-1)%d; }            // modified

// returns the most recent location of the move v in the tabu list tab[] 
// of lenght maxt with entry point determined by element next to lt;
// returns 0 if move v has not been found in the tabu list tab[]
int intaboo(long v, long tab[], int lt, int maxt)
{ int i=lt,j=0,k;
//  for (k=0;k<maxt;k++) if (!(tab[i=pt(i,maxt)]-v)) j=i+1; return j; }         // old
  for (k=0;k<maxt;k++) { i=(i+1)%maxt; if (!(tab[i]-v)) j=i+1; } return j; }    // modified

// copies "manually" an int array s[0..k-1] to t[0..k-1]
// void copy(int s[], int t[], int k) { int i; for (i=0;i<k;i++) t[i]=s[i]; }    // old
void copy(int *s, int *t, int k) { while(k--) *t++=*s++; }                       // modified

// copies "manually" a long array s[0..k-1] to t[0..k-1]
// void lcopy(long s[], long t[], int k) { int i; for (i=0;i<k;i++) t[i]=s[i]; } // old
void lcopy(long *s, long *t, int k) { while(k--) *t++=*s++; }                    // modified

// copies "manually" an unsigned array s[0..k-1] to t[0..k-1]
// void ucopy(unsigned s[], unsigned t[], int k) { int i; for (i=0;i<k;i++) t[i]=s[i]; }    // old
void ucopy(unsigned *s, unsigned *t, int k) { while(k--) *t++=*s++; }                       // modified

// packs the move (x,y) into single long element
long pack(int x, int y) { return (long)x<<MSX|(long)y; }

// returns the x component from the packed move v
int unpx(long v) { return (int)(v>>MSX); }

// returns the y component from the packed move v
int unpy(long v) { return (int)(v-((v>>MSX)<<MSX)); }

unsigned int taboo(int n, int m, int o, int nj[], int p[], int ni[], int sig[],
               int maxt, int maxiter, int maxl, int diter, int maxd, int maxc)
{
  int i,j,k,l,lt=0,dl=0,ad=0,u,x,y,z,w,lsc,lscp,flag=1,kp,m1,m2,m3,m4,
      lopty=0,lp=maxiter,l1,ipam,flagp,ip,lpc=0,lc=0,dlc=0;
  int *nk,*snk,*snj;
  int *a,*b,*c,*sc,*scp,*cand,*cycle,*_sig,*tabp;
  int *_a,*_b;                      // graph extension, predecessors
  long int *tab,v,vp,iter=0,last;
  unsigned csig,csigr,csigmin,csigopty;
  struct list { unsigned csig; int *sig,kp,*cand,*stos,*e; unsigned *d; long *tab; } *stack;

  int *cp,*stos,*e;
  int *_stos,*ep;                  // inverse stos, auxiliary e
  unsigned *d;
  unsigned *dp;                    // auxiliary d
  unsigned xx;
  int nn,le;
  int maxlen=maxd*maxc;

  int li,lj;

  unsigned llb;

  char ch;
  struct time timep;

// get time for time tracing
  gettime(&timep);

// allocate memory
  nk=(int *)calloc(m+1,MS); snk=(int *)calloc(m+1,MS); snj=(int *)calloc(n+1,MS);
  sc=(int *)calloc(o+1,MS); scp=(int *)calloc(o+1,MS);
  a=(int *)calloc(o+1,MS); b=(int *)calloc(o+1,MS);
  c=(int *)calloc(o+1,MS); _sig=(int *)calloc(o+1,MS);
  _a=(int *)calloc(o+1,MS); _b=(int *)calloc(o+1,MS);
  tab=(long *)calloc(maxt,LS); tabp=(int *)calloc(maxt,MS); 
  cand=(int *)calloc(50,MS); // !!! VALID: this is the very high upper bound 
                       // on the number of moves in the set V 
                       // (usually 5..10 is enough, see our paper)
  cycle=(int *)calloc(maxd,MS);
  stack=(struct list*)calloc(maxl,sizeof(struct list));
  for (i=0;i<maxl;i++)
  { (stack+i)->csig=0; (stack+i)->sig=(int *)calloc(o+1,MS);
    (stack+i)->tab=(long *)calloc(maxt,LS); (stack+i)->kp=0;
    (stack+i)->cand=NULL;

    (stack+i)->stos=(int *)calloc(o+1,MS); 
    (stack+i)->d=(unsigned *)calloc(o+1,MS);
    (stack+i)->e=(int *)calloc(o+1,MS);
  }

  cp=(int *)calloc(o+1,MS); stos=(int *)calloc(o+1,MS); 
  e=(int *)calloc(o+1,MS); d=(unsigned *)calloc(o+1,MS);
  _stos=(int *)calloc(o+1,MS); 
  ep=(int *)calloc(o+1,MS); dp=(unsigned *)calloc(o+1,MS);


  llb=lb(n,m,o,nj,p,ni);  // not necessary for the algorithm 
                    // used only for the trace of relative error to lb

// create auxiliary data structures
  for (i=1;i<=m;i++) nk[i]=0;
  for (i=0;i<maxt;i++) tab[i]=0;
  for (i=1;i<=o;i++) nk[ni[i]]++;
  j=1; for (i=1;i<=m;i++) { snk[i]=j; j+=nk[i]; }            
  j=1; for (i=1;i<=n;i++) { snj[i]=j; j+=nj[i]; }           
  for (i=1;i<=o;i++) _sig[sig[i]]=i;                        // inverse sig

// create starting graph and starting makespan
  for (i=0;i<=o;i++) a[i]=b[i]=c[i]=0; 
  for (j=1;j<=n;j++) for (i=snj[j];i<snj[j]+nj[j]-1;i++) _arc_(i,i+1,b,_b,c);
  for (j=1;j<=m;j++) for (i=snk[j];i<snk[j]+nk[j]-1;i++) _arc_(sig[i],sig[i+1],a,_a,c);

// find starting makespan, starting graph order, paths and the critical path
  { register int i,j,*ls=stos,*lp=stos+1; csig=0;

    for (i=1;i<=o;i++) { d[i]=0; if (!(cp[i]=c[i])) { e[i]=0; *++ls=i; }}
    while (lp <= ls)
    { i=*lp++; xx=d[i]+p[i];
      j=a[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
      j=b[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
             else if (csig < xx) { csig=xx; lsc=i; }
    }
      sc[1]=i=lsc; j=1; while (e[i]) { sc[++j]=e[i]; i=e[i]; } lsc=j;
      nn=lsc>>1; j=lsc+1;
      for (i=1;i<=nn;i++) { le=sc[i]; sc[i]=sc[j-i]; sc[j-i]=le; }
  }
  csigopty=csig;

  for (i=1;i<=o;i++) _stos[stos[i]]=i;                      // inverse stos

start:
// trace calculations
//  if (!(iter%100)) printf(    "              (%8ld)   %6.2f'\n",iter,showtime(timep));

  iter++; lopty++; kp=0;

// analyse critical path and blocks; create candidate moves
// each candidate move is stored in array cand[] as the position in pi[]

  for (i=2;i<=lsc-2;i++)
  { m1=ni[sc[i-1]]; m2=ni[sc[i]]; m3=ni[sc[i+1]]; m4=ni[sc[i+2]];
    if ((m2 == m3)&&((m1 != m2)||(m3 != m4))) cand[++kp]=_sig[sc[i]];
  }
  if ((ni[sc[1]] == ni[sc[2]])&&(ni[sc[2]] != ni[sc[3]]))
    cand[++kp]=_sig[sc[1]];
  if ((ni[sc[lsc-1]] == ni[sc[lsc]])&&(ni[sc[lsc-2]] != ni[sc[lsc-1]]))
    cand[++kp]=_sig[sc[lsc-1]];
  if (!kp) { goto fin; }               // optimal solution has been found

rept:
// select move to be performed
  csigmin=MAXUNSINT; flagp=0; l1=0; for (i=0;i<maxt;tabp[i++]=0);
  for (i=1;i<=kp;i++)
  { w=cand[i]; u=sig[w-1]; x=sig[w]; y=sig[w+1];
//    j=intaboo(pack(sig[cand[i]],sig[cand[i]+1]),tab,lt,maxt);
    j=intaboo(pack(x,y),tab,lt,maxt);
    if (!j)
    { 
//      _swap_(w,sig,a,_a,c); 

// make graph swap
//      a[x]=0; c[y]--;
//      if (a[u]) { c[x]--; a[u]=y; c[y]++; } if (a[y]) a[x]=a[y];
//      a[y]=x; c[x]++;
      if (a[u]&&a[y]) { a[u]=y; a[x]=a[y]; a[y]=x; z=a[y]; _a[z]=x; _a[y]=_a[x]; _a[x]=y; }
      else if (a[u]) { a[u]=y; a[y]=x; a[x]=0; }
        else if (a[y]) { c[y]--; c[x]++; a[x]=a[y]; a[y]=x; z=a[y]; _a[z]=x; _a[y]=_a[x]; _a[x]=y; }
          else { c[y]--; c[x]++; a[y]=x; a[x]=0; }
//      sig[w]=y; sig[w+1]=x;  
      
      csigr=0;

      { register int i,j,*ls=stos,*lp=stos+1;

      for (i=1;i<=o;i++) { d[i]=0; if (!(cp[i]=c[i])) { e[i]=0; *++ls=i; }}
      while (lp <= ls)
      { i=*lp++; xx=d[i]+p[i]; 
        j=a[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
        j=b[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
                 else if (csigr < xx) { csigr=xx; lsc=i; }
      }
//      sc[1]=i=lsc; j=1; while (e[i]) { sc[++j]=e[i]; i=e[i]; } lsc=j;
      }

//      _swap_(w,sig,a,_a,c);

// make inverse graph swap
//      a[y]=0; c[x]--;
//      if (a[u]) { c[y]--; a[u]=x; c[x]++; } if (a[x]) a[y]=a[x];
//      a[x]=y; c[y]++;
      if (a[u]&&a[x]) { a[u]=x; a[y]=a[x]; a[x]=y; z=a[x]; _a[z]=y; _a[x]=_a[y]; _a[y]=x; }
      else if (a[u]) { a[u]=x; a[x]=y; a[y]=0; }
        else if (a[x]) { c[x]--; c[y]++; a[y]=a[x]; a[x]=y; z=a[x]; _a[z]=y; _a[x]=_a[y]; _a[y]=x; }
          else { c[x]--; c[y]++; a[x]=y; a[y]=0; }
//      sig[w]=x; sig[w+1]=y;

      if (csigr < csigmin)
      { csigmin=csigr; ipam=i;
      scp[1]=li=lsc; lj=1; while (e[li]) { scp[++lj]=e[li]; li=e[li]; } // extract critical path only if necessary
      lscp=lj;
//      lscp=lsc; copy(sc+1,scp+1,lsc);   // do not copy
      }
      l1=1; flagp=1;
    }
    else tabp[j-1]=i;
  }
  i=lt; for (k=0;k<maxt;k++)
  { j=tabp[i=pt(i,maxt)];
    if(j)
    { w=cand[j]; u=sig[w-1]; x=sig[w]; y=sig[w+1];
//      _swap_(w,sig,a,_a,c); 

// make graph swap
//      a[x]=0; c[y]--;
//      if (a[u]) { c[x]--; a[u]=y; c[y]++; } if (a[y]) a[x]=a[y];
//      a[y]=x; c[x]++;
      if (a[u]&&a[y]) { a[u]=y; a[x]=a[y]; a[y]=x; z=a[y]; _a[z]=x; _a[y]=_a[x]; _a[x]=y; }
      else if (a[u]) { a[u]=y; a[y]=x; a[x]=0; }
        else if (a[y]) { c[y]--; c[x]++; a[x]=a[y]; a[y]=x; z=a[y]; _a[z]=x; _a[y]=_a[x]; _a[x]=y; }
          else { c[y]--; c[x]++; a[y]=x; a[x]=0; }
//      sig[w]=y; sig[w+1]=x;

      csigr=0;
      { register int i,j,*ls=stos,*lp=stos+1;

      for (i=1;i<=o;i++) { d[i]=0; if (!(cp[i]=c[i])) { e[i]=0; *++ls=i; }}
      while (lp <= ls)
      { i=*lp++; xx=d[i]+p[i]; 
        j=a[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
        j=b[i]; if (j) { if (xx > d[j]) { d[j]=xx; e[j]=i; } if (!--cp[j]) *++ls=j; }
                 else if (csigr < xx) { csigr=xx; lsc=i; }
      }
//      sc[1]=i=lsc; j=1; while (e[i]) { sc[++j]=e[i]; i=e[i]; } lsc=j;
      }

//      _swap_(w,sig,a,_a,c);

// make inverse graph swap
//      a[y]=0; c[x]--;
//      if (a[u]) { c[y]--; a[u]=x; c[x]++; } if (a[x]) a[y]=a[x];
//      a[x]=y; c[y]++;
      if (a[u]&&a[x]) { a[u]=x; a[y]=a[x]; a[x]=y; z=a[x]; _a[z]=y; _a[x]=_a[y]; _a[y]=x; }
      else if (a[u]) { a[u]=x; a[x]=y; a[y]=0; }
        else if (a[x]) { c[x]--; c[y]++; a[y]=a[x]; a[x]=y; z=a[x]; _a[z]=y; _a[x]=_a[y]; _a[y]=x; }
          else { c[x]--; c[y]++; a[x]=y; a[y]=0; }
//      sig[w]=x; sig[w+1]=y;

      if ((csigr < csigmin)&&(csigr < csigopty))
      { csigmin=csigr; flagp=1; l1=1; ipam=j;
      scp[1]=li=lsc; lj=1; while (e[li]) { scp[++lj]=e[li]; li=e[li]; } // extract critical path only if necessary
      lscp=lj;
//      lscp=lsc; copy(sc+1,scp+1,lsc);
      }
      if (!flagp)
      { flagp=1; ipam=j; ip=i;
      scp[1]=li=lsc; lj=1; while (e[li]) { scp[++lj]=e[li]; li=e[li]; } // extract critical path only if necessary
      lscp=lj;
//      lscp=lsc; copy(sc+1,scp+1,lsc); // do not copy path
      if (csigr < csigmin) csigmin=csigr;
      }
    }
  }

// invert scp
   nn=lscp>>1; j=lscp+1;
   for (i=1;i<=nn;i++) { le=scp[i]; scp[i]=scp[j-i]; scp[j-i]=le; }

  if (flag) printf(    " %6u -> %2d (%8ld)   %6.2f'   %6.2f\%\n",csigopty,kp-1,iter,
                     showtime(timep),100.0*(csigopty-llb)/(float)llb);

// update list L if necessary
  if ((flag == 1)&&(kp >= 2))
  { (dl < maxl)?dl++:ad++; k=maxl+ad+dl; k%=maxl;
    (stack+k)->csig=csig; (stack+k)->kp=kp-1;
    copy(sig+1,(stack+k)->sig+1,o); lcopy(tab,(stack+k)->tab,maxt);
    (stack+k)->cand=(int *)calloc(kp,MS); copy(cand+1,(stack+k)->cand+1,ipam-1);
    copy(cand+ipam+1,(stack+k)->cand+ipam,kp-ipam);

    copy(stos+1,(stack+k)->stos+1,o);
    ucopy(d+1,(stack+k)->d+1,o); copy(e+1,(stack+k)->e+1,o);
  }

  if (!l1)
  { last=tab[lt]; j=(maxt+ip-lt)%maxt;
    if (kp!=1) for (i=1;i<=j;i++) tab[lt=pt(lt,maxt)]=last;
  }

// cycle detector
  if (flag) { lpc=dlc=lc=0; for (i=0;i<maxd;i++) cycle[i]=-1; }
  lpc=(++lpc)%maxd; cycle[lpc]=csigmin;
  if (dlc)
  { if (!(cycle[(lpc+dlc)%maxd]-csigmin))
    { if (++lc >= maxlen)
      { printf(    "cycle d [%3d] (%8ld)   %6.2f'\n",dlc,iter,showtime(timep));
      goto backjump; }
      goto cont;
    }
  }
  lc=dlc=0; for (i=1;i<maxd;i++) if (!(cycle[(lpc+i)%maxd]-csigmin)) dlc=i;

cont:
// create the new solution
  flag=0; j=cand[ipam]; tab[lt=pt(lt,maxt)]=pack(sig[j+1],sig[j]);
  swap(j,sig,a,c); _sig[sig[j]]=j; _sig[sig[j+1]]=j+1;
  lsc=lscp; copy(scp+1,sc+1,lsc); csig=csigmin;

// update the best found solution and repeat
  if (csigmin < csigopty) { csigopty=csigmin; flag=1; lopty=0; lp=maxiter; }
  if (lopty <= lp) goto start;

backjump:
// perform back jump; restore the last search region
  lopty=1; flag=1; if (!dl) goto fin;
  k=maxl+ad+dl; k%=maxl; dl--;
  csig=(stack+k)->csig; kp=(stack+k)->kp;
  copy((stack+k)->sig+1,sig+1,o); lcopy((stack+k)->tab,tab,maxt);
  copy((stack+k)->cand+1,cand+1,kp); free((stack+k)->cand); (stack+k)->cand=NULL;

  copy((stack+k)->stos+1,stos+1,o); ucopy((stack+k)->d+1,d+1,o);
  copy((stack+k)->e+1,e+1,o);

  for (i=1;i<=o;i++) _sig[sig[i]]=i;
  for (i=1;i<=o;i++) _stos[stos[i]]=i;
  for (i=0;i<=o;i++) a[i]=b[i]=c[i]=0; 
  for (j=1;j<=n;j++) for (i=snj[j];i<snj[j]+nj[j]-1;i++) _arc_(i,i+1,b,_b,c);
  for (j=1;j<=m;j++) for (i=snk[j];i<snk[j]+nk[j]-1;i++) _arc_(sig[i],sig[i+1],a,_a,c);

  printf("back jump\n");
  printf(    ":%6u ->    (%8ld)   %6.2f'\n",csig,iter,showtime(timep));

  lp=maxiter-(maxl-dl)*diter;        // NOTE !!!: for random instances, after
                             // back jump lp is set in another way
  goto rept;

fin:
// end of calculations
  printf(    "\n %6u -%6u (%8ld) %6.2f'   %6.2f   %c\n",csigopty,llb,iter,
               showtime(timep),100.0*(csigopty-llb)/(float)llb);

  free(nk); free(snk); free(snj);
  free(sc); free(scp); free(a); free(b); free(c); free(_a); free(_b); 
  free(tab); free(tabp); free(cand); free(cycle);
  for (i=0;i<maxl;i++)
  { free((stack+i)->sig); free((stack+i)->tab);
    if ((stack+i)->cand!=NULL) free((stack+i)->cand);

    free((stack+i)->stos); free((stack+i)->d); free((stack+i)->e);
  }
  free(stack); 
  free(cp); free(stos); free(d); free(e);
  free(_stos); free(dp); free(ep);
  return csigopty;
}

void main()
{
  FILE *in;
  int n,m,o,i,j,k,iter,diter,mt,ml,mc,md;
  int *nj,*p,*ni,*pi;
  long t_seed,m_seed;
  char name[12];

//  goto rinst;

// this is the sequence of calls for
// fixed instances, see http://mscmga.ms.ic.ac.uk/info.html, file jobshop1.txt
  clrscr();
  printf("\nTABU SEARCH FOR FIXED INSTANCES\n");
  printf("\nTABU PARAMETERS:\n");
  printf("         maxt: "); scanf("%d",&mt);
  printf("         maxl: "); scanf("%d",&ml);
  printf("      maxiter: "); scanf("%d",&iter);
  printf("        diter: "); scanf("%d",&diter);
  printf("         maxd: "); scanf("%d",&md);
  printf("         maxc: "); scanf("%d",&mc);
  printf("instance file: "); scanf("%s",name);
  
  in=fopen(name,"r");
  fscanf(in,"%d%d",&n,&m); o=n*m;

  nj=(int *)calloc(n+1,MS); p=(int *)calloc(o+1,MS);
  ni=(int *)calloc(o+1,MS); pi=(int *)calloc(o+1,MS);

  for (i=1;i<=n;i++) nj[i]=m;
  for (i=1;i<=o;i++) fscanf(in,"%d%d",?[i],&p[i]);
  for (i=1;i<=o;i++) ni[i]++;
  fclose(in);

//  printf("insa ...\n"); insa(n,m,o,nj,p,ni,pi); printf("   end insa\n");
  printf("insa ...\n"); spt(n,m,o,nj,p,ni,pi); printf("   end insa\n");
  printf("\n\nbest: %5u   lb: %5u\n",taboo(n,m,o,nj,p,ni,pi,mt,iter,ml,diter,md,mc),
          lb(n,m,o,nj,p,ni));

  free(nj); free(p); free(ni); free(pi);
  getch();                                     // temporary stop
//  return;

rinst:
// this is the sequence of calls for
// random instances, see http://mscmga.ms.ic.ac.uk/info.html, file jobshop2.txt
  clrscr();
  printf("\nTABU SEARCH FOR GENERATED INSTANCES\n");
  printf("\nTABU PARAMETERS:\n");
  printf("         maxt: "); scanf("%d",&mt);
  printf("         maxl: "); scanf("%d",&ml);
  printf("      maxiter: "); scanf("%d",&iter);
  printf("        diter: "); scanf("%d",&diter);
  printf("         maxd: "); scanf("%d",&md);
  printf("         maxc: "); scanf("%d",&mc);
  printf("    seed file: "); scanf("%s",name);
  printf("  instance no: "); scanf("%d",&k);

  in=fopen(name,"r");
  fscanf(in,"%d%d",&n,&m); o=n*m;

  nj=(int *)calloc(n+1,MS); p=(int *)calloc(o+1,MS);
  ni=(int *)calloc(o+1,MS); pi=(int *)calloc(o+1,MS);

  for (i=1;i<=n;i++) nj[i]=m;
  for (i=1;i<=k;i++) 
  fscanf(in,"%ld%ld",&t_seed,&m_seed);
  printf("getting seeds:  %10ld   %10ld",t_seed,m_seed);
  generator(&t_seed,&m_seed,n,m,p,ni);
  fclose(in);
  printf("\nlower bound: %5u\n",lb(n,m,o,nj,p,ni));

  printf("insa ...\n"); insa(n,m,o,nj,p,ni,pi); printf("   end insa\n");
  printf("\n\nbest: %5u   lb: %5u\n",taboo(n,m,o,nj,p,ni,pi,mt,iter,ml,diter,md,mc),
          lb(n,m,o,nj,p,ni));

  free(nj); free(p); free(ni); free(pi);
  return;
} 