/*
* sim3.c - compare two very similar sequences with affine gap penalties
*
* Version 1 : July 16, 1995.
* Version 2 : Aug. 31, 1995.  (Heuristic phase 1 + New output format)
* Version 2.1 : Nov. 4, 1995. (Rewrite sim3 as a procedure call)
*
* The command syntax is
*	sim3 file1 file2 [dist_limit]
* where file1 and file2 contain arbitrary sequences of characters.  It is
* assumed that the two character strings are very similar.  Sim3 reports
* an alignment of the two sequences.
*
* If dist_limit > 0, a heuristic phase 1 is used.
* If dist_limit <= 0, an exact phase 1 is used.
* Absolute(dist_limit) specifies the tolerable distance limit.
* (Zero means a default distance limit is used.)
*
* For example,
*  sim3 file1 file2      (A heuristic phase 1 is used with default limit.)
*  sim3 file1 file2 500  (A heuristic phase 1 is used with limit=500.)
*  sim3 file1 file2 0    (An exact phase 1 is used with default limit.)
*  sim3 file1 file2 -500 (An exact phase 1 is used with limit=500.)
*/

#include <stdio.h>
#include <ctype.h>

#define max(x,y)  ((x) >= (y) ? (x) : (y))
#define min(x,y)  ((x) <= (y) ? (x) : (y))
#define FALSE 0
#define TRUE 1
#define SUB 0
#define INS 1
#define DEL 2
#define MININT -99999
#define DEF_LIMIT 1000
#define W	8	/* Word size */
#define NW	60	/* Number of words used in the heuristic method */

char *A, *B;
int M, N;
int limit;

/* Edit_script is the internal representation of an alignment. */

typedef struct edit_script {
	char op_type;	/* SUB, INS, or DEL */
	int num;	/* Number of operations */
	struct edit_script *next;
} edit_script;

int get_dist(int, int, int, int);
char locate(int *, int *, int *, int *);
int dist_only(int *, int *, int *, int *);
int snake(int, int, int, int);
int rsnake(int, int, int, int);
void path(int, int, char, int, int, char, int, edit_script **, edit_script **);
void Condense_script(edit_script *);
int Check_script(edit_script *, int, int);
void Reverse_Script(edit_script *);
void Print_Script(edit_script *, int, int);
void S2A(edit_script *, int *);
void DISPLAY(char *, char *, int, int, int *, int, int);
void Script_free(edit_script *);

extern void fatal(char *);
extern void fatalf(char *, char *);
extern FILE *ckopen(char *, char *);
extern char *ckalloc(int);
static int power(int, int);

