/* A PACKAGE FOR SEQUENCE COMPARISON WITH AFFINE WEIGHTS
   AND ROBUSTNESS MEASURE */

/* Globally passed params and macros */

#include <stdio.h>
#define DIGIT 10.0
#define MININT -999999
#define max(x,y)   ((x) >= (y) ? (x) : (y))
static char *AA, *BB;
static long (*w)[128];				/* w = W */
static long g, h, m;				/* g = G, h = H, m = g+h */

#define gap(k)  ((k) <= 0 ? 0 : g+h*(k))	/* k-symbol indel cost */

static int *sapp;				/* Current script append ptr */
static int  last;				/* Last script op appended */

						/* Append "Delete k" op */
#define DEL(k)				\
{ if (last < 0)				\
    last = sapp[-1] -= (k);		\
  else					\
    last = *sapp++ = -(k);		\
}
						/* Append "Insert k" op */
#define INS(k)				\
{ if (last < 0)				\
    { sapp[-1] = (k); *sapp++ = last; }	\
  else					\
    last = *sapp++ = (k);		\
}

#define REP { last = *sapp++ = 0; }		/* Append "Replace" op */

#define FORWARD						\
{ 	if ((c =   c   - m) > (e =   e   - h)) e = c; 	\
       	if ((c = CC[j] - m) > (d = DD[j] - h)) d = c; 	\
       	c = s + wa[BB[j]];				\
       	if (e > c) c = e;				\
       	if (d > c) c = d;				\
       	s = CC[j];					\
       	CC[j] = c;					\
       	DD[j] = d;					\
}
#define BACKWARD			\
{ 	c = s + wa[BB[j+1]];		\
	d = SS[j]-h;			\
	if (e-m > c) c = e-m;		\
	if (d-g > c) c = d-g;		\
	if (c > d) d = c;		\
	if (c > (e = e-h)) e = c;	\
	s = RR[j];			\
	RR[j] = c;			\
	SS[j] = d;			\
}
#define RIGHTWARD					\
{       if ((c =   c   - m) > (d =   d   - h)) d = c;	\
       	if ((c = VC[i] - m) > (e = VE[i] - h)) e = c;	\
       	c = s + w[AA[i]][b];				\
       	if (d > c) c = d;				\
       	if (e > c) c = e;				\
       	s = VC[i];					\
       	VC[i] = c;					\
       	VE[i] = e;					\
}       
#define LEFTWARD			\
{ 	c = s + w[AA[i+1]][b];		\
	e = VT[i]-h;			\
	if (d-m > c) c = d-m;		\
	if (e-g > c) c = e-g;		\
	if (c > (d = d-h)) d = c;	\
	if (c > e) e = c;		\
	s = VR[i];			\
	VR[i] = c;			\
	VT[i] = e;			\
}
static long *CC, *DD;		/* Forward cost-only vectors */
static long *RR, *SS;		/* Reverse cost-only vectors */
static long *ICC, *IDD;		/* Forward cost-only initial vectors */
static long *IRR, *ISS;		/* Reverse cost-only initial vectors */
static long *VC, *VE;		/* Rightward cost-only vectors */
static long *VR, *VT;		/* Leftward cost-only vectors */
static long *IVC, *IVE;		/* Rightward cost-only initial vectors */
static long *IVR, *IVT;		/* Leftward cost-only initial vectors */
static long *UC, *UR;		/* Robustness measure vectors(col & row) */
static long opt;

/* align(A,B,M,N,tb,te) returns the cost of an optimum conversion between
   A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
   and appends such a conversion to the current script.                   */

static long align(A,B,M,N,tb,te) char *A, *B; int M, N; long tb, te;

