diff --git a/Assembly.cpp b/Assembly.cpp
index e0cc482..7d47d6e 100644
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -721,22 +721,7 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
     {
-        /**
-        if (i < 4)
-        {
-            continue;
-        }
 
-        if (i > 4)
-        {
-            break;
-        }
-        fprintf(stderr, "i: %u\n", i);
-        fflush(stderr);
-        **/
-        
-        
-        
         clear_Heap(&heap);
         clear_Candidates_list(&l);
         ///clear_Candidates_list(&debug_l);
diff --git a/Correct.cpp b/Correct.cpp
index c8970c6..6f8565d 100644
--- a/Correct.cpp
+++ b/Correct.cpp
@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
+#include <pthread.h>
 #include "Correct.h"
 #include "Levenshtein_distance.h"
 #include "edlib.h"
@@ -8,6 +9,7 @@
 long long T_total_match=0;
 long long T_total_unmatch=0;
 long long T_total_mis=0;
+pthread_mutex_t debug_statistics ;
 
 
 #define MAX(x, y) ((x >= y)?x:y)
@@ -286,7 +288,6 @@ long long o_len, int threashold, int error, long long* total_mis)
 
     if (result.status == EDLIB_STATUS_OK) {
         if (result.editDistance != error)
-        ///if (result.editDistance != error && result.endLocations[0] > x_len)
         {   
             
             (*total_mis)++;
@@ -310,36 +311,61 @@ long long o_len, int threashold, int error, long long* total_mis)
             char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD);
             int cigar_length = strlen(cigar);
 
-            /**
-            for (int i = cigar_length - 1; i >= 0; i--)
+            int i = cigar_length - 1;
+            int j = 0;
+            char tmp;
+           
+
+            while (i >= 0)
             {
-                switch (cigar[i])
+                if (up_length < 0 || left_length < 0)
                 {
-                case 'I':
-                    
-                    break;
-                
-                default:
                     break;
                 }
-            }
-            **/
-            
-            ///fprintf(stderr,"%s\n", cigar);
-            free(cigar);
-            
-            /**
-            fprintf(stderr, "****\ni: %u, edlib: %d, alignmentLength: %d, startLocations: %d, endLocations: %d\n", 
-            i, result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0]);
-            char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD);
-            fprintf(stderr,"%s\n", cigar);
-            free(cigar);
-            print_string(x_string, x_len);
-            print_string(y_string, o_len);
 
-            fprintf(stderr, "BPM: %d\n", error);
-            **/
+                if (cigar[i] == 'I' || cigar[i] == 'D')
+                {
+                    tmp = cigar[i];
+                    cigar[i] = '\0';
+                    j = i - 1;
+                    while (cigar[j] <= '9' && cigar[j] >= '0' && j >= 0)
+                    {
+                        j--;
+                    }
+                    j++;
+                    int Len = atoi(cigar + j);
+                    cigar[i] = tmp;
+                    i = j - 1;
+
+                    if (tmp == 'I')
+                    {
+                        up_length = up_length - Len;
+                        left_length = left_length + Len;
+                    }
+                    else
+                    {
+                        left_length = left_length - Len;
+                        up_length = up_length + Len;
+                    }
+                }
+                else
+                {
+                    i--;
+                }
+            }
+
             
+            if (up_length >= 0 && left_length >= 0)
+            {
+                fprintf(stderr, "****\nedlib: %d, alignmentLength: %d, startLocations: %d, endLocations: %d\n", 
+                result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0]);
+                fprintf(stderr,"%s\n", cigar);
+                print_string(x_string, x_len);
+                print_string(y_string, o_len);
+                fprintf(stderr, "BPM: %d\n", error);
+            }
+
+            free(cigar);
         }
     }
     edlibFreeAlignResult(result);
@@ -358,10 +384,14 @@ char* r_string)
     long long x_end, x_len;
     int end_site;
     unsigned int error;
-    long long total_match=0;
-    long long total_unmatch=0;
-    long long total_mis=0;
-
+    /************需要注释掉********* */ 
+    // long long total_match=0;
+    // long long total_unmatch=0;
+    // long long total_mis=0;
+    /************需要注释掉********* */
+    int groupLen = 0;
+    int return_sites[GROUP_SIZE];
+	unsigned int return_sites_error[GROUP_SIZE];
 
     ///这些是整个window被完全覆盖的
     for (i = 0; i < dumy->length; i++)
@@ -383,54 +413,136 @@ char* r_string)
         ///不能超过y的剩余长度
         o_len = MIN(Window_Len, currentIDLen - y_start);
         
-        recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
-                R_INF, overlap_list->list[currentID].y_id);
-
-        x_string = r_string + x_start;
-        y_string = dumy->overlap_region;
-
         
-
-
         ///这个长度足够，可以用来一起比
         if (Window_Len == o_len)
         {
-            
-            end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
 
+            recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+                R_INF, overlap_list->list[currentID].y_id);
 
-            if (error!=(unsigned int)-1)
+            x_string = r_string + x_start;
+
+            groupLen++;
+            if (groupLen == GROUP_SIZE)
             {
-                total_match++;
-            }
-            else
-            {
-                total_unmatch++;
-            }
-
             
-            test_edit_distance_by_edlib(x_string, y_string, x_len, 
-            o_len, THRESHOLD, error, &total_mis);
+                Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
+                dumy->overlap_region_group[2], dumy->overlap_region_group[3], o_len, x_string, x_len,
+			        return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
+                groupLen = 0;
+
+
+                // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error);
+                // if (error!=return_sites_error[0])
+                // {
+                //     fprintf(stderr, "error\n");
+                // }
+                // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[1], o_len, x_string, x_len, THRESHOLD, &error); 
+                // if (error!=return_sites_error[1])
+                // {
+                //     fprintf(stderr, "error\n");
+                // }
+                // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[2], o_len, x_string, x_len, THRESHOLD, &error);
+                // if (error!=return_sites_error[2])
+                // {
+                //     fprintf(stderr, "error\n");
+                // } 
+                // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[3], o_len, x_string, x_len, THRESHOLD, &error);
+                // if (error!=return_sites_error[3])
+                // {
+                //     fprintf(stderr, "error\n");
+                // }  
+                
+               
+            //    /************需要注释掉**********/
+            //    for (size_t ijk = 0; ijk < GROUP_SIZE; ijk++)
+            //    {
+            //        if (return_sites_error[ijk]!=(unsigned int)-1)
+            //         {
+            //             total_match++;
+            //         }
+            //         else
+            //         {
+            //             total_unmatch++;
+            //         }
+            //         test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, 
+            //         return_sites_error[ijk], &total_mis);
+            //    }
+            //     /************需要注释掉********* */
+                
+            }
 
         }
         else  ///这个不够，只能单个比 ///不够的地方要置N
         {
-            /**
+
+            recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+                R_INF, overlap_list->list[currentID].y_id);
+
+            x_string = r_string + x_start;
+            y_string = dumy->overlap_region;
+
+            ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
+            memset (y_string + o_len, 0, Window_Len - o_len);
             end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
-
-
-            if (error!=-1)
-            {
-                total_match++;
-            }
-            else
-            {
-                total_unmatch++;
-            }
-            **/
+            
+            // /************需要注释掉**********/
+            // if (error!=(unsigned int)-1)
+            // {
+            //     total_match++;
+            // }
+            // else
+            // {
+            //     total_unmatch++;
+            // }
+            // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
+            // /************需要注释掉********* */
         }
     }
 
+
+
+    if (groupLen == 1)
+    {
+        end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error);
+        // /************需要注释掉**********/
+        // if (error!=(unsigned int)-1)
+        // {
+        //     total_match++;
+        // }
+        // else
+        // {
+        //     total_unmatch++;
+        // }
+        // test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[0], x_len, o_len, THRESHOLD, error, &total_mis);
+        // /************需要注释掉********* */
+    }
+    else if (groupLen > 1)
+    {
+        Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
+                dumy->overlap_region_group[2], dumy->overlap_region_group[3], o_len, x_string, x_len,
+			        return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
+        // /************需要注释掉**********/
+        // for (size_t ijk = 0; ijk < groupLen; ijk++)
+        // {
+        //     if (return_sites_error[ijk]!=(unsigned int)-1)
+        //     {
+        //         total_match++;
+        //     }
+        //     else
+        //     {
+        //         total_unmatch++;
+        //     }
+        //     test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, 
+        //     return_sites_error[ijk], &total_mis);
+        // }
+        // /************需要注释掉**********/
+        groupLen = 0;
+    }
+    
+    
+
     long long reverse_i = dumy->size - 1;
     
     ///这些是整个window被部分覆盖的
@@ -471,41 +583,38 @@ char* r_string)
         y_string = dumy->overlap_region;
 
         ///不够的地方要置N
-
+        if (Window_Len != o_len)
+        {
+            ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
+            memset (y_string + o_len, 0, Window_Len - o_len);
+        }
         end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
 
-        /**
-        if (error!=-1)
-        {
-            total_match++;
-        }
-        else
-        {
-            total_unmatch++;
-        }
-        **/
+
+        // /************需要注释掉********* */
+        // if (error!=(unsigned int)-1)
+        // {
+        //     total_match++;
+        // }
+        // else
+        // {
+        //     total_unmatch++;
+        // }
+        // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
+        // /************需要注释掉********* */
     }
     
