/* A PACKAGE FOR LOCATING A MINIMUM REGION CONTAINING ALL
   THE ALIGNMENTS WHOSE SCORE IS ABOVE THE SPECIFIED VALUE */

/* Globally passed params and macros */

#include <stdio.h>
#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 */
static long cut;
static int *L, *R;

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

#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 */

/* right_locate - compute the right boundary */
right_locate(I1,J1,I2,J2,d1)
int I1, J1, I2, J2;
long d1;
{        int   midi, midi1, midj, midj1;    /* Midpoint and cost */
	 long  tr, ts;
	 long d2;
	 char  type;

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

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

	if (I2 <= I1) {
		R[I1] = J2;
		return 0;
	}
	if (J2 <= J1) {
		for (i = I1; i <= I2; ++i) R[i] = J2;
		return 0;
	}

	/* Initialize CC, DD, RR, SS, VC, VE */

	for (j = J1; j <= J2; ++j) {
		RR[j] = IRR[j];
		SS[j] = ISS[j];
	}
	for (i = I1; i <= I2; ++i) {
		VC[i] = IVC[i];
		VE[i] = IVE[i];
	}
	DD[J1] = d1;
	CC[J1] = t = VC[I1];
	t = max(VE[I1],t-g);
	for (j = J1+1; j <= J2; ++j) {
		CC[j] = t = t-h;
		DD[j] = t-g;
	}

  	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 = t-m;
		e = c;
		t = t-h;
		wa = w[AA[i+1]];
		for (j = J2-1; j >= J1; --j) BACKWARD
	}
	SS[J2] = t;

	/* Determine the rightmost crossing point */
	midj = -1;
	if ((DD[J1]+SS[J1]-h >= cut) || (CC[J1]+SS[J1]-m >= cut)) {
		midj = J1;
		type = 2;
	}
	wa = w[AA[midi1]];
	for (j = J1+1; j <= J2; ++j) {
		if (CC[j-1]+RR[j]+wa[BB[j]] >= cut) {
			midj = j;
			type = 1;
		}
		if ((DD[j]+SS[j]-h >= cut) || (CC[j]+SS[j]-m >= cut)) {
			midj = j;
			type = 2;
		}
	}
	if (midj == -1) return -1;
	if (type == 1) midj1 = midj-1;
	else midj1 = midj;
	/*
	printf("midj=%d, type= %d\n",midj,type);
	*/
	if (midj > J1) {
		for (j = J1+1; j <= midj; ++j) {
			b = BB[j];
			c = CC[j];
			d = DD[j];
			s = CC[j-1];
			for (i = midi1; i <= I2; ++i) RIGHTWARD
		}
	}
	s = RR[J2];
	RR[J2] = c = SS[J2]-m;
	e = c;
	SS[J2] = SS[J2] - h;
	wa = w[AA[midi1]];
	for (j = J2-1; j >= J1; --j) BACKWARD

	d2 = max(CC[midj]-m, DD[midj]-h);
	tr = IRR[midj];
	ts = ISS[midj];
	for (j = J1; j <= midj1; ++j) {
		IRR[j] = RR[j];
		ISS[j] = SS[j];
	}
	for (i = midi1; i <= I2; ++i) {
		IVC[i] = VC[i];
		IVE[i] = VE[i];
	}
   }
	right_locate(I1,J1,midi,midj1,d1);
	IRR[midj] = tr;
	ISS[midj] = ts;
	right_locate(midi1,midj,I2,J2,d2);
	return 0;
}