/* sim3 - a procedure for comparing two very similar sequences */
/* INPUT: sequence1, length1,
	  sequence2, length2,
	  distance limit (to bound the computation)

   OUTPUT: 1. the distance between sequence1 and sequence2;
	   2. an edit script describing the alignment.
*/
edit_script *sim3(char *seq1, int len1, char *seq2, int len2, int LIMIT, int *d)
{
	int dist, check_dist;
	char sflag, hflag;
	char *C;
	int t;
	edit_script *head, *tail, *tp;
	int *S, i;
	int I1, J1, I2, J2;

	A = seq1;
	M = len1;
	B = seq2;
	N = len2;
	limit = LIMIT;

	sflag = FALSE;
	if (M>N) {	/* Switch the sequences */
		C = A;
		A = B;
		B = C;
		t = M;
		M = N;
		N = t;
		sflag = TRUE;
	}

	/* Compute the distance between two sequences A and B */

	if (limit > 0) {
                printf("A heuristic phase 1 is tried.\n");
                hflag = locate(&I1, &J1, &I2, &J2);
                if (hflag == TRUE) dist = get_dist(I1, J1, I2, J2);
                else {
			printf("An exact phase 1 is used.\n");
			dist = dist_only(&I1, &J1, &I2, &J2);
		}
	} else {
		printf("An exact phase 1 is used.\n");
		if (limit == 0) limit = DEF_LIMIT;
		else limit = -limit;
		dist = dist_only(&I1, &J1, &I2, &J2);
	}

	printf("distance limit = %d\n", limit);

#ifdef STATS
	printf("I1=%d, J1=%d, I2=%d, J2=%d, dist=%d\n",I1, J1, I2, J2, dist);
#endif

	/* Deliver the actual alignment.
	   Head is a pointer to an alignment.
	   See the description of the edit_script structure. */

	path(I1, J1, SUB, I2, J2, SUB, dist, &head, &tail);

	/* Merge the contiguous same operations together */
	Condense_script(head);

	/* Reset tail (due to the condensation) */
	tail = head;
	while (tail->next != NULL) tail = tail->next;

#ifdef STATS
	check_dist = Check_script(head, I1, J1);
	printf("check_dist=%d\n", check_dist);
#endif

	/* Add end gaps, if necessary */
	if (J1 > 0) {
	    if (head->op_type == INS) head->num += J1;
	    else {
		tp = (edit_script *) ckalloc(sizeof(edit_script));
		tp->op_type = INS;
		tp->num = J1;
		tp->next = head;
		head = tp;
/*
printf("left end_gaps: num=%d\n",tp->num);
*/
	    }
	}

	if (J2 < N) {
	    if (tail->op_type == INS) tail->num = tail->num + N - J2;
	    else {
/*
printf("tail->op_type=%d, tail->num=%d\n", tail->op_type, tail->num);
*/
                tp = (edit_script *) ckalloc(sizeof(edit_script));
                tp->op_type = INS;
                tp->num = N-J2;
		tp->next = NULL;
                tail->next = tp;
                tail = tp;
/*
printf("right end_gaps: num=%d\n",tp->num);
*/
	    }
	}

	if (sflag == TRUE) {	/* Switch back the sequences */
		C = A;
		A = B;
		B = C;
		t = M;
		M = N;
		N = t;
		Reverse_Script(head);
	}

	/* Transform the edit_script to the format
	   used in DISPLAY routine. */

#ifdef STATS
	S = (int *) ckalloc((M+N)*sizeof(int));
	S2A(head, S);
	DISPLAY(A-1,B-1,M,N,S,1,1);
	free(S);
#endif

	*d = dist;
	return head;
}

