j = lpj->n;		/* calculate sum_rows[j] */		for (lpi=tvalid[0].next ; lpi->n < j ; lpi = lpi->next) {			i = lpi->n;			tmp_sum += tmat[i][j];			/* tmat[j][i] = tmat[i][j]; */		}		sum_rows[j] = tmp_sum;		tmp_sum = 0.0;		/* Set lpi to that lpi->n is greater than j */		if ((lpi != NULL) && (lpi->n == j)) {			lpi = lpi->next;		}		/* calculate sum_cols[j] */		for( ; lpi!=NULL ; lpi = lpi->next) {			i = lpi->n;			tmp_sum += tmat[j][i];			/* tmat[i][j] = tmat[j][i]; */		}		sum_cols[j] = tmp_sum;	}/*********************** Enter The Main Cycle ***************************/	for(nc=1, loop_limit = (last_seq-first_seq+1-3); nc<=loop_limit; ++nc) {		sumd = 0.0;		/* IMPROVEMENT 1, STEP 4 : use sum value */		for(lpj=tvalid[0].next ; lpj!=NULL ; lpj = lpj->next) {			sumd += sum_cols[lpj->n];		}		/* IMPROVEMENT 3, STEP 0 : multiply tmin and 2*fnseqs2 */		fnseqs2 = fnseqs - 2.0;		/* Set fnseqs2 at this point. */		tmin = 99999.0 * 2.0 * fnseqs2;/*.................compute SMATij values and find the smallest one ........*/		mini = minj = 0;		/* jj must starts at least 2 */		if ((tvalid[0].next != NULL) && (tvalid[0].next->n == 1)) {			lpjj = tvalid[0].next->next;		} else {			lpjj = tvalid[0].next;		}		for( ; lpjj != NULL; lpjj = lpjj->next) {			jj = lpjj->n;			for(lpii=tvalid[0].next ; lpii->n < jj ; lpii = lpii->next) {				ii = lpii->n;				diq = djq = 0.0;				/* IMPROVEMENT 1, STEP 4 : use sum value */				diq = sum_cols[ii] + sum_rows[ii];				djq = sum_cols[jj] + sum_rows[jj];				/*				 * always ii < jj in this point. Use upper				 * triangle of score matrix.				 */				dij = tmat[ii][jj];				/*				 * IMPROVEMENT 3, STEP 1 : fnseqs2 is				 * already calculated.				 */				/* fnseqs2 = fnseqs - 2.0 */				/* IMPROVEMENT 4 : transform the equation */  /*-------------------------------------------------------------------*   * OPTIMIZE of expression 'total = d2r + fnseqs2*dij + dr*2.0'       *   * total = d2r + fnseq2*dij + 2.0*dr                                 *   *       = d2r + fnseq2*dij + 2(sumd - dij - d2r)                    *   *       = d2r + fnseq2*dij + 2*sumd - 2*dij - 2*d2r                 *   *       =       fnseq2*dij + 2*sumd - 2*dij - 2*d2r + d2r           *   *       = fnseq2*dij + 2*sumd - 2*dij - d2r                         *   *       = fnseq2*dij + 2*sumd - 2*dij - (diq + djq - 2*dij)         *   *       = fnseq2*dij + 2*sumd - 2*dij - diq - djq + 2*dij           *   *       = fnseq2*dij + 2*sumd - 2*dij + 2*dij - diq - djq           *   *       = fnseq2*dij + 2*sumd  - diq - djq                          *   *-------------------------------------------------------------------*/				total = fnseqs2*dij + 2.0*sumd  - diq - djq;				/* 				 * IMPROVEMENT 3, STEP 2 : abbrevlate				 * the division on comparison between 				 * total and tmin.				 */				/* total = total / (2.0*fnseqs2); */				if(total < tmin) {					tmin = total;					mini = ii;					minj = jj;				}			}		}		/* MEMO: always ii < jj in avobe loop, so mini < minj *//*.................compute branch lengths and print the results ........*/		dio = djo = 0.0;		/* IMPROVEMENT 1, STEP 4 : use sum value */		dio = sum_cols[mini] + sum_rows[mini];		djo = sum_cols[minj] + sum_rows[minj];		dmin = tmat[mini][minj];		dio = (dio - dmin) / fnseqs2;		djo = (djo - dmin) / fnseqs2;		bi = (dmin + dio - djo) * 0.5;		bj = dmin - bi;		bi = bi - av[mini];		bj = bj - av[minj];		if( av[mini] > 0.0 )			typei = 0;		else			typei = 1;		if( av[minj] > 0.0 )			typej = 0;		else			typej = 1;		if(verbose) 	 	    fprintf(tree,"\n Cycle%4d     = ",(pint)nc);/*    set negative branch lengths to zero.  Also set any tiny positive   branch lengths to zero.*/		if( fabs(bi) < 0.0001) bi = 0.0;		if( fabs(bj) < 0.0001) bj = 0.0;	    	if(verbose) {		    if(typei == 0) 			fprintf(tree,"Node:%4d (%9.5f) joins ",(pint)mini,bi);		    else 			fprintf(tree," SEQ:%4d (%9.5f) joins ",(pint)mini,bi);		    if(typej == 0) 			fprintf(tree,"Node:%4d (%9.5f)",(pint)minj,bj);		    else 			fprintf(tree," SEQ:%4d (%9.5f)",(pint)minj,bj);		    fprintf(tree,"\n");	    	}		    	left_branch[nc] = bi;	    	right_branch[nc] = bj;		for(i=1; i<=last_seq-first_seq+1; i++)			tree_description[nc][i] = 0;	     	if(typei == 0) { 			for(i=nc-1; i>=1; i--)				if(tree_description[i][mini] == 1) {					for(j=1; j<=last_seq-first_seq+1; j++)  					     if(tree_description[i][j] == 1)						    tree_description[nc][j] = 1;					break;				}		}		else			tree_description[nc][mini] = 1;		if(typej == 0) {			for(i=nc-1; i>=1; i--) 				if(tree_description[i][minj] == 1) {					for(j=1; j<=last_seq-first_seq+1; j++)  					     if(tree_description[i][j] == 1)						    tree_description[nc][j] = 1;					break;				}		}		else			tree_description[nc][minj] = 1;			/*    Here is where the -0.00005 branch lengths come from for 3 or more   identical seqs.*//*		if(dmin <= 0.0) dmin = 0.0001; */                if(dmin <= 0.0) dmin = 0.000001;		av[mini] = dmin * 0.5;/*........................Re-initialisation................................*/		fnseqs = fnseqs - 1.0;		tkill[minj] = 1;		/* IMPROVEMENT 2, STEP 3 : Remove tvalid[minj] from chain list. */		/* [ Before ]		 *  +---------+        +---------+        +---------+       		 *  |prev     |<-------|prev     |<-------|prev     |<---		 *  |    n    |        | n(=minj)|        |    n    |		 *  |     next|------->|     next|------->|     next|----		 *  +---------+        +---------+        +---------+ 		 *		 * [ After ]		 *  +---------+                           +---------+       		 *  |prev     |<--------------------------|prev     |<---		 *  |    n    |                           |    n    |		 *  |     next|-------------------------->|     next|----		 *  +---------+                           +---------+ 		 *                     +---------+		 *              NULL---|prev     |		 *                     | n(=minj)|		 *                     |     next|---NULL		 *                     +---------+ 		 */		(tvalid[minj].prev)->next = tvalid[minj].next;		if (tvalid[minj].next != NULL) {			(tvalid[minj].next)->prev = tvalid[minj].prev;		}		tvalid[minj].prev = tvalid[minj].next = NULL;		/* IMPROVEMENT 1, STEP 5 : re-calculate sum values. */		for(lpj=tvalid[0].next ; lpj != NULL ; lpj = lpj->next) {			double tmp_di = 0.0;			double tmp_dj = 0.0;			j = lpj->n;			/* 			 * subtrace a score value related with 'minj' from			 * sum arrays .			 */			if (j < minj) {				tmp_dj = tmat[j][minj];				sum_cols[j] -= tmp_dj;			} else if (j > minj) {				tmp_dj = tmat[minj][j];				sum_rows[j] -= tmp_dj;			} /* nothing to do when j is equal to minj. */						/* 			 * subtrace a score value related with 'mini' from			 * sum arrays .			 */			if (j < mini) {				tmp_di = tmat[j][mini];				sum_cols[j] -= tmp_di;			} else if (j > mini) {				tmp_di = tmat[mini][j];				sum_rows[j] -= tmp_di;			} /* nothing to do when j is equal to mini. */			/* 			 * calculate a score value of the new inner node.			 * then, store it temporary to join[] array.			 */			join[j] = (tmp_dj + tmp_di) * 0.5;		}		/* 		 * 1)		 * Set the score values (stored in join[]) into the matrix,		 * row/column position is 'mini'.		 * 2)		 * Add a score value of the new inner node to sum arrays.		 */		for(lpj=tvalid[0].next ; lpj != NULL; lpj = lpj->next) {			j = lpj->n;			if (j < mini) {				tmat[j][mini] = join[j];				sum_cols[j] += join[j];			} else if (j > mini) {				tmat[mini][j] = join[j];				sum_rows[j] += join[j];			} /* nothing to do when j is equal to mini. */		}		/* Re-calculate sum_rows[mini],sum_cols[mini]. */		sum_cols[mini] = sum_rows[mini] = 0.0;		/* calculate sum_rows[mini] */		da = 0.0;		for(lpj=tvalid[0].next ; lpj->n < mini ; lpj = lpj->next) {                      da += join[lpj->n];		}		sum_rows[mini] = da;		/* skip if 'lpj->n' is equal to 'mini' */		if ((lpj != NULL) && (lpj->n == mini)) {			lpj = lpj->next;		}		/* calculate sum_cols[mini] */		da = 0.0;		for( ; lpj != NULL; lpj = lpj->next) {                      da += join[lpj->n];		}		sum_cols[mini] = da;		/*		 * Clean up sum_rows[minj], sum_cols[minj] and score matrix		 * related with 'minj'.		 */		sum_cols[minj] = sum_rows[minj] = 0.0;		for(j=1; j<=last_seq-first_seq+1; ++j)			tmat[minj][j] = tmat[j][minj] = join[j] = 0.0;/****/	}						/**end main cycle**//******************************Last Cycle (3 Seqs. left)********************/	nude = 1;	for(lpi=tvalid[0].next; lpi != NULL; lpi = lpi->next) {		l[nude] = lpi->n;		++nude;	}	b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;	b2 =  tmat[l[1]][l[2]] - b1;	b3 =  tmat[l[1]][l[3]] - b1; 	branch[1] = b1 - av[l[1]];	branch[2] = b2 - av[l[2]];	branch[3] = b3 - av[l[3]];/* Reset tiny negative and positive branch lengths to zero */	if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;	if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;	if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;	left_branch[last_seq-first_seq+1-2] = branch[1];	left_branch[last_seq-first_seq+1-1] = branch[2];	left_branch[last_seq-first_seq+1]   = branch[3];	for(i=1; i<=last_seq-first_seq+1; i++)		tree_description[last_seq-first_seq+1-2][i] = 0;	if(verbose)		fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",(pint)nc);	for(i=1; i<=3; ++i) {	   if( av[l[i]] > 0.0) {	      	if(verbose)	      	    fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",(pint)l[i],branch[i]);		for(k=last_seq-first_seq+1-3; k>=1; k--)			if(tree_description[k][l[i]] == 1) {				for(j=1; j<=last_seq-first_seq+1; j++)				 	if(tree_description[k][j] == 1)					    tree_description[last_seq-first_seq+1-2][j] = i;				break;			}	   }	   else  {	      	if(verbose)	   	    fprintf(tree,"\n\t\t  SEQ:%4d (%9.5f) ",(pint)l[i],branch[i]);		tree_description[last_seq-first_seq+1-2][l[i]] = i;	   }	   if(i < 3) {	      	if(verbose)	            fprintf(tree,"joins");	   }	}	if(verbose)		fprintf(tree,"\n");		/* IMPROVEMENT 1, STEP 6 : release memory area */	ckfree(sum_cols);	ckfree(sum_rows);	ckfree(join);	/* IMPROVEMENT 2, STEP 4 : release memory area */	ckfree(tvalid);}#endif /* ORIGINAL_NJ_TREE */void bootstrap_tree(char *phylip_name,char *clustal_name, char *nexus_name){	sint i,j;	int ranno;	char path[MAXLINE+1];    char dummy[10];	char err_mess[1024];	static char **sample_tree;	static char **standard_tree;	static char **save_tree;	sint total_dists, overspill = 0, total_overspill = 0;	sint nfails = 0;	if(empty) {		error("You must load an alignment first");		return;	}        if(nseqs<4) {                error("Alignment has only %d sequences",nseqs);                return;        }	if(!output_tree_clustal && !output_tree_phylip && !output_tree_nexus) {		error("You must select either clustal or phylip or nexus tree output format");		return;	}	get_path(seqname, path);		if (output_tree_clustal) {        if (clustal_name[0]!=EOS) {                if((clustal_phy_tree_file = open_explicit_file(                clustal_name))==NULL) return;        }        else {		if((clustal_phy_tree_file = open_output_file(		"\nEnter name for bootstrap output file  ",path,		clustal_name,"njb")) == NULL) return;        }	}	first_seq=1;	last_seq=nseqs;	if (output_tree_phylip) {        if (phylip_name[0]!=EOS) {                if((phylip_phy_tree_file = open_explicit_file(                phylip_name))==NULL) return;        }	else {		if((phylip_phy_tree_file = open_output_file(