-    /**
-    if (total_match!=0)
-    {
+
+
+
+    // /************需要注释掉********* */
+    // pthread_mutex_lock(&debug_statistics);
+    // T_total_match = T_total_match + total_match;
+    // T_total_unmatch = T_total_unmatch + total_unmatch;
+    // T_total_mis = T_total_mis + total_mis;
+	// pthread_mutex_unlock(&debug_statistics);
+    // /************需要注释掉********* */
         
-        fprintf(stderr, "total_match: %u\n", total_match);
-        fprintf(stderr, "total_unmatch: %u\n", total_unmatch);
-        fprintf(stderr, "total_mis: %u\n", total_mis);
-        
-    }
-    **/
-
-    T_total_match = T_total_match + total_match;
-    T_total_unmatch = T_total_unmatch + total_unmatch;
-    T_total_mis = T_total_mis + total_mis;
-    
-    
-
-
-    
-    
 }
 
 void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy)
@@ -522,90 +631,6 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re
 
     window_start = 0;
     window_end = WINDOW - 1;
-
-
-
-    /**
-    UC_Read debug_read;
-    init_UC_Read(&debug_read);
-
-    for (i = 0; i < overlap_list->length; i++)
-    {
-        
-        if (overlap_list->list[i].y_pos_strand)
-        {
-            recover_UC_Read_RC(&debug_read, R_INF, overlap_list->list[i].y_id);
-        }
-        else
-        {
-            recover_UC_Read(&debug_read, R_INF, overlap_list->list[i].y_id);
-        }
-        
-        
-        long long x_length = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1;
-        long long y_length = overlap_list->list[i].y_pos_e - overlap_list->list[i].y_pos_s + 1;
-
-        EdlibAlignResult result = edlibAlign(g_read->seq+overlap_list->list[i].x_pos_s, 
-        x_length, 
-        debug_read.seq + overlap_list->list[i].y_pos_s, 
-        y_length, 
-            edlibNewAlignConfig(-1, EDLIB_MODE_NW, EDLIB_TASK_DISTANCE, NULL, 0));
-            
-
-        if (result.status == EDLIB_STATUS_OK) {
-            if (result.editDistance>0 && result.editDistance < x_length*0.02)
-            {
-                fprintf(stderr, "inner_i: %u, x_length: %u, y_length: %u, editDistance: %u\n", 
-                i, x_length, y_length, result.editDistance);
-                fprintf(stderr, "x_pos_s: %u, x_pos_e: %u\n", 
-                overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e);
-                fprintf(stderr, "y_pos_s: %u, y_pos_e: %u, y_pos_strand: %u\n", 
-                overlap_list->list[i].y_pos_s, overlap_list->list[i].y_pos_e, overlap_list->list[i].y_pos_strand);
-
-
-                EdlibAlignResult n_result = edlibAlign(g_read->seq+overlap_list->list[i].x_pos_s, 
-                350, 
-                debug_read.seq + overlap_list->list[i].y_pos_s - 15, 
-                380, 
-                    edlibNewAlignConfig(-1, EDLIB_MODE_HW, EDLIB_TASK_DISTANCE, NULL, 0));
-                
-                fprintf(stderr, "n_editDistance: %u\n", 
-                n_result.editDistance);
-
-
-                
-
-                edlibFreeAlignResult(n_result);
-            }
-            
-            
-        }
-        edlibFreeAlignResult(result);
-        
-
-    }
-
-    destory_UC_Read(&debug_read);
-    **/
-    ///return;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
     int flag;
 
     for (i = 0; i < window_num; i++)
@@ -641,18 +666,16 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re
         {
             window_end = g_read->length - 1;
         }
-
-        ///break;
-
     }
 
+
     
-    fprintf(stderr, "total_match: %u, total_unmatch: %u, total_mis: %u\n", 
-    T_total_match, T_total_unmatch, T_total_mis);
-    
-
-
-
+    // /************需要注释掉********* */
+    // pthread_mutex_lock(&debug_statistics);
+    // fprintf(stderr, "total_match: %u, total_unmatch: %u, total_mis: %u\n", 
+    // T_total_match, T_total_unmatch, T_total_mis);
+	// pthread_mutex_unlock(&debug_statistics);
+    // /************需要注释掉********* */
 }
 
 
@@ -663,6 +686,11 @@ void init_Correct_dumy(Correct_dumy* list)
     list->lengthNT = 0;
     list->start_i = 0;
     list->overlapID = NULL;
+    int i;
+    for (i = 0; i < 256; i++)
+	{
+		list->Peq_SSE[i] = _mm_setzero_si128();
+	}
 }
 
 void destory_Correct_dumy(Correct_dumy* list)
diff --git a/Correct.h b/Correct.h
index 3990ac3..225a4f6 100644
--- a/Correct.h
+++ b/Correct.h
@@ -2,9 +2,11 @@
 #define __CORRECT__
 #include <stdint.h>
 #include "Hash_Table.h"
+#include "Levenshtein_distance.h"
 
 #define WINDOW 350
 #define THRESHOLD  15
+#define GROUP_SIZE 4
 
 typedef struct
 {
@@ -14,6 +16,8 @@ typedef struct
     uint64_t size;
     uint64_t start_i;
     char overlap_region[WINDOW + THRESHOLD*2 + 10];
+    char overlap_region_group[GROUP_SIZE][WINDOW + THRESHOLD*2 + 10];
+    __m128i Peq_SSE[256];
 } Correct_dumy;
 
 
diff --git a/Hash_Table.cpp b/Hash_Table.cpp
index 02fda01..b13bfc2 100644
--- a/Hash_Table.cpp
+++ b/Hash_Table.cpp
@@ -3,6 +3,7 @@
 #include <stdint.h>
 #include "Hash_Table.h"
 #include "Process_Read.h"
+#include "Correct.h"
 #include <pthread.h>
 pthread_mutex_t output_mutex;
 
@@ -715,6 +716,19 @@ inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType*
     candidates->length++;
 }
 