char locate(int *i1, int *j1, int *i2, int *j2)
{
	int encoding[128];
	int mask;
	int BUCKETS;
	char *bucket;
	int i, j;
	char *p, *q;
	long ecode, ecode1;
	int fcount, bcount, fmax, bmax;
	float trust_ratio;
	int paired_fj, paired_bj, best_range, t;
	char pflag;

	typedef struct max_chain {
		int j;	/* j position */
		struct max_chain *next;
	} max_chain;
	max_chain *fj, *bj, *tj, *tj1, *tj2;

	if (M < W+NW-1) {
		printf("The sequences are too short for the heuristic method.\n");
		return FALSE;
	}

	/* Set up the hash table (derived from Sequence A) */
	for (i=0; i<128; ++i) encoding[i] = 0;
	encoding['C'] = encoding['c'] = 1;
	encoding['G'] = encoding['g'] = 2;
	encoding['T'] = encoding['t'] = 3;
	mask = (1 <<(W+W-2)) - 1;
	BUCKETS = power(4, W);
	bucket = (char *) ckalloc(BUCKETS * sizeof(char));
	for (i=0; i<BUCKETS; ++i) bucket[i] = 0;

	/* front end */
	p = A;
	ecode = 0L;
	for (i=1; i<W; ++i)
		ecode = (ecode << 2) + encoding[*p++];
	for (i=0; i<NW; ++i) {
		ecode = ((ecode & mask) << 2) + encoding[*p++];
		bucket[ecode] = 1;
	}

	/* back end */
	p = A+M-W-NW+1;
	ecode = 0L;
        for (i=1; i<W; ++i)
                ecode = (ecode << 2) + encoding[*p++];
        for (i=0; i<NW; ++i) {
                ecode = ((ecode & mask) << 2) + encoding[*p++];
                if (bucket[ecode] == 1 || bucket[ecode] == 3) bucket[ecode] = 3;
		else bucket[ecode] = 2;
        }

	/* Locate the interval with the maximum number of the matched 8-tuples */
	trust_ratio = 0.6;
	fmax = bmax = NW * trust_ratio;
	fj = bj = NULL;
	fcount = bcount = 0;
	p = B;
	ecode = 0L;
	for (j=1; j<W; ++j)
                ecode = (ecode << 2) + encoding[*p++];
	ecode1 = ecode;
	q = p;
	for (j=0; j<NW; ++j) {
                ecode = ((ecode & mask) << 2) + encoding[*p++];
		if (bucket[ecode] == 1 || bucket[ecode] == 3) ++fcount;
		if (bucket[ecode] == 2 || bucket[ecode] == 3) ++bcount;
	}

	if (fcount >= fmax) {
		fj = (max_chain *) ckalloc(sizeof(max_chain));
		fj->j = 0;
		fj->next = NULL;
		fmax = fcount;
	}
	if (bcount >= bmax) {
		bj = (max_chain *) ckalloc(sizeof(max_chain));
		bj->j = W+NW-2;
		bj->next = NULL;
		bmax = bcount;
	}

	for (j=1; j <= N-W-NW+1; ++j) {
		ecode1 = ((ecode1 & mask) << 2) + encoding[*q++];
		if (bucket[ecode1] == 1 || bucket[ecode1] == 3) --fcount;
                if (bucket[ecode1] == 2 || bucket[ecode1] == 3) --bcount;
                ecode = ((ecode & mask) << 2) + encoding[*p++];
		if (bucket[ecode] == 1 || bucket[ecode] == 3) ++fcount;
                if (bucket[ecode] == 2 || bucket[ecode] == 3) ++bcount;

		if (fcount > fmax) {
			while (fj != NULL) {
				tj = fj->next;
				free(fj);
				fj = tj;
			}
			fj = (max_chain *) ckalloc(sizeof(max_chain));
			fj->j = j;
			fj->next = NULL;
			fmax = fcount;
		} else if (fcount == fmax) {
			tj = (max_chain *) ckalloc(sizeof(max_chain));
			tj->j = j;
			tj->next = fj;
			fj = tj;
		}

		if (bcount > bmax) {
			while (bj != NULL) {
				tj = bj->next;
				free(bj);
				bj = tj;
			}
			bj = (max_chain *) ckalloc(sizeof(max_chain));
			bj->j = j+W+NW-2;
			bj->next = NULL;
			bmax = bcount;
		} else if (bcount == bmax) {
			tj = (max_chain *) ckalloc(sizeof(max_chain));
			tj->j = j+W+NW-2;
			tj->next = bj;
			bj = tj;
		}
	}

	free(bucket);
	
#ifdef STATS
	printf("NW=%d, fmax=%d, bmax=%d\n", NW, fmax, bmax);
#endif

	if (fj == NULL || bj == NULL) {
		printf("The end-intervals are too different for the heuristic method.\n");
		return FALSE;
	}

	best_range = limit+1;
	for (tj1 = fj; tj1 != NULL; tj1 = tj1->next) {
		for (tj2 = bj; tj2 != NULL; tj2 = tj2->next) {
			t = tj2->j - tj1->j + 1 - M;
			if (t<0) t = -t;
/*
printf("tj1=%d, tj2=%d, t=%d\n", tj1->j, tj2->j, t);
*/
			if (t<best_range) {
				best_range = t;
				paired_fj = tj1->j;
				paired_bj = tj2->j;
			}
		}
	}

#ifdef STATS
	printf("paired_fj=%d, paired_bj=%d, best_range=%d\n", paired_fj, paired_bj, best_range);
#endif

	if (best_range > limit) {
		printf("Heuristic phase 1 failed\n");
		return FALSE;
	}

	*i1 = 0;
	*i2 = M;
	*j1 = paired_fj;
	*j2 = paired_bj+1;
	return TRUE;

}