{        int   midi, midj, type;	/* Midpoint, type, and cost */
         long midc;
	 long c1, c2;

{ register int   i, j;
  register long c, e, d, s;
           long t, *wa;

/* Boundary cases: M <= 1 or N == 0 */

  if (N <= 0)
    { if (M > 0) DEL(M)
      return -gap(M);
    }
  if (M <= 1)
    { if (M <= 0)
        { INS(N);
          return -gap(N);
        }
      if (tb < te) tb = te;
      midc = (tb-h) - gap(N);
      midj = 0;
      wa = w[A[1]];
      for (j = 1; j <= N; j++)
        { c = -gap(j-1) + wa[B[j]] - gap(N-j);
          if (c > midc)
            { midc = c;
              midj = j;
            }
        }
      if (midj == 0)
        { INS(N) DEL(1) }
      else
        { if (midj > 1) INS(midj-1)
          REP
          if (midj < N) INS(N-midj)
        }
      return midc;
    }

/* Divide: Find optimum midpoint (midi,midj) of cost midc */

  midi = M/2;			/* Forward phase:                          */
  CC[0] = 0;			/*   Compute C(M/2,k) & D(M/2,k) for all k */
  t = -g;
  for (j = 1; j <= N; j++)
    { CC[j] = t = t-h;
      DD[j] = t-g;
    }
  t = tb;
  for (i = 1; i <= midi; i++)
    { s = CC[0];
      CC[0] = c = t = t-h;
      e = t-g;
      wa = w[A[i]];
      for (j = 1; j <= N; j++) 
	{  if ((c =   c   - m) > (e =   e   - h)) e = c; 	
       	   if ((c = CC[j] - m) > (d = DD[j] - h)) d = c; 	
       	   c = s + wa[B[j]];				
       	   if (e > c) c = e;				
       	   if (d > c) c = d;				
       	   s = CC[j];					
       	   CC[j] = c;					
       	   DD[j] = d;			
	}
    }
  DD[0] = CC[0];

  RR[N] = 0;			/* Reverse phase:                          */
  t = -g;			/*   Compute R(M/2,k) & S(M/2,k) for all k */
  for (j = N-1; j >= 0; j--)
    { RR[j] = t = t-h;
      SS[j] = t-g;
    }
  t = te;
  for (i = M-1; i >= midi; i--)
    { s = RR[N];
      RR[N] = c = t = t-h;
      e = t-g;
      wa = w[A[i+1]];
      for (j = N-1; j >= 0; j--)
        { if ((c =   c   - m) > (e =   e   - h)) e = c;
          if ((c = RR[j] - m) > (d = SS[j] - h)) d = c;
          c = s + wa[B[j+1]];
          if (e > c) c = e;
          if (d > c) c = d;
          s = RR[j];
          RR[j] = c;
          SS[j] = d;
        }
    }
  SS[N] = RR[N];

  midc = CC[0]+RR[0];		/* Find optimal midpoint */
  midj = 0;
  type = 1;
  for (j = 0; j <= N; j++)
    if ((c = CC[j] + RR[j]) >= midc)
      if (c > midc || CC[j] != DD[j] && RR[j] == SS[j])
        { midc = c;
          midj = j;
        }
  for (j = N; j >= 0; j--)
    if ((c = DD[j] + SS[j] + g) > midc)
      { midc = c;
        midj = j;
        type = 2;
      }
}

/* Conquer: recursively around midpoint */

  if (type == 1)
    { c1 = align(A,B,midi,midj,tb,-g);
      c2 = align(A+midi,B+midj,M-midi,N-midj,-g,te);
    }
  else
    { align(A,B,midi-1,midj,tb,0);
      DEL(2);
      align(A+midi+1,B+midj,M-midi-1,N-midj,0,te);
    }
  return midc;
}