+void pre_filter_by_nearby(k_mer_pos* n_list, uint64_t n_length, uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read)
+{
+    char* x_string = NULL;
+    char* y_string = NULL;
+    
+    long long i = 0;
+    for (i = 0; i < n_length; i++)
+    {
+        /* code */
+    }
+    
+}
+
 
 void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand)
 {
diff --git a/Levenshtein_distance.cpp b/Levenshtein_distance.cpp
index 4622c9d..6468c31 100644
--- a/Levenshtein_distance.cpp
+++ b/Levenshtein_distance.cpp
@@ -1 +1,806 @@
 #include "Levenshtein_distance.h"
+
+void output_bit_myers(Word x, int length)
+{
+	int i = 0;
+	while (i < length)
+	{
+		fprintf(stderr, "%u", (x >> i) & ((Word)1));
+		i++;
+	}
+	fprintf(stderr, "\n");
+}
+
+void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i)
+{
+
+	i++;
+	int k = 0;
+	int j = i;
+	int diff;
+	while (k < length)
+	{
+		
+		int x_p = (VP >> k) & ((Word)1);
+		int x_n = (VN >> k) & ((Word)1);
+		if (x_p == 1 && x_n == 1)
+		{
+			fprintf(stderr, "error\n");
+		}
+		
+		
+		
+		if (x_p == 1)
+		{
+			diff = 1;
+			///fprintf(stderr, "[+1]");
+		}
+		
+		if (x_n == 1)
+		{
+			diff = -1;
+			///fprintf(stderr, "[-1]");
+		}
+
+		if (x_p == 0 && x_n == 0)
+		{
+			diff = 0;
+			///fprintf(stderr, "[+0]");
+		}
+		j++;
+		if (matrix[i][j] - matrix[i][j - 1] != diff)
+		{
+			fprintf(stderr, "*************V(k): %u\n", k);
+			///return;
+		}
+		k++;
+	}
+	///fprintf(stderr, "\n");
+}
+
+void prase_D0(Word D0, int length, int matrix[1000][1000], int i)
+{
+
+	i++;
+	int k = 0;
+	int j = i;
+	int diff;
+	while (k < length)
+	{
+		
+		int diff = (D0 >> k) & ((Word)1);
+
+		if(diff==matrix[i][j] - matrix[i-1][j-1])
+		{
+			fprintf(stderr, "*************D(k): %u\n", k);
+			fprintf(stderr, "diff: %u, matrix[i][j]: %u, matrix[i-1][j-1]: %u\n", diff, matrix[i][j], matrix[i-1][j-1]);
+			///return;
+		}
+		j++;
+		k++;
+	}
+	///fprintf(stderr, "\n");
+}
+
+void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i)
+{
+
+	i++;
+	int k = 0;
+	int j = i;
+	int diff;
+	while (k < length)
+	{
+		
+		int x_p = (HP >> k) & ((Word)1);
+		int x_n = (HN >> k) & ((Word)1);
+		if (x_p == 1 && x_n == 1)
+		{
+			fprintf(stderr, "error\n");
+		}
+		if (x_p == 1)
+		{
+			diff = 1;
+			///fprintf(stderr, "[+1]");
+		}
+		if (x_n == 1)
+		{
+			diff = -1;
+			///fprintf(stderr, "[-1]");
+		}
+		if (x_p == 0 && x_n == 0)
+		{
+			diff = 0;
+			///fprintf(stderr, "[+0]");
+		}
+
+		if(diff!=matrix[i][j] - matrix[i-1][j])
+		{
+			fprintf(stderr, "*************H(k): %u\n", k);
+			///return;
+		}
+
+		j++;
+		k++;
+	}
+	///fprintf(stderr, "\n");
+}
+
+/**
+ pattern是长的那个，是y
+ p_length是长的那个的长度,  p_length实际没用
+ text是短的那个，是x
+ t_length是短的那个的长度
+ errthold是阈值
+ return_err是编辑距离
+ 返回值是结束位置
+ **/
+int Reserve_Banded_BPM_debug
+(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err, int matrix[1000][1000])
+{
+	(*return_err) = (unsigned int)-1;
+
+	Word Peq[256];
+
+	int band_length = (errthold << 1) + 1;
+	int i = 0;
+	Word tmp_Peq_1 = (Word)1;
+
+	Peq['A'] = (Word)0;
+	Peq['T'] = (Word)0;
+	Peq['G'] = (Word)0;
+	Peq['C'] = (Word)0;
+
+
+	Word Peq_A;
+	Word Peq_T;
+	Word Peq_C;
+	Word Peq_G;
+
+	///band_length = 2k + 1
+	for (i = 0; i<band_length; i++)
+	{
+		Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
+		tmp_Peq_1 = tmp_Peq_1 << 1;
+	}
+
+	///Peq['T'] = Peq['T'] | Peq['C'];
+
+	Peq_A = Peq['A'];
+	Peq_C = Peq['C'];
+	Peq_T = Peq['T'];
+	Peq_G = Peq['G'];
+
+
+	memset(Peq, 0, sizeof(Word)* 256);
+
+
+	Peq['A'] = Peq_A;
+	Peq['C'] = Peq_C;
+	Peq['T'] = Peq_T;
+	Peq['G'] = Peq_G;
+
+
+	
+
+	Word Mask = ((Word)1 << (errthold << 1));
+
+	Word VP = 0;
+	Word VN = 0;
+	Word X = 0;
+	Word D0 = 0;
+	Word HN = 0;
+	Word HP = 0;
+
+
+	i = 0;
+
+	int err = 0;
+
+	Word err_mask = (Word)1;
+
+
+	///band_down = 2k
+	///i_bd = 2k
+	///int i_bd = i + band_down;
+	int i_bd = (errthold << 1);
+
+
+	int last_high = (errthold << 1);
+
+
+	/// t_length_1 = SEQ_LENGTH - 1
+	int t_length_1 = t_length - 1;
+	//while(i<t_length)
+
+	while (i<t_length_1)
+	{
+		fprintf(stderr, "i: %u, j_s: %u, j_e: %u\n", i + 1, i + 1, i + band_length + 1);
+		///if (i >= 6)
+		if (i >= 0)
+		{
+			/**
+			fprintf(stderr, "VP:\n");
+			output_bit_myers(VP, band_length);
+
+			fprintf(stderr, "VN:\n");
+			output_bit_myers(VN, band_length);
+
+			fprintf(stderr, "HP:\n");
+			output_bit_myers(HP, band_length);
+
+			fprintf(stderr, "HN:\n");
+			output_bit_myers(HN, band_length);
+
+			fprintf(stderr, "D0:\n");
+			output_bit_myers(D0, band_length);
+
+			fprintf(stderr, "text[i]: %c\n", text[i]);
+
+			fprintf(stderr, "Peq[text[i]]:\n");
+			output_bit_myers(Peq[text[i]], band_length);
+
+			for (size_t j = i; j < i + band_length; j++)
+			{
+				fprintf(stderr, "%c", pattern[j]);
+			}
+			fprintf(stderr, "\n");
+			
+			fprintf(stderr, "Previous begin.\n");
+			prase_vertical(VP, VN, band_length, matrix, i-1);
+			prase_D0(D0, band_length, matrix, i-1);
+			prase_H(HP, HN, band_length, matrix, i-1);
+			fprintf(stderr, "Previous test done.\n");
+			**/
+
+			X = Peq[text[i]] | VN;
+			/**
+			fprintf(stderr, "#X:\n");
+			output_bit_myers(X, band_length);
+			**/
+
+			D0 = ((VP + (X&VP)) ^ VP) | X;
+			/**
+			fprintf(stderr, "#(X&VP):\n");
+			output_bit_myers((X&VP), band_length);
+
+			fprintf(stderr, "#(VP + (X&VP)):\n");
+			output_bit_myers((VP + (X&VP)), band_length);
+
+			fprintf(stderr, "#((VP + (X&VP)) ^ VP):\n");
+			output_bit_myers(((VP + (X&VP)) ^ VP), band_length);
+
+			fprintf(stderr, "#D0:\n");
+			output_bit_myers(D0, band_length);
+			**/
+			HN = VP&D0;
+			HP = VN | ~(VP | D0);
+
+			X = D0 >> 1;
+			VN = X&HP;
+			VP = HN | ~(X | HP);
+		}
+		else
+		{
+			///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
+			X = Peq[text[i]] | VN;
+
+			D0 = ((VP + (X&VP)) ^ VP) | X;
+
+			HN = VP&D0;
+			HP = VN | ~(VP | D0);
+
+			X = D0 >> 1;
+			VN = X&HP;
+			VP = HN | ~(X | HP);
+		}
+		
+
+		
+
+		
+		
+		
+		/**
+		for (size_t j = i + 1; j <= i + band_length + 1; j++)
+		{
+			fprintf(stderr, "[%u]", matrix[i + 1][j]);
+		}
+		fprintf(stderr, "\n");
+		**/
+
+		prase_vertical(VP, VN, band_length, matrix, i);
+		prase_D0(D0, band_length, matrix, i);
+		prase_H(HP, HN, band_length, matrix, i);
+		/**
+		fprintf(stderr, "VP:\n");
+		output_bit_myers(VP, band_length);
+
+		fprintf(stderr, "VN:\n");
+		output_bit_myers(VN, band_length);
+		**/
+		
+
+		if (!(D0&err_mask))
+		{
+			++err;
+
+			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
+			if ((err - last_high)>errthold)
+			{
+				///fprintf(stderr, "0 ######, i: %u\n", i);
+				return -1;
+			}
+				
+		}
+
+
+		Peq['A'] = Peq['A'] >> 1;
+		Peq['C'] = Peq['C'] >> 1;
+		Peq['G'] = Peq['G'] >> 1;
+		Peq['T'] = Peq['T'] >> 1;
+
+
+		++i;
+		++i_bd;
+		Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
+
+
+		///Peq['T'] = Peq['T'] | Peq['C'];
+	}
+
+
+
+
+
+	X = Peq[text[i]] | VN;
+	D0 = ((VP + (X&VP)) ^ VP) | X;
+	HN = VP&D0;
+	HP = VN | ~(VP | D0);
+	X = D0 >> 1;
+	VN = X&HP;
+	VP = HN | ~(X | HP);
+	if (!(D0&err_mask))
+	{
+		++err;
+		if ((err - last_high)>errthold)
+			return -1;
+	}
+
+	prase_vertical(VP, VN, band_length, matrix, i);
+	prase_D0(D0, band_length, matrix, i);
+	prase_H(HP, HN, band_length, matrix, i);
+
+
+	fprintf(stderr, "err: %d, matrix[][]: %d\n", err, matrix[i+1][i+1]);
+	fprintf(stderr, "VP:\n");
+	output_bit_myers(VP, band_length);
+
+	fprintf(stderr, "VN:\n");
+	output_bit_myers(VN, band_length);
+
+	////fprintf(stderr, "sucess(2)\n");
+
+	/// last_high = 2k
+	/// site = (SEQ_LENGTH + 2k) - 2k -1
+	/// site = SEQ_LENGTH - 1
+	///int site = p_length - last_high - 1;
+	int site = t_length - 1;
+	int return_site = -1;
+	if ((err <= errthold) && (err<=*return_err))
+	{
+		*return_err = err;
+		return_site = site;
+	}
+	int i_last = i;
+	i = 0;
+
+
+
+
+	while (i<errthold)
+	{
+		err = err + ((VP >> i)&(Word)1);
+		err = err - ((VN >> i)&(Word)1);
+		++i;
+
+		fprintf(stderr, "*i: %u, err: %d\n", i, err);
+
+		if ((err <= errthold) && (err <= *return_err))
+		{
+			*return_err = err;
+			return_site = site + i;
+		}
+	}
+
+
+	unsigned int ungap_err;
+	ungap_err = err;
+
+
+	while (i<last_high)
+	{
+		err = err + ((VP >> i)&(Word)1);
+		err = err - ((VN >> i)&(Word)1);
+		++i;
+
+		fprintf(stderr, "*i: %u, err: %d\n", i, err);
+
+		if ((err <= errthold) && (err<=*return_err))
+		{
+			*return_err = err;
+			return_site = site + i;
+		}
+
+
+
+
+	}
+
+
+	if ((ungap_err <= errthold) && (ungap_err == *return_err))
+	{
+		return_site = site + errthold;
+	}
+
+	return return_site;
+
+}
+
+int BS_Reserve_Banded_BPM
+(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err)
+{
+	(*return_err) = (unsigned int)-1;
+
+	///Õâ¸öÊÇÄǸöÐèÒªÔ¤´¦ÀíµÄÏòÁ¿
+	Word Peq[256];
+
+	int band_length = (errthold << 1) + 1;
+	int i = 0;
+	Word tmp_Peq_1 = (Word)1;
+
+	Peq['A'] = (Word)0;
+	Peq['T'] = (Word)0;
+	Peq['G'] = (Word)0;
+	Peq['C'] = (Word)0;
+
+
+	Word Peq_A;
+	Word Peq_T;
+	Word Peq_C;
+	Word Peq_G;
+
+	///band_length = 2k + 1
+	///ÕâÊÇ°ÑpatternµÄÇ°2k + 1¸ö×Ö·ûÔ¤´¦Àí
+	///pattern[0]¶ÔÓ¦Peq[0]
+	///pattern[2k]¶ÔÓ¦Peq[2k]
+	for (i = 0; i<band_length; i++)
+	{
+		Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
+		tmp_Peq_1 = tmp_Peq_1 << 1;
+	}
+
+	Peq['T'] = Peq['T'] | Peq['C'];
+
+	Peq_A = Peq['A'];
+	Peq_C = Peq['C'];
+	Peq_T = Peq['T'];
+	Peq_G = Peq['G'];
+
+
+	memset(Peq, 0, sizeof(Word)* 256);
+
+
+	Peq['A'] = Peq_A;
+	Peq['C'] = Peq_C;
+	Peq['T'] = Peq_T;
+	Peq['G'] = Peq_G;
+
+
+	
+
+	Word Mask = ((Word)1 << (errthold << 1));
+
+	Word VP = 0;
+	Word VN = 0;
+	Word X = 0;
+	Word D0 = 0;
+	Word HN = 0;
+	Word HP = 0;
+
+
+	i = 0;
+
+	int err = 0;
+
+	Word err_mask = (Word)1;
+
+
+	///band_down = 2k
+	///i_bd = 2k
+	///int i_bd = i + band_down;
+	int i_bd = (errthold << 1);
+
+
+	int last_high = (errthold << 1);
+
+
+	/// t_length_1 = SEQ_LENGTH - 1
+	int t_length_1 = t_length - 1;
+	//while(i<t_length)
+
+	while (i<t_length_1)
+	{
+		///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
+		X = Peq[text[i]] | VN;
+
+		D0 = ((VP + (X&VP)) ^ VP) | X;
+
+		HN = VP&D0;
+		HP = VN | ~(VP | D0);
+
+		X = D0 >> 1;
+		VN = X&HP;
+		VP = HN | ~(X | HP);
+		///Èç¹ûб¶Ô½ÇÏß·½ÏòÆ¥ÅäÔòD0ÊÇ1
+		///Èç¹û²»Æ¥ÅäÔòD0ÊÇ0
+		///Õâ¸öÒâ˼ÊÇÈç¹û×îÉÏÃæÄÇÌõ¶Ô½ÇÏßÉϵÄб¶Ô½ÇÏß·½Ïò·¢ÉúÎóÅä,ÔòÖ´ÐÐÄÚ²¿³ÌÐò
+		///
+		if (!(D0&err_mask))
+		{
+			++err;
+
+			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
+			if ((err - last_high)>errthold)
+				return -1;
+		}
+
+		///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
+		//ÓÒÒÆʵ¼ÊÉÏÊÇ°Ñpattern[0]ÒƵôÁË
+		Peq['A'] = Peq['A'] >> 1;
+		Peq['C'] = Peq['C'] >> 1;
+		Peq['G'] = Peq['G'] >> 1;
+		Peq['T'] = Peq['T'] >> 1;
+
+
+		++i;
+		++i_bd;
+		///ÕâÊÇ°ÑеÄpattern[2k]¼Ó½øÀ´, ÕâòËÆÊǼӵ½Peq[2k]ÉÏÁË
+		Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
+
+
+		Peq['T'] = Peq['T'] | Peq['C'];
+	}
+
+
+
+
+	///fprintf(stderr, "sucess(1)\n");
+
+
+	///Õâ¸öÑ­»·ÄóöÀ´ÊÇΪÁË·ÀÖ¹ÄÚ´æй¶
+	///ÆäʵҲ¾ÍÊÇÑ­»·ÀïµÄ×îºóÒ»ÐÐÓï¾ä°É
+	///ÍêÈ«¿ÉÒÔ°ÑpatternÔö´óһλ
+	///²»¹ýÕâÑùÒ²ºÃ£¬¿ÉÒÔ¼õÉÙ¼ÆË㿪Ïú
+	X = Peq[text[i]] | VN;
+	D0 = ((VP + (X&VP)) ^ VP) | X;
+	HN = VP&D0;
+	HP = VN | ~(VP | D0);
+	X = D0 >> 1;
+	VN = X&HP;
+	VP = HN | ~(X | HP);
+	if (!(D0&err_mask))
+	{
+		++err;
+		if ((err - last_high)>errthold)
+			return -1;
+	}
+
+
+
+
+
+	////fprintf(stderr, "sucess(2)\n");
+
+	/// last_high = 2k
+	/// site = (SEQ_LENGTH + 2k) - 2k -1
+	/// site = SEQ_LENGTH - 1
+	///´ËʱÕâ¸ösiteòËÆÊÇ×îÉÏÃæÄÇÌõ¶Ô½ÇÏßµÄλÖÃ
+	///int site = p_length - last_high - 1;
+	int site = t_length - 1;
+	int return_site = -1;
+	if ((err <= errthold) && (err<=*return_err))
+	{
+		*return_err = err;
+		return_site = site;
+	}
+	int i_last = i;
+	i = 0;
+
+
+
+
+	while (i<errthold)
+	{
+		err = err + ((VP >> i)&(Word)1);
+		err = err - ((VN >> i)&(Word)1);
+		++i;
+
+		if ((err <= errthold) && (err <= *return_err))
+		{
+			*return_err = err;
+			return_site = site + i;
+		}
+	}
+
+
+	unsigned int ungap_err;
+	ungap_err = err;
+
+
+	while (i<last_high)
+	{
+		err = err + ((VP >> i)&(Word)1);
+		err = err - ((VN >> i)&(Word)1);
+		++i;
+
+		if ((err <= errthold) && (err<=*return_err))
+		{
+			*return_err = err;
+			return_site = site + i;
+		}
+
+
+
+
+	}
+
+
+	if ((ungap_err <= errthold) && (ungap_err == *return_err))
+	{
+		return_site = site + errthold;
+	}
+
+	return return_site;
+
+}
+
+int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length,
+                unsigned short errthold,unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id)
+{
+
+    Word Peq[128];
+    char Peq_index[4]= {'A','C','G','T'};
+    int symbol = 0;
+    int r;
+    Word tmp_Peq_1=(Word)1;
+
+
+    Peq['A']=(Word)0;
+    Peq['T']=(Word)0;
+    Peq['G']=(Word)0;
+    Peq['C']=(Word)0;
+    Word Peq_A;
+    Word Peq_T;
+    Word Peq_C;
+    Word Peq_G;
+
+    for (r =0; r<band_length; r++)
+    {
+            Peq[pattern[r]]=Peq[pattern[r]]|tmp_Peq_1;
+        tmp_Peq_1=tmp_Peq_1<<1;
+    }
+    Peq_A=Peq['A'];
+    Peq_C=Peq['C'];
+    Peq_T=Peq['T'];
+    Peq_G=Peq['G'];
+
+    for(symbol = 0; symbol < 128; symbol++)
+    {
+        Peq[symbol]=(Word)0;
+        //jump_c[symbol]=128;
+    }
+    Peq['A']=Peq_A;
+    Peq['C']=Peq_C;
+    Peq['T']=Peq_T;
+    Peq['G']=Peq_G;
+
+
+    Word Mask_Pre=(Word)1<<(band_length-2);
+    Word Mask=(Word)1<<(band_length-1);
+    Word VP=0;
+    Word VN=0;
+    Word X=0;
+    Word D0=0;
+    Word HN=0;
+    Word HP=0;
+
+
+    int s=0;
+    int i = 0;
+    int j=0;
+
+    int bound=band_length-2-band_down;
+    int err=0;
+
+    Word err_mask=(Word)1;
+    int s1=band_length-2;
+    int i_bd=i+band_down;
+    int last_high=band_length-t_length+p_length-band_down-1;
+   int t_length_1=t_length-1;
+    //while(i<t_length)
+    while(i<t_length_1)
+    {
+
+        X=Peq[text[i]]|VN;
+
+        D0=((VP+(X&VP))^VP)|X;
+
+        HN=VP&D0;
+        HP=VN|~(VP|D0);
+
+        X=D0>>1;
+        VN=X&HP;
+        VP=HN|~(X|HP);
+        if(!(D0&err_mask))
+         {
+            ++err;
+            if((err-last_high)>errthold)
+                return -1;
+        }
+
+        Peq['A']=Peq['A']>>1;
+        Peq['C']=Peq['C']>>1;
+        Peq['G']=Peq['G']>>1;
+        Peq['T']=Peq['T']>>1;
+
+
+        ++i;
+        ++i_bd;
+        Peq[pattern[i_bd]]=Peq[pattern[i_bd]]|Mask;
+    }
+
+
+    ///这个循环拿出来是为了防止内存泄露
+
+       X=Peq[text[i]]|VN;
+        D0=((VP+(X&VP))^VP)|X;
+        HN=VP&D0;
+        HP=VN|~(VP|D0);
+        X=D0>>1;
+        VN=X&HP;
+        VP=HN|~(X|HP);
+        if(!(D0&err_mask))
+         {
+            ++err;
+            if((err-last_high)>errthold)
+                return -1;
+        }
+
+
+    int site=p_length-last_high-1;
+    int return_site=-1;
+    if((err<=errthold)&&(err<*return_err))
+    {
+            *return_err=err;
+            return_site=site;
+    }
+    int i_last=i;
+    i=0;
+    while(i<last_high)
+    {
+        err=err+((VP>>i)&(Word)1);
+        err=err-((VN>>i)&(Word)1);
+        ++i;
+
+        if((err<=errthold)&&(err<*return_err))
+        {
+                *return_err=err;
+                return_site=site+i;
+        }
+
+
+    }
+    return return_site;
+
+}
diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h
index cd8d239..2bdc58a 100644
--- a/Levenshtein_distance.h
+++ b/Levenshtein_distance.h
@@ -10,462 +10,11 @@
 #include <stdio.h>
 
 typedef uint64_t Word;
