ta=ckfree((void *)ta);   return(mask);}static lint prfscore(sint n, sint m){   sint    ix;   lint  score;   score = 0.0;   for (ix=0; ix<=max_aa; ix++)     {         score += (profile1[n][ix] * profile2[m][ix]);     }   score += (profile1[n][gap_pos1] * profile2[m][gap_pos1]);   score += (profile1[n][gap_pos2] * profile2[m][gap_pos2]);   return(score/10);   }static void ptracepath(sint *alen){    sint i,j,k,pos,to_do;    pos = 0;    to_do=print_ptr-1;    for(i=1;i<=to_do;++i) {if (debug>1) fprintf(stdout,"%d ",(pint)displ[i]);            if(displ[i]==0) {                    aln_path1[pos]=2;                    aln_path2[pos]=2;                    ++pos;            }            else {                    if((k=displ[i])>0) {                            for(j=0;j<=k-1;++j) {                                    aln_path2[pos+j]=2;                                    aln_path1[pos+j]=1;                            }                            pos += k;                    }                    else {                            k = (displ[i]<0) ? displ[i] * -1 : displ[i];                            for(j=0;j<=k-1;++j) {                                    aln_path1[pos+j]=2;                                    aln_path2[pos+j]=1;                            }                            pos += k;                    }            }    }if (debug>1) fprintf(stdout,"\n");   (*alen) = pos;}static void pdel(sint k){        if(last_print<0)                last_print = displ[print_ptr-1] -= k;        else                last_print = displ[print_ptr++] = -(k);}static void padd(sint k){        if(last_print<0) {                displ[print_ptr-1] = k;                displ[print_ptr++] = last_print;        }        else                last_print = displ[print_ptr++] = k;}static void palign(void){        displ[print_ptr++] = last_print = 0;}static lint pdiff(sint A,sint B,sint M,sint N,sint go1, sint go2){        sint midi,midj,type;        lint midh;        static lint t, tl, g, h;{		static sint i,j;        static lint hh, f, e, s;/* Boundary cases: M <= 1 or N == 0 */if (debug>2) fprintf(stdout,"A %d B %d M %d N %d midi %d go1 %d go2 %d\n", (pint)A,(pint)B,(pint)M,(pint)N,(pint)M/2,(pint)go1,(pint)go2);/* if sequence B is empty....                                            */        if(N<=0)  {/* if sequence A is not empty....                                        */                if(M>0) {/* delete residues A[1] to A[M]                                          */                        pdel(M);                }                return(-gap_penalty1(A,B,M));        }/* if sequence A is empty....                                            */        if(M<=1) {                if(M<=0) {/* insert residues B[1] to B[N]                                          */                        padd(N);                        return(-gap_penalty2(A,B,N));                }/* if sequence A has just one residue....                                */                if (go1 == 0)                	midh =  -gap_penalty1(A+1,B+1,N);                else                	midh =  -gap_penalty2(A+1,B,1)-gap_penalty1(A+1,B+1,N);                midj = 0;                for(j=1;j<=N;j++) {                        hh = -gap_penalty1(A,B+1,j-1) + prfscore(A+1,B+j)                            -gap_penalty1(A+1,B+j+1,N-j);                        if(hh>midh) {                                midh = hh;                                midj = j;                        }                }                if(midj==0) {                        padd(N);                        pdel(1);                }                else {                        if(midj>1) padd(midj-1);                        palign();                        if(midj<N) padd(N-midj);                }                return midh;        }/* Divide sequence A in half: midi */        midi = M / 2;/* In a forward phase, calculate all HH[j] and HH[j] */        HH[0] = 0.0;        t = -open_penalty1(A,B+1);        tl = -ext_penalty1(A,B+1);        for(j=1;j<=N;j++) {                HH[j] = t = t+tl;                DD[j] = t-open_penalty2(A+1,B+j);        }		if (go1 == 0) t = 0;		else t = -open_penalty2(A+1,B);        tl = -ext_penalty2(A+1,B);        for(i=1;i<=midi;i++) {                s = HH[0];                HH[0] = hh = t = t+tl;                f = t-open_penalty1(A+i,B+1);                for(j=1;j<=N;j++) {                	g = open_penalty1(A+i,B+j);                	h = ext_penalty1(A+i,B+j);                        if ((hh=hh-g-h) > (f=f-h)) f=hh;                	g = open_penalty2(A+i,B+j);                	h = ext_penalty2(A+i,B+j);                        if ((hh=HH[j]-g-h) > (e=DD[j]-h)) e=hh;                        hh = s + prfscore(A+i, B+j);                        if (f>hh) hh = f;                        if (e>hh) hh = e;                        s = HH[j];                        HH[j] = hh;                        DD[j] = e;                }        }        DD[0]=HH[0];/* In a reverse phase, calculate all RR[j] and SS[j] */        RR[N]=0.0;        tl = 0.0;        for(j=N-1;j>=0;j--) {                g = -open_penalty1(A+M,B+j+1);                tl -= ext_penalty1(A+M,B+j+1);                RR[j] = g+tl;                SS[j] = RR[j]-open_penalty2(A+M,B+j);                gS[j] = open_penalty2(A+M,B+j);        }        tl = 0.0;        for(i=M-1;i>=midi;i--) {                s = RR[N];                if (go2 == 0) g = 0;                else g = -open_penalty2(A+i+1,B+N);                tl -= ext_penalty2(A+i+1,B+N);                RR[N] = hh = g+tl;                t = open_penalty1(A+i,B+N);                f = RR[N]-t;                for(j=N-1;j>=0;j--) {                	g = open_penalty1(A+i,B+j+1);                	h = ext_penalty1(A+i,B+j+1);                        if ((hh=hh-g-h) > (f=f-h-g+t)) f=hh;                        t = g;                	g = open_penalty2(A+i+1,B+j);                	h = ext_penalty2(A+i+1,B+j);                        hh=RR[j]-g-h;                        if (i==(M-1)) {				 e=SS[j]-h;			}                        else {				e=SS[j]-h-g+open_penalty2(A+i+2,B+j);				gS[j] = g;			}                        if (hh > e) e=hh;                        hh = s + prfscore(A+i+1, B+j+1);                        if (f>hh) hh = f;                        if (e>hh) hh = e;                        s = RR[j];                        RR[j] = hh;                        SS[j] = e;                }        }        SS[N]=RR[N];        gS[N] = open_penalty2(A+midi+1,B+N);/* find midj, such that HH[j]+RR[j] or DD[j]+SS[j]+gap is the maximum */        midh=HH[0]+RR[0];        midj=0;        type=1;        for(j=0;j<=N;j++) {                hh = HH[j] + RR[j];                if(hh>=midh)                        if(hh>midh || (HH[j]!=DD[j] && RR[j]==SS[j])) {                                midh=hh;                                midj=j;                        }        }        for(j=N;j>=0;j--) {                hh = DD[j] + SS[j] + gS[j];                if(hh>midh) {                        midh=hh;                        midj=j;                        type=2;                }        }}/* Conquer recursively around midpoint                                   */        if(type==1) {             /* Type 1 gaps  */if (debug>2) fprintf(stdout,"Type 1,1: midj %d\n",(pint)midj);                pdiff(A,B,midi,midj,go1,1);if (debug>2) fprintf(stdout,"Type 1,2: midj %d\n",(pint)midj);                pdiff(A+midi,B+midj,M-midi,N-midj,1,go2);        }        else {if (debug>2) fprintf(stdout,"Type 2,1: midj %d\n",(pint)midj);                pdiff(A,B,midi-1,midj,go1, 0);                pdel(2);if (debug>2) fprintf(stdout,"Type 2,2: midj %d\n",(pint)midj);                pdiff(A+midi+1,B+midj,M-midi-1,N-midj,0,go2);        }        return midh;       /* Return the score of the best alignment */}/* calculate the score for opening a gap at residues A[i] and B[j]       */static sint open_penalty1(sint i, sint j){   sint g;   if (!endgappenalties &&(i==0 || i==prf_length1)) return(0);   g = profile2[j][GAPCOL] + profile1[i][GAPCOL];   return(g);}/* calculate the score for extending an existing gap at A[i] and B[j]    */static sint ext_penalty1(sint i, sint j){   sint h;   if (!endgappenalties &&(i==0 || i==prf_length1)) return(0);   h = profile2[j][LENCOL];   return(h);}/* calculate the score for a gap of length k, at residues A[i] and B[j]  */static sint gap_penalty1(sint i, sint j, sint k){   sint ix;   sint gp;   sint g, h = 0;   if (k <= 0) return(0);   if (!endgappenalties &&(i==0 || i==prf_length1)) return(0);   g = profile2[j][GAPCOL] + profile1[i][GAPCOL];   for (ix=0;ix<k && ix+j<prf_length2;ix++)      h += profile2[ix+j][LENCOL];   gp = g + h;   return(gp);}/* calculate the score for opening a gap at residues A[i] and B[j]       */static sint open_penalty2(sint i, sint j){   sint g;   if (!endgappenalties &&(j==0 || j==prf_length2)) return(0);   g = profile1[i][GAPCOL] + profile2[j][GAPCOL];   return(g);}/* calculate the score for extending an existing gap at A[i] and B[j]    */static sint ext_penalty2(sint i, sint j){   sint h;   if (!endgappenalties &&(j==0 || j==prf_length2)) return(0);   h = profile1[i][LENCOL];   return(h);}/* calculate the score for a gap of length k, at residues A[i] and B[j]  */static sint gap_penalty2(sint i, sint j, sint k){   sint ix;   sint gp;   sint g, h = 0;   if (k <= 0) return(0);   if (!endgappenalties &&(j==0 || j==prf_length2)) return(0);   g = profile1[i][GAPCOL] + profile2[j][GAPCOL];   for (ix=0;ix<k && ix+i<prf_length1;ix++)      h += profile1[ix+i][LENCOL];   gp = g + h;   return(gp);}