static int power(int base, int n)
{
	int i, p;

	for (i=1, p=1; i<=n; ++i) p = p*base;
	return p;
}

int get_dist(int i1, int j1, int i2, int j2)
{
	int *SS, *DD, *II;
	int goal_diag;
	int c, k, t1, t2, t;
	int start, lower, upper;

	/* Compute the boundary diagonals */
        start = j1 - i1;
        lower = max(j1-i2, start-limit);
        upper = min(j2-i1, start+limit);
	goal_diag = j2-i2;

	if (goal_diag > upper || goal_diag < lower) {
		printf("Two sequences are not really similar\n.");
		printf("Please try exact phase 1 method\n.");
		exit(1);
	}

        /* Allocate space for forward vectors */
        SS = (int *)ckalloc((upper-lower+1)*sizeof(int)) - lower;
        DD = (int *)ckalloc((upper-lower+2)*sizeof(int)) - lower;
        II = (int *)ckalloc((upper-lower+2)*sizeof(int)) - lower + 1;

	/* Initialization */
	for (k=lower; k<=upper; ++k) SS[k] = MININT;
        for (k=lower; k<=upper+1; ++k) DD[k] = MININT;
        for (k=lower-1; k<=upper; ++k) II[k] = MININT;
	SS[start] = snake(start, i1, i2, j2);

	if (SS[goal_diag] >= i2) {
#ifdef STATS
		printf("get_dist = %d\n", 0);
#endif

		/* Free working vectors */
        	free(SS+lower);
        	free(DD+lower);
        	free(II+lower-1);
		return 0;
	}

	for (c=1; c<=limit; ++c) {
                t = max(lower, start-c);
                t1 = II[t-1];
                for (k=t; k<=min(upper, start+c); ++k) {
                        t2 = II[k];
                        II[k] = max(t1, SS[k]);
                        t1 = t2;
                        DD[k] = max(DD[k+1]+1, SS[k]);
                        SS[k] = snake(k, min(j2-k,max(max(SS[k]+1, II[k]), DD[k]
)), i2, j2);
                }

		if (SS[goal_diag] >= i2) {
#ifdef STATS
			printf("get_dist = %d\n", c);
#endif

			/* Free working vectors */
        		free(SS+lower);
        		free(DD+lower);
        		free(II+lower-1);
			return c;
		}
        }

	/* Ran out of distance limit */
	printf("Two sequences are not really similar.\n");
	printf("Please try exact phase 1.\n");
	exit(1);
}