+typedef uint32_t Word_32;
 
 
 
 
-inline void output_bit_myers(Word x, int length)
-{
-	int i = 0;
-	while (i < length)
-	{
-		fprintf(stderr, "%u", (x >> i) & ((Word)1));
-		i++;
-	}
-	fprintf(stderr, "\n");
-}
-
-
-inline void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i)
-{
-
-	i++;
-	int k = 0;
-	int j = i;
-	int diff;
-	while (k < length)
-	{
-		
-		int x_p = (VP >> k) & ((Word)1);
-		int x_n = (VN >> k) & ((Word)1);
-		if (x_p == 1 && x_n == 1)
-		{
-			fprintf(stderr, "error\n");
-		}
-		
-		
-		
-		if (x_p == 1)
-		{
-			diff = 1;
-			///fprintf(stderr, "[+1]");
-		}
-		
-		if (x_n == 1)
-		{
-			diff = -1;
-			///fprintf(stderr, "[-1]");
-		}
-
-		if (x_p == 0 && x_n == 0)
-		{
-			diff = 0;
-			///fprintf(stderr, "[+0]");
-		}
-		j++;
-		if (matrix[i][j] - matrix[i][j - 1] != diff)
-		{
-			fprintf(stderr, "*************V(k): %u\n", k);
-			///return;
-		}
-		k++;
-	}
-	///fprintf(stderr, "\n");
-}
-
-
-
-
-inline void prase_D0(Word D0, int length, int matrix[1000][1000], int i)
-{
-
-	i++;
-	int k = 0;
-	int j = i;
-	int diff;
-	while (k < length)
-	{
-		
-		int diff = (D0 >> k) & ((Word)1);
-
-		if(diff==matrix[i][j] - matrix[i-1][j-1])
-		{
-			fprintf(stderr, "*************D(k): %u\n", k);
-			fprintf(stderr, "diff: %u, matrix[i][j]: %u, matrix[i-1][j-1]: %u\n", diff, matrix[i][j], matrix[i-1][j-1]);
-			///return;
-		}
-		j++;
-		k++;
-	}
-	///fprintf(stderr, "\n");
-}
-
-inline void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i)
-{
-
-	i++;
-	int k = 0;
-	int j = i;
-	int diff;
-	while (k < length)
-	{
-		
-		int x_p = (HP >> k) & ((Word)1);
-		int x_n = (HN >> k) & ((Word)1);
-		if (x_p == 1 && x_n == 1)
-		{
-			fprintf(stderr, "error\n");
-		}
-		if (x_p == 1)
-		{
-			diff = 1;
-			///fprintf(stderr, "[+1]");
-		}
-		if (x_n == 1)
-		{
-			diff = -1;
-			///fprintf(stderr, "[-1]");
-		}
-		if (x_p == 0 && x_n == 0)
-		{
-			diff = 0;
-			///fprintf(stderr, "[+0]");
-		}
-
-		if(diff!=matrix[i][j] - matrix[i-1][j])
-		{
-			fprintf(stderr, "*************H(k): %u\n", k);
-			///return;
-		}
-
-		j++;
-		k++;
-	}
-	///fprintf(stderr, "\n");
-}
-
-
-/**
- pattern是长的那个，是y
- p_length是长的那个的长度,  p_length实际没用
- text是短的那个，是x
- t_length是短的那个的长度
- errthold是阈值
- return_err是编辑距离
- 返回值是结束位置
- **/
-inline int Reserve_Banded_BPM_debug
-(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err, int matrix[1000][1000])
-{
-	(*return_err) = (unsigned int)-1;
-
-	Word Peq[256];
-
-	int band_length = (errthold << 1) + 1;
-	int i = 0;
-	Word tmp_Peq_1 = (Word)1;
-
-	Peq['A'] = (Word)0;
-	Peq['T'] = (Word)0;
-	Peq['G'] = (Word)0;
-	Peq['C'] = (Word)0;
-
-
-	Word Peq_A;
-	Word Peq_T;
-	Word Peq_C;
-	Word Peq_G;
-
-	///band_length = 2k + 1
-	for (i = 0; i<band_length; i++)
-	{
-		Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
-		tmp_Peq_1 = tmp_Peq_1 << 1;
-	}
-
-	///Peq['T'] = Peq['T'] | Peq['C'];
-
-	Peq_A = Peq['A'];
-	Peq_C = Peq['C'];
-	Peq_T = Peq['T'];
-	Peq_G = Peq['G'];
-
-
-	memset(Peq, 0, sizeof(Word)* 256);
-
-
-	Peq['A'] = Peq_A;
-	Peq['C'] = Peq_C;
-	Peq['T'] = Peq_T;
-	Peq['G'] = Peq_G;
-
-
-	
-
-	Word Mask = ((Word)1 << (errthold << 1));
-
-	Word VP = 0;
-	Word VN = 0;
-	Word X = 0;
-	Word D0 = 0;
-	Word HN = 0;
-	Word HP = 0;
-
-
-	i = 0;
-
-	int err = 0;
-
-	Word err_mask = (Word)1;
-
-
-	///band_down = 2k
-	///i_bd = 2k
-	///int i_bd = i + band_down;
-	int i_bd = (errthold << 1);
-
-
-	int last_high = (errthold << 1);
-
-
-	/// t_length_1 = SEQ_LENGTH - 1
-	int t_length_1 = t_length - 1;
-	//while(i<t_length)
-
-	while (i<t_length_1)
-	{
-		fprintf(stderr, "i: %u, j_s: %u, j_e: %u\n", i + 1, i + 1, i + band_length + 1);
-		///if (i >= 6)
-		if (i >= 0)
-		{
-			/**
-			fprintf(stderr, "VP:\n");
-			output_bit_myers(VP, band_length);
-
-			fprintf(stderr, "VN:\n");
-			output_bit_myers(VN, band_length);
-
-			fprintf(stderr, "HP:\n");
-			output_bit_myers(HP, band_length);
-
-			fprintf(stderr, "HN:\n");
-			output_bit_myers(HN, band_length);
-
-			fprintf(stderr, "D0:\n");
-			output_bit_myers(D0, band_length);
-
-			fprintf(stderr, "text[i]: %c\n", text[i]);
-
-			fprintf(stderr, "Peq[text[i]]:\n");
-			output_bit_myers(Peq[text[i]], band_length);
-
-			for (size_t j = i; j < i + band_length; j++)
-			{
-				fprintf(stderr, "%c", pattern[j]);
-			}
-			fprintf(stderr, "\n");
-			
-			fprintf(stderr, "Previous begin.\n");
-			prase_vertical(VP, VN, band_length, matrix, i-1);
-			prase_D0(D0, band_length, matrix, i-1);
-			prase_H(HP, HN, band_length, matrix, i-1);
-			fprintf(stderr, "Previous test done.\n");
-			**/
-
-			X = Peq[text[i]] | VN;
-			/**
-			fprintf(stderr, "#X:\n");
-			output_bit_myers(X, band_length);
-			**/
-
-			D0 = ((VP + (X&VP)) ^ VP) | X;
-			/**
-			fprintf(stderr, "#(X&VP):\n");
-			output_bit_myers((X&VP), band_length);
-
-			fprintf(stderr, "#(VP + (X&VP)):\n");
-			output_bit_myers((VP + (X&VP)), band_length);
-
-			fprintf(stderr, "#((VP + (X&VP)) ^ VP):\n");
-			output_bit_myers(((VP + (X&VP)) ^ VP), band_length);
-
-			fprintf(stderr, "#D0:\n");
-			output_bit_myers(D0, band_length);
-			**/
-			HN = VP&D0;
-			HP = VN | ~(VP | D0);
-
-			X = D0 >> 1;
-			VN = X&HP;
-			VP = HN | ~(X | HP);
-		}
-		else
-		{
-			///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
-			X = Peq[text[i]] | VN;
-
-			D0 = ((VP + (X&VP)) ^ VP) | X;
-
-			HN = VP&D0;
-			HP = VN | ~(VP | D0);
-
-			X = D0 >> 1;
-			VN = X&HP;
-			VP = HN | ~(X | HP);
-		}
-		
-
-		
-
-		
-		
-		
-		/**
-		for (size_t j = i + 1; j <= i + band_length + 1; j++)
-		{
-			fprintf(stderr, "[%u]", matrix[i + 1][j]);
-		}
-		fprintf(stderr, "\n");
-		**/
-
-		prase_vertical(VP, VN, band_length, matrix, i);
-		prase_D0(D0, band_length, matrix, i);
-		prase_H(HP, HN, band_length, matrix, i);
-		/**
-		fprintf(stderr, "VP:\n");
-		output_bit_myers(VP, band_length);
-
-		fprintf(stderr, "VN:\n");
-		output_bit_myers(VN, band_length);
-		**/
-		
-
-		if (!(D0&err_mask))
-		{
-			++err;
-
-			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
-			if ((err - last_high)>errthold)
-			{
-				///fprintf(stderr, "0 ######, i: %u\n", i);
-				return -1;
-			}
-				
-		}
-
-
-		Peq['A'] = Peq['A'] >> 1;
-		Peq['C'] = Peq['C'] >> 1;
-		Peq['G'] = Peq['G'] >> 1;
-		Peq['T'] = Peq['T'] >> 1;
-
-
-		++i;
-		++i_bd;
-		Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
-
-
-		///Peq['T'] = Peq['T'] | Peq['C'];
-	}
-
-
-
-
-
-	X = Peq[text[i]] | VN;
-	D0 = ((VP + (X&VP)) ^ VP) | X;
-	HN = VP&D0;
-	HP = VN | ~(VP | D0);
-	X = D0 >> 1;
-	VN = X&HP;
-	VP = HN | ~(X | HP);
-	if (!(D0&err_mask))
-	{
-		++err;
-		if ((err - last_high)>errthold)
-			return -1;
-	}
-
-	prase_vertical(VP, VN, band_length, matrix, i);
-	prase_D0(D0, band_length, matrix, i);
-	prase_H(HP, HN, band_length, matrix, i);
-
-
-	fprintf(stderr, "err: %d, matrix[][]: %d\n", err, matrix[i+1][i+1]);
-	fprintf(stderr, "VP:\n");
-	output_bit_myers(VP, band_length);
-
-	fprintf(stderr, "VN:\n");
-	output_bit_myers(VN, band_length);
-
-	////fprintf(stderr, "sucess(2)\n");
-
-	/// last_high = 2k
-	/// site = (SEQ_LENGTH + 2k) - 2k -1
-	/// site = SEQ_LENGTH - 1
-	///int site = p_length - last_high - 1;
-	int site = t_length - 1;
-	int return_site = -1;
-	if ((err <= errthold) && (err<=*return_err))
-	{
-		*return_err = err;
-		return_site = site;
-	}
-	int i_last = i;
-	i = 0;
-
-
-
-
-	while (i<errthold)
-	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
-		++i;
-
-		fprintf(stderr, "*i: %u, err: %d\n", i, err);
-
-		if ((err <= errthold) && (err <= *return_err))
-		{
-			*return_err = err;
-			return_site = site + i;
-		}
-	}
-
-
-	unsigned int ungap_err;
-	ungap_err = err;
-
-
-	while (i<last_high)
-	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
-		++i;
-
-		fprintf(stderr, "*i: %u, err: %d\n", i, err);
-
-		if ((err <= errthold) && (err<=*return_err))
-		{
-			*return_err = err;
-			return_site = site + i;
-		}
-
-
-
-
-	}
-
-
-	if ((ungap_err <= errthold) && (ungap_err == *return_err))
-	{
-		return_site = site + errthold;
-	}
-
-	return return_site;
-
-}
-
 
 /**
  pattern是长的那个，是y
@@ -679,372 +228,362 @@ inline int Reserve_Banded_BPM
 
 }
 
+inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *pattern3, char *pattern4, int p_length, char *text, int t_length,
+	int* return_sites, unsigned int* return_sites_error, unsigned short errthold, __m128i* Peq_SSE)
 
-
-
-
-
-
-
-inline int BS_Reserve_Banded_BPM
-(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err)
 {
-	(*return_err) = (unsigned int)-1;
-
-	///Õâ¸öÊÇÄǸöÐèÒªÔ¤´¦ÀíµÄÏòÁ¿
-	Word Peq[256];
+	memset(return_sites, -1, sizeof(int)* 4);
+	memset(return_sites_error, -1, sizeof(unsigned int)* 4);
 
+	Word_32 Peq[256][4];
 	int band_length = (errthold << 1) + 1;
-	int i = 0;
-	Word tmp_Peq_1 = (Word)1;
-
-	Peq['A'] = (Word)0;
-	Peq['T'] = (Word)0;
-	Peq['G'] = (Word)0;
-	Peq['C'] = (Word)0;
 
 
-	Word Peq_A;
-	Word Peq_T;
-	Word Peq_C;
-	Word Peq_G;
+	int i;
+
+	Word_32 tmp_Peq_1 = 1;
+
+
+	memset(Peq['A'], 0, sizeof(Word_32)* 4);
+	memset(Peq['C'], 0, sizeof(Word_32)* 4);
+	memset(Peq['G'], 0, sizeof(Word_32)* 4);
+	memset(Peq['T'], 0, sizeof(Word_32)* 4);
 
-	///band_length = 2k + 1
-	///ÕâÊÇ°ÑpatternµÄÇ°2k + 1¸ö×Ö·ûÔ¤´¦Àí
-	///pattern[0]¶ÔÓ¦Peq[0]
-	///pattern[2k]¶ÔÓ¦Peq[2k]
 	for (i = 0; i<band_length; i++)
 	{
-		Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
+		Peq[pattern1[i]][0] = Peq[pattern1[i]][0] | tmp_Peq_1;
+		Peq[pattern2[i]][1] = Peq[pattern2[i]][1] | tmp_Peq_1;
+		Peq[pattern3[i]][2] = Peq[pattern3[i]][2] | tmp_Peq_1;
+		Peq[pattern4[i]][3] = Peq[pattern4[i]][3] | tmp_Peq_1;
+
 		tmp_Peq_1 = tmp_Peq_1 << 1;
 	}
 
-	Peq['T'] = Peq['T'] | Peq['C'];
+	Peq_SSE['A'] = _mm_set_epi32(Peq['A'][3], Peq['A'][2], Peq['A'][1], Peq['A'][0]);
+	Peq_SSE['C'] = _mm_set_epi32(Peq['C'][3], Peq['C'][2], Peq['C'][1], Peq['C'][0]);
+	Peq_SSE['G'] = _mm_set_epi32(Peq['G'][3], Peq['G'][2], Peq['G'][1], Peq['G'][0]);
+	Peq_SSE['T'] = _mm_set_epi32(Peq['T'][3], Peq['T'][2], Peq['T'][1], Peq['T'][0]);
+	///Peq_SSE['T'] = _mm_or_si128(Peq_SSE['T'], Peq_SSE['C']);
 
-	Peq_A = Peq['A'];
-	Peq_C = Peq['C'];
-	Peq_T = Peq['T'];
-	Peq_G = Peq['G'];
+	Word_32 Mask = ((Word_32)1 << (errthold << 1));
+
+	__m128i Mask1 = _mm_set_epi32(0, 0, 0, Mask);
+	__m128i Mask2 = _mm_set_epi32(0, 0, Mask, 0);
+	__m128i Mask3 = _mm_set_epi32(0, Mask, 0, 0);
+	__m128i Mask4 = _mm_set_epi32(Mask, 0, 0, 0);
 
 
-	memset(Peq, 0, sizeof(Word)* 256);
+	__m128i VP = _mm_setzero_si128();
+	__m128i VN = _mm_setzero_si128();
+	__m128i X = _mm_setzero_si128();
+	__m128i D0 = _mm_setzero_si128();
+	__m128i HN = _mm_setzero_si128();
+	__m128i HP = _mm_setzero_si128();
+	__m128i tmp_process;
+	__m128i tmp_process1;
 
 
-	Peq['A'] = Peq_A;
-	Peq['C'] = Peq_C;
-	Peq['T'] = Peq_T;
-	Peq['G'] = Peq_G;
+
+	__m128i Err_4 = _mm_setzero_si128();
+	__m128i err_mask = _mm_set_epi32(1, 1, 1, 1);
+	///for_not li quan shi 1
+	__m128i for_not = _mm_set1_epi32(-1);
+	__m128i err_arry;
+	__m128i cmp_result;
 
 
-	
-
-	Word Mask = ((Word)1 << (errthold << 1));
-
-	Word VP = 0;
-	Word VN = 0;
-	Word X = 0;
-	Word D0 = 0;
-	Word HN = 0;
-	Word HP = 0;
-
-
-	i = 0;
-
-	int err = 0;
-
-	Word err_mask = (Word)1;
-
-
-	///band_down = 2k
-	///i_bd = 2k
-	///int i_bd = i + band_down;
 	int i_bd = (errthold << 1);
-
-
 	int last_high = (errthold << 1);
-
-
-	/// t_length_1 = SEQ_LENGTH - 1
 	int t_length_1 = t_length - 1;
-	//while(i<t_length)
+	int err1;
+	int err2;
+	int err3;
+	int err4;
+
+	__m128i pre_end = _mm_set_epi32
+		(last_high + errthold, last_high + errthold, last_high + errthold, last_high + errthold);
+	i = 0;
 
 	while (i<t_length_1)
 	{
-		///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
-		X = Peq[text[i]] | VN;
+		///X = Peq[text[i]] | VN;
+		X = _mm_or_si128(Peq_SSE[text[i]], VN);
 
-		D0 = ((VP + (X&VP)) ^ VP) | X;
 
-		HN = VP&D0;
-		HP = VN | ~(VP | D0);
 
-		X = D0 >> 1;
-		VN = X&HP;
-		VP = HN | ~(X | HP);
-		///Èç¹ûб¶Ô½ÇÏß·½ÏòÆ¥ÅäÔòD0ÊÇ1
-		///Èç¹û²»Æ¥ÅäÔòD0ÊÇ0
-		///Õâ¸öÒâ˼ÊÇÈç¹û×îÉÏÃæÄÇÌõ¶Ô½ÇÏßÉϵÄб¶Ô½ÇÏß·½Ïò·¢ÉúÎóÅä,ÔòÖ´ÐÐÄÚ²¿³ÌÐò
-		///
-		if (!(D0&err_mask))
-		{
-			++err;
+		/*************D0 = ((VP + (X&VP)) ^ VP) | X*********************/
+		///X&VP
+		tmp_process1 = _mm_and_si128(X, VP);
+		///(VP + (X&VP))
+		tmp_process = _mm_add_epi32(tmp_process1, VP);
+		///((VP + (X&VP)) ^ VP)
+		tmp_process = _mm_xor_si128(tmp_process, VP);
+		///((VP + (X&VP)) ^ VP) | X
+		D0 = _mm_or_si128(tmp_process, X);
+		/*************D0 = ((VP + (X&VP)) ^ VP) | X*********************/
 