/* robust - compute robustness measure */
long robust(I1,J1,I2,J2) int I1, J1, I2, J2;
{        int   midi, midi1, midj, midj1;    /* Midpoint and cost */
	 long  tc, td, tr, ts;
         long  opt_num;
	 char  type;

  { register int  i, j, M, N;
    register long c, e, d, s, c1, c2, c3, c4;
	     long ls, ld, rs, rd;
             long t, *wa;
	     char b;

/*
printf("I1=%d, J1=%d, I2=%d, J2=%d\n",I1,J1,I2,J2);
fflush(stdout);
*/

	M = I2 - I1;
	N = J2 - J1;

	/* boundary cases: M <= 0 or N == 0 */
	if (N <= 0) return;
	if (M <= 0) return;

	/* Initialize CC, DD, RR, SS, VC, VE, VR, VS */
	for (j = J1; j <= J2; ++j) {
		CC[j] = ICC[j];
		DD[j] = IDD[j];
		RR[j] = IRR[j];
		SS[j] = ISS[j];
	}
	for (i = I1; i <= I2; ++i) {
		VC[i] = IVC[i];
		VE[i] = IVE[i];
		VR[i] = IVR[i];
		VT[i] = IVT[i];
	}

/* Divide: Find optimum midpoint (midi,midj) of cost midc */

  	midi = (I1+I2)/2;
	midi1 = midi+1;

	/* Phase 1 : horozontal computation */
	/*           Find the rightmost crossing edge */
	/* Forward computation */
	t = DD[J1];
  	for (i = I1+1; i <= midi; ++i) {
     		s = CC[J1];
      		CC[J1] = c = VC[i];
		t = max(s-m, t-h);
      		e = VE[i];
      		wa = w[AA[i]];
      		for (j = J1+1; j <= J2; ++j) FORWARD
	}
	DD[J1] = t;

	/* Backward computation */
	t = SS[J2];
	for ( i = I2-1; i > midi; --i) {
		s = RR[J2];
		RR[J2] = c = VR[i];
		e = VT[i];
		t = max(c, t-h);
		wa = w[AA[i+1]];
		for (j = J2-1; j >= J1; --j) BACKWARD
	}
	SS[J2] = t;

	/* Determine the rightmost crossing point */
	if ((DD[J1]+SS[J1]-h == opt) || (CC[J1]+SS[J1]-m == opt)) {
		midj = J1;
		type = 2;
	}
	wa = w[AA[midi1]];
	for (j = J1+1; j <= J2; ++j) {
		if (CC[j-1]+RR[j]+wa[BB[j]] == opt) {
			midj = j;
			type = 1;
		}
		if ((DD[j]+SS[j]-h == opt) || (CC[j]+SS[j]-m == opt)) {
			midj = j;
			type = 2;
		}
	}
	if (type == 1) midj1 = midj-1;
	else midj1 = midj;
	/*
	printf("midj = %d, type = %d\n",midj,type);
	*/
	/* End of phase 1 */

	/* Phase 2 : determine the largest score outside the subproblems
	   	     for each row and column */
	/* Determine the ... */
	c1 = c2 = MININT;
	if (midj > J1) { /* step 1 */
	   c1 = DD[J1]+SS[J1]-h;
	   if ((t = CC[J1]+SS[J1]-m) > c1) c1 = t;
	   wa = w[AA[midi1]];
	   for (j = J1+1; j < midj; ++j) {
		if ((t = CC[j-1]+RR[j]+wa[BB[j]]) > c1) c1 = t;
		if (c1 > UC[j]) UC[j] = c1;
		if ((t = DD[j]+SS[j]-h) > c1) c1 = t;
		if ((t = CC[j]+SS[j]-m) > c1) c1 = t;
	   }
	   if (type == 2 && (t = CC[midj-1]+RR[midj]+wa[BB[midj]]) > c1) c1 = t;
	   if (c1 > UC[midj]) UC[midj] = c1;
	   /*
	   printf("**** c1 = %d\n",c1);
	   */
	}
	if (midj1 < J2) { /* step 2 */
	   wa = w[AA[midi1]];
	   for (j = J2; j > midj; --j) {
		if ((t = CC[j-1]+RR[j]+wa[BB[j]]) > c2) c2 = t;
		if ((t = DD[j]+SS[j]-h) > c2) c2 = t;
		if ((t = CC[j]+SS[j]-m) > c2) c2 = t;
		if (c2 > UC[j]) UC[j] = c2;
	   }
	   if (type == 1) {
		if ((t = DD[midj]+SS[midj]-h) > c2) c2 = t;
		if ((t = CC[midj]+SS[midj]-m) > c2) c2 = t;
	        if (c2 > UC[midj]) UC[midj] = c2;
	   }
	   /*
	   printf("**** c2 = %d\n",c2);
	   */
	}
	c = max(c1,c2);
	if (c > UR[midi1]) UR[midi1] = c;

	/* Move one row down for whole CC and DD */
	/* Keep initial scores of CC and DD for subproblem 2 unchanged */
	if (midi1 < I2) {
	   t = DD[J1];
     	   s = CC[J1];
      	   CC[J1] = c = VC[midi1];
	   t = max(s-m, t-h);
      	   e = VE[midi1];
      	   wa = w[AA[midi1]];
	   if (midj <= J1)
	   	for (j = J1+1; j <= J2; ++j) FORWARD
	   else {
	   	for (j = J1+1; j < midj; ++j) FORWARD
	   	VC[midi1] = c;
	   	VE[midi1] = e;
	   	for (j = midj; j <= J2; ++j) FORWARD
	   }
	   DD[J1] = t;
	}
	   /*
	   printf("midi1=%d, I2=%d, midj=%d, J1=%d\n",midi1,I2,midj,J1);
	   */
	if (midj > J1) { /* step 3 */
	   /* Move CC[J1..midj-1] and DD[J1..midj-1] down to row I2-1 */
	   t = DD[J1];
  	   for (i = midi1+1; i <= I2-1; ++i) {
     		s = CC[J1];
      		CC[J1] = c = VC[i];
		t = max(s-m, t-h);
      		e = VE[i];
      		wa = w[AA[i]];
      		for (j = J1+1; j < midj; ++j) FORWARD
		VC[i] = c;
		VE[i] = e;
	   }
	   DD[J1] = t;
	   c3 = DD[J1]+ISS[J1]-h;
	   if ((t = CC[J1]+ISS[J1]-m) > c3) c3 = t;
	   wa = w[AA[I2]];
	   for (j = J1+1; j < midj; ++j) {
		if ((t = CC[j-1]+IRR[j]+wa[BB[j]]) > c3) c3 = t;
		if ((t = DD[j]+ISS[j]-h) > c3) c3 = t;
		if ((t = CC[j]+ISS[j]-m) > c3) c3 = t;
	   }
	   if (type == 2 && (t = CC[midj-1]+IRR[midj]+wa[BB[midj]]) > c3) c3 = t;
	   /*
	   printf("**** c3 = %d\n",c3);
	   */
	   /* Move one row down to row I2 for columns J1+1 .. midj-1 */
	   if (midi1 < I2) {
     	      s = CC[J1];
      	      CC[J1] = c = VC[I2];
      	      e = VE[I2];
      	      wa = w[AA[I2]];
      	      for (j = J1+1; j < midj; ++j) FORWARD
	      VC[I2] = c;
	      VE[I2] = e;
	   }
	}
	/* Move one row up for whole RR and SS */
	/* Keep initial scores of RR and SS for subproblem 1 unchanged */
	if (midi > I1) {
	   t = SS[J2];
	   s = RR[J2];
	   RR[J2] = c = VR[midi];
	   e = VT[midi];
	   t = max(c, t-h);
	   wa = w[AA[midi1]];
	   if (midj1 >= J2) {
	   	for (j = J2-1; j >= J1; --j) BACKWARD
	   } else {
		for (j = J2-1; j > midj1; --j) BACKWARD
		VR[midi] = c;
		VT[midi] = e;
		for (j = midj1; j >= J1; --j) BACKWARD
	   }
	   SS[J2] = t;
	}

	/* Move RR[midj1+1..J2] and SS[midj1+1..J2] up to row I1+1 */
	if (midj1 < J2) { /* step 4 */
	   t = SS[J2];
	   for (i = midi-1; i > I1; --i) {
		s = RR[J2];
		RR[J2] = c = VR[i];
		e = VT[i];
		t = max(c, t-h);
		wa = w[AA[i+1]];
		for (j = J2-1; j > midj1; --j) BACKWARD
		VR[i] = c;
		VT[i] = e;
	   }
	   SS[J2] = t;
	   c4 = IDD[J2]+SS[J2]-h;
	   wa = w[AA[I1+1]];
	   for (j = J2; j > midj1; --j) {
		if ((t = ICC[j-1]+RR[j]+wa[BB[j]]) > c4) c4 = t;
		if ((t = IDD[j]+SS[j]-h) > c4) c4 = t;
		if ((t = ICC[j]+SS[j]-m) > c4) c4 = t;
	   }
	   /*
	   printf("**** c4 = %d\n",c4);
	   */
	   /* Move one row up to row I1 for colums midj1+1 .. J2-1 */
	   if (midi > I1) {
		s = RR[J2];
		RR[J2] = c = VR[I1];
		e = VT[I1];
		wa = w[AA[I1+1]];
		for (j = J2-1; j > midj1; --j) BACKWARD
		VR[I1] = c;
		VT[I1] = e;
	   }
	}

	/* Rightward computation */
	if (midj > J1) {
	   t = VE[I1];
	   for (j = J1+1; j < midj; ++j) {
		b = BB[j];
		c = ICC[j];
		d = IDD[j];
		s = ICC[j-1];
		t = max(s-m, t-h);
      		for (i = I1+1; i <= midi; ++i) RIGHTWARD
	   }
	   VC[I1] = ICC[midj-1];
	   rs = VC[midi];
	   if (type == 2) {
		b = BB[midj];
		c = ICC[midj];
		d = IDD[midj];
		s = ICC[midj-1];
		t = max(s-m, t-h);
		for (i = I1+1; i <= midi; ++i) RIGHTWARD
		VC[I1] = ICC[midj];
		rd = d;
	   }
	   VE[I1] = t;
	}

	/* Leftward computation */
	if (midj1 < J2) {
	   if (type == 1) {
	      t = VT[I2];
	      for (j = J2-1; j >= midj; --j) {
		   b = BB[j+1];
		   d = ISS[j];
		   s = IRR[j+1];
		   t = max(IRR[j],t-h);
		   for (i = I2-1; i > midi; --i) LEFTWARD
	      }
	   } else {
	      t = VT[I2];
	      for (j = J2-1; j > midj; --j) {
		   b = BB[j+1];
		   d = ISS[j];
		   s = IRR[j+1];
		   t = max(IRR[j],t-h);
		   for (i = I2-1; i > midi; --i) LEFTWARD
	      }
	      ls = VR[midi1];
	      b = BB[midj+1];
	      d = ISS[midj];
	      s = IRR[midj+1];
	      t = max(IRR[midj],t-h);
	      for (i = I2-1; i > midi; --i) LEFTWARD
	      ld = d;
	   }
	   VT[I2] = t;
	   VR[I2] = IRR[midj];
	}

	/* steps 5 & 6 */
	if (midj > J1) { /* step 5 */
	   b = BB[midj];
	   for (i = I2; i > midi1; --i) {
		if ((t = VE[i]+VT[i]-h) > c3) c3 = t;
		if ((t = VC[i]+VT[i]-m) > c3) c3 = t;
		if ((t = VC[i-1]+VR[i]+w[AA[i]][b]) > c3) c3 = t;
		if (c3 > UR[i]) UR[i] = c3;
	   }
	   /*
	   printf("**** c5 = %d\n",c3);
	   */
	}
	if (midj1 < J2) { /* step 6 */
	   b = BB[midj1+1];
	   for (i = I1; i < midi; ++i) {
		if ((t = VE[i]+VT[i]-h) > c4) c4 = t;
		if ((t = VC[i]+VT[i]-m) > c4) c4 = t;
		if ((t = VC[i]+VR[i+1]+w[AA[i+1]][b]) > c4) c4 = t;
		if (c4 > UR[i+1]) UR[i+1] = c4;
	   }
	   /*
	   printf("**** c6 = %d\n",c4);
	   */
	}

	/* Compute columewise initial scores */
	/* Rightward */
	if (midj > J1 && midj < J2) {
	   if (type == 1) {
	  	b = BB[midj];
	   	c = ICC[midj];
	   	d = IDD[midj];
	   	s = ICC[midj-1];
	   	VE[I1] = max(s-m, VE[I1]-h);
	   	for (i = I1+1; i <= I2; ++i) RIGHTWARD
	   	VC[I1] = ICC[midj];
	   } else {
		b = BB[midj];
		c = VC[midi];
		d = rd;
		s = rs;
		for (i = midi1; i <= I2; ++i) RIGHTWARD
	   }
	}
	/* Leftward */
	if (midj1 < J2 && midj1 > J1) {
	   if (type == 1) {
	      b = BB[midj];
	      d = ISS[midj1];
	      s = IRR[midj];
	      VT[I2] = max(IRR[midj1],VT[I2]-h);
	      for (i = I2-1; i >= I1; --i) LEFTWARD
	      VR[I2] = IRR[midj1];
	   } else {
	      b = BB[midj+1];
	      d = ld;
	      s = ls;
	      for (i = midi; i >= I1; --i) LEFTWARD;
	   }
	}
	
	/* Steps 7 & 8 */
	/* Leftward */
	if (midj > J1) { /* step 7 */
	   t = VT[I2];
	   for (j = midj1-1; j > J1; --j) {
	   	b = BB[j+1];
	   	d = ISS[j];
	   	s = IRR[j+1];
	   	t = max(IRR[j],t-h);
	   	for (i = I2-1; i >= midi1; --i) LEFTWARD
       	   }
	   VT[I2] = t;
	   VR[I2] = IRR[J1+1];

	   b = BB[J1+1];
	   c = IVE[midi1]+VT[midi1]-h;
	   if ((t = IVC[midi1]+VT[midi1]-m) > c) c = t;
	   for (i = I2; i > midi1; --i) {
		if ((t = IVE[i]+VT[i]-h) > c) c = t;
		if ((t = IVC[i]+VT[i]-m) > c) c = t;
		if ((t = IVC[i-1]+VR[i]+w[AA[i]][b]) > c) c = t;
	   }
	   /*
	   printf("**** c7 = %d\n",c);
	   */
	   for (j = J1+1; j <= midj;  ++j)
		if (c > UC[j]) UC[j] = c;
	}
	/* Rightward */
	if (midj1 < J2) { /* step 8 */
	   t = VE[I1];
	   for (j = midj+1; j < J2; ++j) {
		b = BB[j];
		c = ICC[j];
		d = IDD[j];
		s = ICC[j-1];
		t = max(s-m, t-h);
      		for (i = I1+1; i <= midi; ++i) RIGHTWARD
	   }
	   VE[I1] = t;
	   VC[I1] = ICC[J2-1];

	   b = BB[J2];
	   c = VE[I1]+IVT[I1]-h;
	   if ((t = VC[I1]+IVT[I1]-m) > c) c = t;
	   for (i = I1+1; i <= midi; ++i) {
		if ((t = VE[i]+IVT[i]-h) > c) c = t;
		if ((t = VC[i]+IVT[i]-m) > c) c = t;
		if ((t = VC[i-1]+IVR[i]+w[AA[i]][b]) > c) c = t;
	   }
	   /*
	   printf("**** c8 = %d\n",c);
	   */
	   for (j = J2; j > midj1; --j )
		if (c > UC[j]) UC[j] = c;
	}
/*
printf("CC: %d %d %d %d\n",CC[J1],ICC[J1],CC[J2],ICC[J2]);
printf("DD: %d %d %d %d\n",DD[J1],IDD[J1],DD[J2],IDD[J2]);
printf("RR: %d %d %d %d\n",RR[J1],IRR[J1],RR[J2],IRR[J2]);
printf("SS: %d %d %d %d\n",SS[J1],ISS[J1],SS[J2],ISS[J2]);
printf("VC: %d %d %d %d\n",VC[I1],IVC[I1],VC[I2],IVC[I2]);
printf("VE: %d %d %d %d\n",VE[I1],IVE[I1],VE[I2],IVE[I2]);
printf("VR: %d %d %d %d\n",VR[I1],IVR[I1],VR[I2],IVR[I2]);
printf("VT: %d %d %d %d\n",VT[I1],IVT[I1],VT[I2],IVT[I2]);
*/
	/* Setup ICC, IDD, IRR, ISS */
	tc = CC[midj];
	td = DD[midj];
	tr = IRR[midj];
	ts = ISS[midj];
	for (j = J1; j <= midj; ++j) { /* Subproblem 1 */
		IRR[j] = RR[j];
		ISS[j] = SS[j];
	}
	for (j = midj+1; j <= J2; ++j) { /* Subproblem 2 */
		ICC[j] = CC[j];
		IDD[j] = DD[j];
	}
	/* Setup IVC, IVE, IVR, IVT */
	for (i = I1; i <= midi; ++i) { /* Subproblem 1 */
		IVR[i] = VR[i];
		IVT[i] = VT[i];
	}
	for (i = midi1; i <= I2; ++i) { /* Subproblem 2 */
		IVC[i] = VC[i];
		IVE[i] = VE[i];
	}
  }

  /* Phase 3: recursively solve subproblems */
  robust(I1,J1,midi,midj1);
  ICC[midj] = tc;
  IDD[midj] = td;
  IRR[midj] = tr;
  ISS[midj] = ts;
  robust(midi1,midj,I2,J2);
}