int dist_only(int *i1, int *j1, int *i2, int *j2)
{
	int *S, *D, *I;
	int *SP, *DP, *IP;
	int c, k;
	int t1, t2;
	int tp1, tp2;

	/* Initial space allocation */
	S = (int *)ckalloc((N+limit)*sizeof(int)) + limit - 1;
	SP = (int *)ckalloc((N+limit)*sizeof(int)) + limit - 1;
	D = (int *)ckalloc((N+limit+1)*sizeof(int)) + limit - 1;
	DP = (int *)ckalloc((N+limit+1)*sizeof(int)) + limit - 1;
	I = (int *)ckalloc((N+limit+1)*sizeof(int)) + limit;
	IP = (int *)ckalloc((N+limit+1)*sizeof(int)) + limit;

	/* Free end-gaps for the shorter sequence */
	for (k=0; k<=N; ++k) {
			S[k] = snake(k,0,M,N);
			SP[k] = k;
			if (S[k] >= M) {
#ifdef STATS
				printf("k=%d, S[k]=%d, SP[k]=%d, distance = %d\n", k, S[k], SP[k], c);
#endif
				*i1 = 0;
				*j1 = SP[k];
				*i2 = M;
				*j2 = min(N, k+M);
				free(S-limit+1);
				free(SP-limit+1);
				free(D-limit+1);
				free(DP-limit+1);
				free(I-limit);
				free(IP-limit);
				return 0;
			}
	}
	for (k=1-limit; k<=-1; ++k) {
		S[k] = MININT;
		SP[k] = 0;
	}
	for (k=-limit; k<=N; ++k) I[k] = MININT;
	for (k=1-limit; k<=N+1; ++k) D[k] = MININT;

	for (c=1; c<=limit; ++c) {
		t1 = I[-c];
		tp1 = IP[-c];
		for (k=1-c; k<=N; ++k) {
			t2 = I[k];
			tp2 = IP[k];
			/* I[k] = max(t1, S[k]); */
			if (t1 > S[k]) {
				I[k] = t1;
				IP[k] = tp1;
			} else {
				I[k] = S[k];
				IP[k] = SP[k];
			}
			t1 = t2;
			tp1 = tp2;

			/* D[k] = max(D[k+1]+1,S[k]); */
			if (D[k+1]+1 > S[k]) {
				D[k] = D[k+1]+1;
				DP[k] = DP[k+1];
			} else {
				D[k] = S[k];
				DP[k] = SP[k];
			}

			/* S[k] = snake(k, min(N-k,max(max(S[k]+1, I[k]), D[k])), M, N); */
			if (S[k]+1 >= I[k]) {
				S[k] = min(N-k, S[k]+1);
			} else {
				S[k] = I[k];
				SP[k] = IP[k];
			}
			if (D[k] > S[k]) {
				S[k] = D[k];
				SP[k] = DP[k];
			}

			S[k] = snake(k, S[k], M, N);

/*
printf("k=%d, S[k]=%d, D[k]=%d, I[k]=%d, distance = %d\n", k, S[k], D[k], I[k], c);
*/

			if (S[k] >= M) {
#ifdef STATS
				printf("k=%d, S[k]=%d, SP[k]=%d, distance = %d\n", k, S[k], SP[k], c);
#endif
				*i1 = 0;
				*j1 = SP[k];
				*i2 = M;
				*j2 = min(N, k+M);
				free(S-limit+1);
				free(SP-limit+1);
				free(D-limit+1);
				free(DP-limit+1);
				free(I-limit);
				free(IP-limit);
				return c;
			}
		}
	}

	/* Ran out of distance limit */
	printf("Two sequences are not really similar\n");
	exit(1);
}

int snake(int k, int x, int endx, int endy)
{
	int y;

	if (x<0) return x;
	y = x+k;
	while (x<endx && y<endy && A[x]==B[y]) {
		++x; ++y;
	}
	return x;
}

int rsnake(int k, int x, int startx, int starty)
{
	int y;

	if (x>M) return x;
	y = x+k;
	while (x>startx && y>starty && A[x-1]==B[y-1]) {
		--x; --y;
	}
	return x;
}