-			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
-			if ((err - last_high)>errthold)
-				return -1;
-		}
+		///HN = VP&D0;
+		HN = _mm_and_si128(D0, VP);
 
-		///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
-		//ÓÒÒÆʵ¼ÊÉÏÊÇ°Ñpattern[0]ÒƵôÁË
-		Peq['A'] = Peq['A'] >> 1;
-		Peq['C'] = Peq['C'] >> 1;
-		Peq['G'] = Peq['G'] >> 1;
-		Peq['T'] = Peq['T'] >> 1;
+		///HP = VN | ~(VP | D0);
+		tmp_process = _mm_or_si128(D0, VP);
+		tmp_process = _mm_andnot_si128(tmp_process, for_not);
+		HP = _mm_or_si128(tmp_process, VN);
 
+		///X = D0 >> 1;
+		X = _mm_srli_epi32(D0, 1);
+		///VN = X&HP;
+		VN = _mm_and_si128(X, HP);
+		///VP = HN | ~(X | HP);
+		tmp_process = _mm_or_si128(X, HP);
+		tmp_process = _mm_andnot_si128(tmp_process, for_not);
+		VP = _mm_or_si128(HN, tmp_process);
+
+		///D0&err_mask
+		err_arry = _mm_and_si128(D0, err_mask);
+		Err_4 = _mm_add_epi32(Err_4, err_mask);
+		Err_4 = _mm_sub_epi32(Err_4, err_arry);
+
+		/**************** */
+		///shi ji shang zhe ge zhi hen xiao d
+		cmp_result = _mm_cmpgt_epi32(Err_4, pre_end);
+
+		///jian zhi
+		if (_mm_extract_epi32(cmp_result, 0) && _mm_extract_epi32(cmp_result, 1)
+			&& _mm_extract_epi32(cmp_result, 2) && _mm_extract_epi32(cmp_result, 3))
+			return 1;
+		/**************** */
+
+
+		Peq_SSE['A'] = _mm_srli_epi32(Peq_SSE['A'], 1);
+		Peq_SSE['T'] = _mm_srli_epi32(Peq_SSE['T'], 1);
+		Peq_SSE['G'] = _mm_srli_epi32(Peq_SSE['G'], 1);
+		Peq_SSE['C'] = _mm_srli_epi32(Peq_SSE['C'], 1);
 
 		++i;
 		++i_bd;