/* robust_init - initializations for robust() */
robust_init(M,N) int M, N;
{
	register int i, j;
	long t;

	ICC[0] = 0;
	IDD[0] = t = -g;
	for (j = 1; j <= N; ++j) {
		ICC[j] = t = t-h;
		IDD[j] = t-g;
	}
	IVC[0] = 0;
	IVE[0] = t = -g;
	for (i = 1; i <= M; ++i) {
		IVC[i] = t = t-h;
		IVE[i] = t-g;
	}
	IRR[N] = ISS[N] = 0;
	t = -g;
	for (j = N-1; j >= 0; --j) {
		IRR[j] = ISS[j] = t = t-h;
	}
	IVR[M] = IVT[M] = 0;
	t = -g;
	for (i = M-1; i >= 0; --i) {
		IVR[i] = IVT[i] = t = t-h;
	}
	t = -gap(M)-gap(N);
	for (i = 1; i <= M; ++i) UR[i] = t;
	for (j = 1; j <= N; ++j) UC[j] = t;
}

/* Interface and top level of comparator */

long ALIGN(A,B,M,N,W,G,H,S) char A[],B[]; int M,N; long W[][128],G,H; int S[];

{ 
  int i, j;
  long ck;
  char *ckalloc();


  AA = A;			/* Setup global parameters */
  BB = B;
  w = W;
  g = G;
  h = H;
  m = g+h;
  sapp = S;
  last = 0;
  j = (N+1) * sizeof(long);
  CC = (long *) ckalloc(j);
  DD = (long *) ckalloc(j);
  RR = (long *) ckalloc(j);
  SS = (long *) ckalloc(j);
  ICC = (long *) ckalloc(j);
  IDD = (long *) ckalloc(j);
  IRR = (long *) ckalloc(j);
  ISS = (long *) ckalloc(j);
  UC = (long *) ckalloc(j);
  j = (M+1) * sizeof(long);
  VC = (long *) ckalloc(j);
  VE = (long *) ckalloc(j);
  VR = (long *) ckalloc(j);
  VT = (long *) ckalloc(j);
  IVC = (long *) ckalloc(j);
  IVE = (long *) ckalloc(j);
  IVR = (long *) ckalloc(j);
  IVT = (long *) ckalloc(j);
  UR = (long *) ckalloc(j);

  opt = align(A,B,M,N,-g,-g);	/* OK, do it */
  ck = CHECK_SCORE(A,B,M,N,S);
  if (opt != ck) printf("Check_score error.\n");
  robust_init(M,N);
  robust(0,0,M,N);
  return opt;
}