void path(int i1, int j1, char type1, int i2, int j2, char type2, int dist, edit_script **head, edit_script **tail)
{
	int *SS, *DD, *II;	 /* Forward vectors */
	int *RS, *RD, *RI;	 /* Backward vectors */

	edit_script *head1, *tail1, *head2, *tail2;
	int midc, rmidc;
	int start, lower, upper;
	int rstart, rlower, rupper;
	int c, k, t1, t2, t;
	int maxint;
	int mi, mj, mtype;
	char flag;

/*
printf("i1=%d,j1=%d,type1=%d,i2=%d,j2=%d,type2=%d,dist=%d\n",i1,j1,type1,i2,j2,type2,dist);
*/

	/* Boundary cases */
	if (i1 == i2) {
	   if (j1 == j2) *head = NULL;
	   else {
		head1 = (edit_script *) ckalloc(sizeof(edit_script));
		head1->op_type = INS;
		head1->num = j2-j1;
		head1->next = NULL;
		*head = *tail = head1;
	   }
	   return;
	}

	if (j1 == j2) {
		head1 = (edit_script *) ckalloc(sizeof(edit_script));
		head1->op_type = DEL;
		head1->num = i2-i1;
		head1->next = NULL;
		*head = *tail = head1;
		return;
	}

	if (dist <= 1) {
	   if (j2-i2 == j1-i1) {
		head1 = (edit_script *) ckalloc(sizeof(edit_script));
		head1->op_type = SUB;
		head1->num = i2-i1;
		head1->next = NULL;
		*head = *tail = head1;
	   } else if (j2-i2 > j1-i1) {
		if (type1 == INS) {
			head1 = (edit_script *) ckalloc(sizeof(edit_script));
                	head1->op_type = INS;
                	head1->num = 1;
			head2 = (edit_script *) ckalloc(sizeof(edit_script));
                	head2->op_type = SUB;
                	head2->num = i2-i1;
		} else {
			head1 = (edit_script *) ckalloc(sizeof(edit_script));
                	head1->op_type = SUB;
                	head1->num = i2-i1;
			head2 = (edit_script *) ckalloc(sizeof(edit_script));
                	head2->op_type = INS;
                	head2->num = 1;
		}
		head1->next = head2;
		head2->next = NULL;
		*head = head1;
		*tail = head2;
           } else if (j2-i2 < j1-i1) {
		if (type1 == DEL) {
			head1 = (edit_script *) ckalloc(sizeof(edit_script));
                	head1->op_type = DEL;
                	head1->num = 1;
			head2 = (edit_script *) ckalloc(sizeof(edit_script));
                	head2->op_type = SUB;
                	head2->num = j2-j1;
		} else {
			head1 = (edit_script *) ckalloc(sizeof(edit_script));
                	head1->op_type = SUB;
                	head1->num = j2-j1;
			head2 = (edit_script *) ckalloc(sizeof(edit_script));
                	head2->op_type = DEL;
                	head2->num = 1;
		}
		head1->next = head2;
		head2->next = NULL;
		*head = head1;
		*tail = head2;
           }
	   return;
	}

	/* Divide the problem at the middle cost */
	midc = dist/2;
	rmidc = dist - midc;

	/* Compute the boundary diagonals */
	start = j1 - i1;
	lower = max(j1-i2, start-midc);
	upper = min(j2-i1, start+midc);
	rstart = j2-i2;
	rlower = max(j1-i2, rstart-rmidc);
	rupper = min(j2-i1, rstart+rmidc);

	/* Allocate space for forward vectors */
	SS = (int *)ckalloc((upper-lower+1)*sizeof(int)) - lower;
	DD = (int *)ckalloc((upper-lower+2)*sizeof(int)) - lower;
	II = (int *)ckalloc((upper-lower+2)*sizeof(int)) - lower + 1;

	/* Forward computation */
	for (k=lower; k<=upper; ++k) SS[k] = MININT;
	for (k=lower; k<=upper+1; ++k) DD[k] = MININT;
	for (k=lower-1; k<=upper; ++k) II[k] = MININT;
	if (type1 == SUB) SS[start] = snake(start, i1, i2, j2);
	else if (type1 == DEL) {
		DD[start] = i1;
		SS[start] = snake(start,i1,i2,j2);
	} else {
		II[start] = i1;
		SS[start] = snake(start,i1,i2,j2);
	}

	for (c=1; c<=midc; ++c) {
		t = max(lower, start-c);
		t1 = II[t-1];
		for (k=t; k<=min(upper, start+c); ++k) {
			t2 = II[k];
			II[k] = max(t1, SS[k]);
			t1 = t2;
			DD[k] = max(DD[k+1]+1, SS[k]);
			SS[k] = snake(k, min(j2-k,max(max(SS[k]+1, II[k]), DD[k])), i2, j2);
		}
	}

	/* Allocate space for backward vectors */
	RS = (int *)ckalloc((rupper-rlower+1)*sizeof(int)) - rlower;
	RD = (int *)ckalloc((rupper-rlower+2)*sizeof(int)) - rlower + 1;
	RI = (int *)ckalloc((rupper-rlower+2)*sizeof(int)) - rlower;

	/* Backward computation */
	maxint = i2 + dist + N;
	for (k=rlower; k<=rupper; ++k) RS[k] = maxint;
	for (k=rlower-1; k<=rupper; ++k) RD[k] = maxint;
	for (k=rlower; k<=rupper+1; ++k) RI[k] = maxint;
	if (type2 == SUB)
	   RI[rstart] = RD[rstart] = RS[rstart] = rsnake(rstart, i2, i1, j1);
	else if (type2 == DEL) RD[rstart] = i2;
	else RI[rstart] = i2;

	for (c=1; c<=rmidc; ++c) {
		t = max(rlower, rstart-c);
		t1 = RD[t-1];
		for (k=t; k<=min(rupper, rstart+c); ++k) {
			RS[k] = rsnake(k, max(j1-k, min(min(RS[k]-1,RD[k]),RI[k])),i1,j1);
			t2 = RD[k];
			RD[k] = min(t1-1, RS[k]);
			t1 = t2;
			RI[k] = min(RI[k+1], RS[k]);
		}
	}

	/* Find (mi, mj, mtype) such that
	       the distance from (i1, j1, type1) to (mi, mj, mtype) is midc
	   and the distance from (mi, mj, mtype) to (i2, j2, type2) is rmidc.
	*/

	flag = FALSE;
	for (k=max(lower,rlower); k<=min(upper,rupper);++k) {

/*
printf("k=%d, SS=%d, RS=%d, DD=%d, RD=%d, II=%d, RI=%d\n",k,SS[k],RS[k],DD[k],RD[k],II[k],RI[k]);
*/

		if (SS[k]>=RS[k] || DD[k]>=RD[k] || II[k]>=RI[k]) {
			if (DD[k]>=RD[k]) {
				mi = DD[k];
				mj = k+mi;
				mtype = DEL;
			} else if (II[k] >= RI[k]) {
				mi = II[k];
				mj = k+mi;
				mtype = INS;
			} else {
				mi = SS[k];
				mj = k+mi;
				mtype = SUB;
			}

/*
			printf("mi=%d, mj=%d, mtype=%d\n", mi, mj, mtype);
*/
			flag = TRUE;
			break;
		}
	}

	/* Free working vectors */
	free(SS+lower);
	free(DD+lower);
	free(II+lower-1);
	free(RS+rlower);
	free(RD+rlower-1);
	free(RI+rlower);

	if (flag) {
		/* Find a path from (i1,j1,type1) to (mi,mj,mtype) */
		path(i1,j1,type1,mi,mj,mtype,midc,&head1,&tail1);

		/* Find a path from (mi,mj,mtype) to (i2,j2,type2) */
		path(mi,mj,mtype,i2,j2,type2,rmidc,&head2,&tail2);

		/* Join these two paths together */
		if (head1) tail1->next = head2;
		else head1 = head2;
	} else {
		printf("Something wrong when dividing\n");
		head1 = NULL;
	}
	*head = head1;
	if (head2) *tail = tail2;
	else *tail = tail1;
}