-		///ÕâÊÇ°ÑеÄpattern[2k]¼Ó½øÀ´, ÕâòËÆÊǼӵ½Peq[2k]ÉÏÁË
-		Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
 
-
-		Peq['T'] = Peq['T'] | Peq['C'];
+		Peq_SSE[pattern1[i_bd]] = _mm_or_si128(Mask1, Peq_SSE[pattern1[i_bd]]);
+		Peq_SSE[pattern2[i_bd]] = _mm_or_si128(Mask2, Peq_SSE[pattern2[i_bd]]);
+		Peq_SSE[pattern3[i_bd]] = _mm_or_si128(Mask3, Peq_SSE[pattern3[i_bd]]);
+		Peq_SSE[pattern4[i_bd]] = _mm_or_si128(Mask4, Peq_SSE[pattern4[i_bd]]);
+		///Peq_SSE['T'] = _mm_or_si128(Peq_SSE['T'], Peq_SSE['C']);
 	}
 
 
 
+		///X = Peq[text[i]] | VN;
+	X = _mm_or_si128(Peq_SSE[text[i]], VN);
 
-	///fprintf(stderr, "sucess(1)\n");
+	/*************D0 = ((VP + (X&VP)) ^ VP) | X*********************/
+	///X&VP
+	tmp_process1 = _mm_and_si128(X, VP);
+	///(VP + (X&VP))
+	tmp_process = _mm_add_epi32(tmp_process1, VP);
+	///((VP + (X&VP)) ^ VP)
+	tmp_process = _mm_xor_si128(tmp_process, VP);
+	///((VP + (X&VP)) ^ VP) | X
+	D0 = _mm_or_si128(tmp_process, X);
+	/*************D0 = ((VP + (X&VP)) ^ VP) | X*********************/
+
+	///HN = VP&D0;
+	HN = _mm_and_si128(D0, VP);
+
+	///HP = VN | ~(VP | D0);
+	tmp_process = _mm_or_si128(D0, VP);
+	tmp_process = _mm_andnot_si128(tmp_process, for_not);
+	HP = _mm_or_si128(tmp_process, VN);
 
 
-	///Õâ¸öÑ­»·ÄóöÀ´ÊÇΪÁË·ÀÖ¹ÄÚ´æй¶
-	///ÆäʵҲ¾ÍÊÇÑ­»·ÀïµÄ×îºóÒ»ÐÐÓï¾ä°É
-	///ÍêÈ«¿ÉÒÔ°ÑpatternÔö´óһλ
-	///²»¹ýÕâÑùÒ²ºÃ£¬¿ÉÒÔ¼õÉÙ¼ÆË㿪Ïú
-	X = Peq[text[i]] | VN;
-	D0 = ((VP + (X&VP)) ^ VP) | X;
-	HN = VP&D0;
-	HP = VN | ~(VP | D0);
-	X = D0 >> 1;
-	VN = X&HP;
-	VP = HN | ~(X | HP);
-	if (!(D0&err_mask))
-	{
-		++err;
-		if ((err - last_high)>errthold)
-			return -1;
-	}
+	///X = D0 >> 1;
+	X = _mm_srli_epi32(D0, 1);
+	///VN = X&HP;
+	VN = _mm_and_si128(X, HP);
+	///VP = HN | ~(X | HP);
+	tmp_process = _mm_or_si128(X, HP);
+	tmp_process = _mm_andnot_si128(tmp_process, for_not);
+	VP = _mm_or_si128(HN, tmp_process);
 