/* Alignment display routine */

static char ALINE[51], BLINE[51], CLINE[51];

long DISPLAY(A,B,M,N,S,AP,BP) char A[], B[]; long M, N; long S[], AP, BP;
{ register char *a, *b, *c;
  register long   i,  j, op;
           long   lines, ap, bp;

  i = j = op = lines = 0;
  ap = AP;
  bp = BP;
  a = ALINE;
  b = BLINE;
  c = CLINE;
  while (i < M || j < N)
    { if (op == 0 && *S == 0)
        { op = *S++;
          *a = A[++i];
          *b = B[++j];
	  if (UR[i] == opt)
          	*c++ = (*a++ == *b++) ? '|' : ' ';
	  else
          	*c++ = (*a++ == *b++) ? '*' : '!';
        }
      else
        { if (op == 0)
            op = *S++;
          if (op > 0)
            { *a++ = ' ';
              *b++ = B[++j];
	      *c++ = (UC[j] == opt) ? '-' : '=';
              op--;
            }
          else
            { *a++ = A[++i];
              *b++ = ' ';
	      *c++ = (UR[i] == opt) ? '-' : '=';
              op++;
            }
        }
      if (a >= ALINE+50 || i >= M && j >= N)
        { *a = *b = *c = '\0';
          printf("\n%5d ",50*lines++);
          for (b = ALINE+10; b <= a; b += 10)
            printf("    .    :");
          if (b <= a+5)
            printf("    .");
          printf("\n%5d %s\n      %s\n%5d %s\n",ap,ALINE,CLINE,bp,BLINE);
	  ap = AP + i;
	  bp = BP + j;
          a = ALINE;
          b = BLINE;
          c = CLINE;
        }
    }
}