/* Condense_script - merge contiguous operations of the same type together */
void Condense_script(edit_script *head)
{
	edit_script *tp, *tp1;

	tp = head;
	while (tp != NULL) {
	   while (((tp1 = tp->next) != NULL) && (tp->op_type == tp1->op_type)) {
		 tp->num = tp->num + tp1->num;
		 tp->next = tp1->next;
		 free(tp1);
	   }
	   tp = tp->next;
	}
}

/* Check_script - check the total distance described in the script */
int Check_script(edit_script *head, int I1, int J1)
{
	int i, j, t, dist;
	edit_script *tp;

	dist = 0;
	i = I1;
	j = J1;
	tp = head;
	while (tp != NULL) {

/*
printf("tp->op_type=%d, tp->num=%d\n", tp->op_type, tp->num);
*/

	   if (tp->op_type == SUB) {
		for (t=0; t<tp->num; ++t) {
		   if (A[i] != B[j]) ++dist;
		   ++i; ++j;	
		}
	   } else if (tp->op_type == INS) {
		dist += 1+tp->num;
		j += tp->num;
	   } else {	/* DEL */
		dist += 1+tp->num;
		i += tp->num;
	   }
	   tp = tp->next;
	}

	return dist;
}

void Reverse_Script(edit_script *head)
{
	edit_script *tp;

	tp = head;
	while (tp != NULL) {
	   if (tp->op_type == INS) tp->op_type = DEL;
	   else if (tp->op_type == DEL) tp->op_type = INS;
	   tp = tp->next;
	}
}