+	///D0&err_mask
+	err_arry = _mm_and_si128(D0, err_mask);
+	Err_4 = _mm_add_epi32(Err_4, err_mask);
+	Err_4 = _mm_sub_epi32(Err_4, err_arry);
 
+	///shi ji shang zhe ge zhi hen xiao d
+	cmp_result = _mm_cmpgt_epi32(Err_4, pre_end);
 
+	///jian zhi
+	if (_mm_extract_epi32(cmp_result, 0) && _mm_extract_epi32(cmp_result, 1)
+		&& _mm_extract_epi32(cmp_result, 2) && _mm_extract_epi32(cmp_result, 3))
+		return 1;
 
-
-	////fprintf(stderr, "sucess(2)\n");
-
-	/// last_high = 2k
-	/// site = (SEQ_LENGTH + 2k) - 2k -1
-	/// site = SEQ_LENGTH - 1
-	///´ËʱÕâ¸ösiteòËÆÊÇ×îÉÏÃæÄÇÌõ¶Ô½ÇÏßµÄλÖÃ
-	///int site = p_length - last_high - 1;
 	int site = t_length - 1;
-	int return_site = -1;
-	if ((err <= errthold) && (err<=*return_err))
-	{
-		*return_err = err;
-		return_site = site;
-	}
-	int i_last = i;
-	i = 0;
+	err1 = _mm_extract_epi32(Err_4, 0);
+	err2 = _mm_extract_epi32(Err_4, 1);
+	err3 = _mm_extract_epi32(Err_4, 2);
+	err4 = _mm_extract_epi32(Err_4, 3);
 
 
+	if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
+	{
+		return_sites[0] = site;
+		return_sites_error[0] = err1;
+	}
+	if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
+	{
+		return_sites[1] = site;
+		return_sites_error[1] = err2;
+	}
+	if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
+	{
+		return_sites[2] = site;
+		return_sites_error[2] = err3;
+	}
+	if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
+	{
+		return_sites[3] = site;
+		return_sites_error[3] = err4;
+	}
+
+
+	i = 0;
+
 
 
 	while (i<errthold)
 	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