/* CHECK_SCORE - return the score of the alignment stored in S */

static long CHECK_SCORE(A,B,M,N,S) char A[], B[]; long M, N; long S[];
{ 
  register long   i,  j, op;
  long score;

  score = i = j = op = 0;
  while (i < M || j < N) {
	op = *S++;
	if (op == 0) 
		score = w[A[++i]][B[++j]] + score;
	else if (op > 0) {
		score = score - (g+op*h);
		j = j+op;
	} else {
		score = score - (g-op*h);
		i = i-op;
	}
  }
  return(score);
}

/* ROBUST - return robustness measure */

long ROBUST(A,B,M,N,S,OPT) char A[], B[]; long M, N; long S[]; long OPT;
{ 
  register long   i, j, k, op;
  long score, rc, ac;

  printf("\n\n");
  ac = rc = score = i = j = op = 0;
  while (i < M || j < N) {
	op = *S++;
	if (op == 0) {
		ac++;
		++i;
		++j;
		if (OPT-UR[i] != 0) rc++;
		if (UR[i] != UC[j]) {
			score = ralign(A,B,M,N,i,j,0);
			printf("Spair(%6d,%6d): robustness measure= %.1f,   ",
				i,j,((float)(OPT-score))/DIGIT);
			printf("rm = %.1f,  ",((float)(OPT-UR[i]))/DIGIT);
			printf("%.1f\n",((float)(OPT-UC[j]))/DIGIT);
		}
		/*
		score = ralign(A,B,M,N,i,j,0);
		if (OPT - score != 0) {
			rc++;
			printf("Spair(%6d,%6d): robustness measure= %.1f,   ",
				i,j,((float)(OPT-score))/DIGIT);
			printf("rm = %.1f,  ",((float)(OPT-UR[i]))/DIGIT);
			printf("%.1f\n",((float)(OPT-UC[j]))/DIGIT);
		} else {
			printf("Spair(%6d,%6d): not robust,   ",i,j);
			printf("rm = %.1f,  ",((float)(OPT-UR[i]))/DIGIT);
			printf("%.1f\n",((float)(OPT-UC[j]))/DIGIT);
		}
		*/
	} else if (op > 0) {
		for (k = 1; k <= op; ++k) {
		   ac++;
		   ++j;
		   if (UC[j] != OPT) rc++;
		   /*
		   score = ralign(A,B,M,N,i,j,2);
		   if (OPT - score != 0) {
			rc++;
			printf("Ipair(%6d,%6d): robustness measure= %.1f,   ",
				i,j,((float)(OPT-score))/DIGIT);
			printf("rm = %.1f\n",((float)(OPT-UC[j]))/DIGIT);
		   } else {
			printf("Ipair(%6d,%6d): not robust,   ", i,j);
			printf("rm = %.1f\n",((float)(OPT-UC[j]))/DIGIT);
		   }
		   */
		}
	} else {
		for (k = 1; k <= -op; ++k) {
		   ac++;
		   ++i;
		   if (UR[i] != OPT) rc++;
		   /*
		   score = ralign(A,B,M,N,i,j,1);
		   if (OPT - score != 0) {
			rc++;
			printf("Dpair(%6d,%6d): robustness measure= %.1f,   ",
				i,j,((float)(OPT-score))/DIGIT);
			printf("rm = %.1f\n",((float)(OPT-UR[i]))/DIGIT);
		   } else {
			printf("Dpair(%6d,%6d): not robust,   ", i,j);
			printf("rm = %.1f\n",((float)(OPT-UR[i]))/DIGIT);
		   }
		   */
		}
	}
  }
  if (ac == 0)
     printf("No diagonal aligned pairs\n"); 
  else {
     printf("Percentage of robust pairs is %3d %% \n",(int)(((float)rc/(float)ac)*100.0));
  }
}