/* left_locate - compute the left boundary */
left_locate(I1,J1,I2,J2,s2)
int I1, J1, I2, J2;
long s2;
{        int   midi, midi1, midj, midj1;    /* Midpoint and cost */
	 long  tc, td;
	 long  s1;
	 char  type;

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

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

	if (I2 <= I1) {
	/*
	printf("I1=%d, J1=%d\n",I1,J1);
	*/
		L[I1] = J1;
		return 0;
	}
	if (J2 <= J1) {
	/*
	printf("I1=%d, I2=%d, J1=%d\n",I1,I2,J1);
	*/
		for (i = I1; i <= I2; ++i) L[i] = J1;
		return 0;
	}

	/* Initialize CC, DD, RR, SS, VR, VT */

	SS[J2] = s2;
	RR[J2] = IVR[I2];
	t = IVT[I2]-g;
	for (j = J2-1; j >= J1; --j) {
		SS[j] = RR[j] = t = t-h;
	}
	for (j = J1; j <= J2; ++j) {
		CC[j] = ICC[j];
		DD[j] = IDD[j];
	}
	for (i = I1; i <= I2; ++i) {
		VR[i] = IVR[i];
		VT[i] = IVT[i];
	}

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

	/* Phase 1 : horozontal computation */
	/*           Find the leftmost crossing edge */
	/* Forward computation */
	t = DD[J1];
  	for (i = I1+1; i <= midi; ++i) {
     		s = CC[J1];
      		CC[J1] = c = t = max(s-m,t-h);
      		e = c-g;
      		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;

	/*
	for (j=J1; j <= J2; ++j)
		printf("%d %d %d %d\n",CC[j],DD[j],RR[j],SS[j]);
	*/
	/* Determine the leftmost crossing point */
	midj = -1;
	wa = w[AA[midi1]];
	for (j = J2; j > J1; --j) {
		if ((DD[j]+SS[j]-h >= cut) || (CC[j]+SS[j]-m >= cut)) {
			midj = j;
			type = 2;
		}
		if (CC[j-1]+RR[j]+wa[BB[j]] >= cut) {
			midj = j;
			type = 1;
		}
	}
	if ((DD[J1]+SS[J1]-h >= cut) || (CC[J1]+SS[J1]-m >= cut)) {
		midj = J1;
		type = 2;
	}
	if (midj == -1) return -1;
	if (type == 1) midj1 = midj-1;
	else midj1 = midj;
	if (midj1 < J2) {
		for (j = J2-1; j >= midj1; --j) {
			b = BB[j+1];
			d = SS[j];
			s = RR[j+1];
			for (i = midi; i >= I1; --i) LEFTWARD
		}
	}
     	s = CC[J1];
      	CC[J1] = DD[J1] = c = max(s-m,DD[J1]-h);
      	e = c-g;
      	wa = w[AA[midi1]];
      	for (j = J1+1; j <= J2; ++j) FORWARD

	s1 = max(VR[midi],SS[midj1]-h);
	tc = CC[midj];
	td = DD[midj];
	for (j = midj+1; j <= J2; ++j) {
		ICC[j] = CC[j];
		IDD[j] = DD[j];
	}
	for (i = I1; i <= midi; ++i) {
		IVR[i] = VR[i];
		IVT[i] = VT[i];
	}
   }
	left_locate(I1,J1,midi,midj1,s1);
	ICC[midj] = tc;
	IDD[midj] = td;
	left_locate(midi1,midj,I2,J2,s2);
	return 0;
}

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

	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;
	}
}

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

	IVR[M] = IVT[M] = 0;
	t = -g;
	for (i = M-1; i >= 0; --i) {
		IVR[i] = IVT[i] = t = t-h;
	}
	ICC[0] = 0;
	IDD[0] = t = -g;
	for (j = 1; j <= N; ++j) {
		ICC[j] = t = t-h;
		IDD[j] = t-g;
	}
}

/* Interface and top level of comparator */

long LOCATE(A,B,M,N,W,G,H,LB,RB,CUT_OFF)
char A[],B[]; int M,N; long W[][128],G,H; int LB[],RB[]; long CUT_OFF;

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

  AA = A;			/* Setup global parameters */
  BB = B;
  w = W;
  g = G;
  h = H;
  m = g+h;
  cut = CUT_OFF;
  L = LB;
  R = RB;
  j = (N+1) * sizeof(long);
  CC = (long *) ckalloc(j);
  DD = (long *) ckalloc(j);
  RR = (long *) ckalloc(j);
  SS = (long *) ckalloc(j);
  ICC = IRR = (long *) ckalloc(j);
  IDD = ISS = (long *) ckalloc(j);
  j = (M+1) * sizeof(long);
  VC = VR = (long *) ckalloc(j);
  VE = VT = (long *) ckalloc(j);
  IVC = IVR = (long *) ckalloc(j);
  IVE = IVT = (long *) ckalloc(j);

  right_init(M,N);
  c = right_locate(0,0,M,N,-g);
  if ( c == 0) {
  	left_init(M,N);
  	left_locate(0,0,M,N,0);
  }
  /*
  for (i=0; i<=M; ++i) printf("i=%d, l=%d, r=%d\n",i,L[i],R[i]);
  */
  return c;
}

/* 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);
}