+		///err = err + ((VP >> i)&(Word_32)1);
+		tmp_process = _mm_srli_epi32(VP, i);
+		tmp_process = _mm_and_si128(tmp_process, err_mask);
+		Err_4 = _mm_add_epi32(Err_4, tmp_process);
+
+		///err = err - ((VN >> i)&(Word_32)1);
+		tmp_process1 = _mm_srli_epi32(VN, i);
+		tmp_process1 = _mm_and_si128(tmp_process1, err_mask);
+		Err_4 = _mm_sub_epi32(Err_4, tmp_process1);
 		++i;
 
-		if ((err <= errthold) && (err <= *return_err))
+		err1 = _mm_extract_epi32(Err_4, 0);
+		err2 = _mm_extract_epi32(Err_4, 1);
+		err3 = _mm_extract_epi32(Err_4, 2);
+		err4 = _mm_extract_epi32(Err_4, 3);
+
+
+		if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
 		{
-			*return_err = err;
-			return_site = site + i;
+			return_sites[0] = site + i;
+			return_sites_error[0] = err1;
+		}
+		if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
+		{
+			return_sites[1] = site + i;
+			return_sites_error[1] = err2;
+		}
+		if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
+		{
+			return_sites[2] = site + i;
+			return_sites_error[2] = err3;
+		}
+		if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
+		{
+			return_sites[3] = site + i;
+			return_sites_error[3] = err4;
 		}
 	}
 
 
-	unsigned int ungap_err;
-	ungap_err = err;
+	unsigned int ungap_err1;
+	unsigned int ungap_err2;
+	unsigned int ungap_err3;
+	unsigned int ungap_err4;
+	ungap_err1 = err1;
+	ungap_err2 = err2;
+	ungap_err3 = err3;
+	ungap_err4 = err4;
+
+
 
 
 	while (i<last_high)
 	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
+		///err = err + ((VP >> i)&(Word_32)1);
+		tmp_process = _mm_srli_epi32(VP, i);
+		tmp_process = _mm_and_si128(tmp_process, err_mask);
+		Err_4 = _mm_add_epi32(Err_4, tmp_process);
+
+		///err = err - ((VN >> i)&(Word_32)1);
+		tmp_process1 = _mm_srli_epi32(VN, i);
+		tmp_process1 = _mm_and_si128(tmp_process1, err_mask);
+		Err_4 = _mm_sub_epi32(Err_4, tmp_process1);
 		++i;
 
-		if ((err <= errthold) && (err<=*return_err))
+		err1 = _mm_extract_epi32(Err_4, 0);
+		err2 = _mm_extract_epi32(Err_4, 1);
+		err3 = _mm_extract_epi32(Err_4, 2);
+		err4 = _mm_extract_epi32(Err_4, 3);
+
+
+		if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
 		{
-			*return_err = err;
-			return_site = site + i;
+			return_sites[0] = site + i;
+			return_sites_error[0] = err1;
+		}
+		if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
+		{
+			return_sites[1] = site + i;
+			return_sites_error[1] = err2;
+		}
+		if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
+		{
+			return_sites[2] = site + i;
+			return_sites_error[2] = err3;
+		}
+		if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
+		{
+			return_sites[3] = site + i;
+			return_sites_error[3] = err4;
 		}
-
-
-
-
 	}
 
 
-	if ((ungap_err <= errthold) && (ungap_err == *return_err))
+	if ((ungap_err1 <= errthold) && ungap_err1 == return_sites_error[0])
 	{
-		return_site = site + errthold;
+		return_sites[0] = site + errthold;
 	}
 
-	return return_site;
+
+	if ((ungap_err2 <= errthold) && ungap_err2 == return_sites_error[1])
+	{
+		return_sites[1] = site + errthold;
+	}
+
+
+	if ((ungap_err3 <= errthold) && ungap_err3 == return_sites_error[2])
+	{
+		return_sites[2] = site + errthold;
+	}
+
+
+	if ((ungap_err4 <= errthold) && ungap_err4 == return_sites_error[3])
+	{
+		return_sites[3] = site + errthold;
+	}
+
+	return 1;
+	
 
 }
 
 
-
-inline int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length,
-                unsigned short errthold,unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id)
-{
-
-    Word Peq[128];
-    char Peq_index[4]= {'A','C','G','T'};
-    int symbol = 0;
-    int r;
-    Word tmp_Peq_1=(Word)1;
-
-
-    Peq['A']=(Word)0;
-    Peq['T']=(Word)0;
-    Peq['G']=(Word)0;
-    Peq['C']=(Word)0;
-    Word Peq_A;
-    Word Peq_T;
-    Word Peq_C;
-    Word Peq_G;
-
-    for (r =0; r<band_length; r++)
-    {
-            Peq[pattern[r]]=Peq[pattern[r]]|tmp_Peq_1;
-        tmp_Peq_1=tmp_Peq_1<<1;
-    }
-    Peq_A=Peq['A'];
-    Peq_C=Peq['C'];
-    Peq_T=Peq['T'];
-    Peq_G=Peq['G'];
-
-    for(symbol = 0; symbol < 128; symbol++)
-    {
-        Peq[symbol]=(Word)0;
-        //jump_c[symbol]=128;
-    }
-    Peq['A']=Peq_A;
-    Peq['C']=Peq_C;
-    Peq['T']=Peq_T;
-    Peq['G']=Peq_G;
-
-
-    Word Mask_Pre=(Word)1<<(band_length-2);
-    Word Mask=(Word)1<<(band_length-1);
-    Word VP=0;
-    Word VN=0;
-    Word X=0;
-    Word D0=0;
-    Word HN=0;
-    Word HP=0;
-
-
-    int s=0;
-    int i = 0;
-    int j=0;
-
-    int bound=band_length-2-band_down;
-    int err=0;
-
-    Word err_mask=(Word)1;
-    int s1=band_length-2;
-    int i_bd=i+band_down;
-    int last_high=band_length-t_length+p_length-band_down-1;
-   int t_length_1=t_length-1;
-    //while(i<t_length)
-    while(i<t_length_1)
-    {
-
-        X=Peq[text[i]]|VN;
-
-        D0=((VP+(X&VP))^VP)|X;
-
-        HN=VP&D0;
-        HP=VN|~(VP|D0);
-
-        X=D0>>1;
-        VN=X&HP;
-        VP=HN|~(X|HP);
-        if(!(D0&err_mask))
-         {
-            ++err;
-            if((err-last_high)>errthold)
-                return -1;
-        }
-
-        Peq['A']=Peq['A']>>1;
-        Peq['C']=Peq['C']>>1;
-        Peq['G']=Peq['G']>>1;
-        Peq['T']=Peq['T']>>1;
-
-
-        ++i;
-        ++i_bd;
-        Peq[pattern[i_bd]]=Peq[pattern[i_bd]]|Mask;
-    }
-
-
-    ///这个循环拿出来是为了防止内存泄露
-
-       X=Peq[text[i]]|VN;
-        D0=((VP+(X&VP))^VP)|X;
-        HN=VP&D0;
-        HP=VN|~(VP|D0);
-        X=D0>>1;
-        VN=X&HP;
-        VP=HN|~(X|HP);
-        if(!(D0&err_mask))
-         {
-            ++err;
-            if((err-last_high)>errthold)
-                return -1;
-        }
-
-
-    int site=p_length-last_high-1;
-    int return_site=-1;
-    if((err<=errthold)&&(err<*return_err))
-    {
-            *return_err=err;
-            return_site=site;
-    }
-    int i_last=i;
-    i=0;
-    while(i<last_high)
-    {
-        err=err+((VP>>i)&(Word)1);
-        err=err-((VN>>i)&(Word)1);
-        ++i;
-
-        if((err<=errthold)&&(err<*return_err))
-        {
-                *return_err=err;
-                return_site=site+i;
-        }
-
-
-    }
-    return return_site;
-
-}
-
+void output_bit_myers(Word x, int length);
+void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i);
+void prase_D0(Word D0, int length, int matrix[1000][1000], int i);
+void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i);
+int Reserve_Banded_BPM_debug(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, 
+unsigned int* return_err, int matrix[1000][1000]);
+int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length,unsigned short errthold,
+unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id);
+int BS_Reserve_Banded_BPM
+(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err);
 
 #endif
\ No newline at end of file