static long ralign(A,B,M,N,I,J,type) char *A, *B; int M, N, I, J; char type;

{
  register int   i, j;
  register long c, e, d, s;
           long t, *wa;

  CC[0] = 0;
  t = -g;
  for (j = 1; j <= N; j++)
    { 
      if (type == 2 && I == 0 && J == j) t = MININT;
      CC[j] = t = t-h;
      DD[j] = t-g;
    }
  t = -g;
  for (i = 1; i <= M; i++)
    { s = CC[0];
      if (type == 1 && J == 0 && I == i) t = MININT;
      CC[0] = c = t = t-h;
      e = t-g;
      wa = w[A[i]];
      for (j = 1; j <= N; j++)
        { if (type == 2 && i == I && j == J) c = e = MININT;
	  else if ((c =   c   - m) > (e =   e   - h)) e = c;
	  if (type == 1 && i == I && j == J) d = MININT;
          else if ((c = CC[j] - m) > (d = DD[j] - h)) d = c;
	  if (type == 0 && i == I && j == J) c = e;
          else c = s + wa[B[j]];
          if (e > c) c = e;
          if (d > c) c = d;
          s = CC[j];
          CC[j] = c;
          DD[j] = d;
	}
    }
  return CC[N];
}

/* lib.c - library of C procedures. */

/* fatal - print message and die */
fatal(msg)
char *msg;
{
	fprintf(stderr, "%s\n", msg);
	exit(1);
}

/* fatalf - format message, print it, and die */
fatalf(msg, val)
char *msg, *val;
{
	fprintf(stderr, msg, val);
	putc('\n', stderr);
	exit(1);
}
	
/* ckopen - open file; check for success */
FILE *ckopen(name, mode)
char *name, *mode;
{
	FILE *fopen(), *fp;

	if ((fp = fopen(name, mode)) == NULL)
		fatalf("Cannot open %s.", name);
	return(fp);
}

/* ckalloc - allocate space; check for success */
char *ckalloc(amount)
long amount;
{
	char *malloc(), *p;

	if ((p = malloc( (unsigned) amount)) == NULL)
		fatal("Ran out of memory.");
	return(p);
}