void Print_Script(edit_script *head, int M, int N)
{
	edit_script *tp;
	int i, j, k, count;

	i = j = 0;
	tp = head;
	while (tp != NULL) {
	   if (tp->op_type == SUB) {
		k = 0;
		while (k < tp->num) {
			count = 0;
			while (A[i] == B[j]) {
				++i;
				++j;
				++count;
				++k;
			}
			if (count > 0) printf("copy %d\n", count);
			if (k < tp->num) {
				printf("replace %c by %c\n", A[i], B[j]);
				++i;
				++j;
				++k;
			}
		}
		/*
		if (tp->num > 1) printf("%d substitutions\n", tp->num);
		else printf("%d substitution\n", tp->num);
		*/
	   } else if (tp->op_type == INS) {
		if ((tp==head || tp->next==NULL) && (M <= N))
			printf("skip (second sequence) %d\n", tp->num);
		else {
			/*
			printf("insert %d\n", tp->num);
			*/
			printf("insert ");
			for (k=j; k<j+tp->num; ++k) printf("%c", B[k]);
			printf("\n");
		}
		j += tp->num;
	   } else {	/* DEL */
		if ((tp==head || tp->next==NULL) && (M > N))
			printf("skip (first sequence) %d\n", tp->num);
		else {
			/*
			printf("delete %d\n", tp->num);
			*/
			printf("delete ");
			for (k=i; k<i+tp->num; ++k) printf("%c", A[k]);
			printf("\n");
		}
		i += tp->num;
	   }
	   tp = tp->next;
	}
}

void S2A(edit_script *head, int *S)
{
	edit_script *tp;
	int *lastS, i;

	tp = head;
	lastS = S;
	while (tp != NULL) {
/*
	printf("tp->op_type=%d, tp->num=%d\n",tp->op_type, tp->num);
*/
	   if (tp->op_type == SUB) {
		for (i=0; i<tp->num; ++i) *lastS++ = 0;
	   } else if (tp->op_type == INS) {
		*lastS++ = tp->num;
	   } else {	/* DEL */
		*lastS++ = 0 - tp->num;
	   }
	   tp = tp->next;
	}
}

/* Alignment display routine */

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

void DISPLAY(char A[], char B[], int M, int N, int S[], int AP, int BP)
{ register char *a, *b, *c;
  register int   i,  j, op;
           int   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];
          *c++ = (*a++ == *b++) ? '|' : ' ';
        }
      else
        { if (op == 0)
            op = *S++;
          if (op > 0)
            { *a++ = ' ';
              *b++ = B[++j];
              op--;
            }
          else
            { *a++ = A[++i];
              *b++ = ' ';
              op++;
            }
          *c++ = '-';
        }
      if (a >= ALINE+50 || i >= M && j >= N)
        { *a = *b = *c = '\0';
          printf("\n%6d ",50*lines++);
          for (b = ALINE+10; b <= a; b += 10)
            printf("    .    :");
          if (b <= a+5)
            printf("    .");
          printf("\n%6d %s\n       %s\n%6d %s\n",ap,ALINE,CLINE,bp,BLINE);
	  ap = AP + i;
	  bp = BP + j;
          a = ALINE;
          b = BLINE;
          c = CLINE;
        }
    }
}

void Script_free(edit_script *head)
{
	edit_script *tp, *tp1;

	tp = head;
	while (tp != NULL) {
	   tp1 = tp->next;
	   free(tp);
	   tp = tp1;
	}
}