DNA雙序列對比

#include<cstdio>
#include<cstring>
#include<bitset>
using std::bitset;

#define BOTH 0
#define S1 1
#define S2 2
const int Standard[4][4]=
{{10,-1,-3,-4},
 {-1,7,-5,-3},
 {-3,-5,9,0},
 {-4,-3,0,8}};
#define GAP -5
int Translate(char a)
{
	switch(a)
	{
		case 'A':	return 0;break;
		case 'G':	return 1;break;
		case 'C':	return 2;break;
		case 'T':	return 3;break;
	}
	return -1;
}

char io1[5000],io2[5000],*Seq1,*Seq2,Ans1[5000],Ans2[5000];
int Point[5000][5000];
bitset<5000> From[3][5000];//From[Seq1/Seq2?][x][y]
int Len1,Len2;
void Update(int x,int y,int&temp,int mode,int val)
{
	if(val>=temp)
		temp=val,
		From[BOTH][x][y]=false,
		From[S1][x][y]=false,
		From[S2][x][y]=false,
		From[mode][x][y]=true;
	return;
}
int GetPoint(int s1n,int s2n)
	{return Standard[Translate(Seq1[s1n])][Translate(Seq2[s2n])];}

int main()
{
	int temp=0;
	scanf("%s %s",io1,io2);
	Len1=strlen(io1);
	Len2=strlen(io2);
	if(Len1>Len2)	Seq2=io1,Seq1=io2,Len2=Len1,Len1=strlen(io2);
	else			Seq1=io1,Seq2=io2;
	//Seq1 short, Seq2 long

	Point[0][0]=0;
	for(int i=0;i<=Len1;i++)	Point[i][0]=0;
	for(int i=0;i<=Len2;i++)	Point[0][i]=0;
	for(int i=0;i<5000;i++)	Ans1[i]=Ans2[i]='\0';

	for(int i=1;i<=Len1;i++)
		for(int j=1;j<=Len2;j++)
		{
			temp=0;
			Update(i,j,temp,BOTH,Point[i-1][j-1]+GetPoint(i-1,j-1));
			Update(i,j,temp,S1,Point[i][j-1]+GAP);
			Update(i,j,temp,S2,Point[i-1][j]+GAP);
			Point[i][j]=temp;
		}

	int x=Len1,y=Len2,Pnt=0;
	while(x!=0&&y!=0)
	{
		if(From[BOTH][x][y])	x--,y--,Ans1[Pnt]=Seq1[x],Ans2[Pnt]=Seq2[y];
		else if(From[S2][x][y])	x--,Ans1[Pnt]=Seq1[x],Ans2[Pnt]='_';
		else if(From[S1][x][y])	y--,Ans1[Pnt]='_',Ans2[Pnt]=Seq2[y];
		else break;
		Pnt++;
	}
	while(x!=0||y!=0)
	{
		if(x==0)Ans1[Pnt]='_';
		else	Ans1[Pnt]=Seq1[--x];
		if(y==0)Ans2[Pnt]='_';
		else	Ans2[Pnt]=Seq2[--y];
		Pnt++;
	}

	for(int i=Pnt-1;i>=0;i--)	putchar(Ans1[i]);
	putchar('\n');
	for(int i=Pnt-1;i>=0;i--)	putchar(Ans2[i]);
	return 0;
}