From cc0b4bed17cb388a24c5ecdcc54ee71cb09e58fd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?P=C3=A4r=20Larsson?= <par.larsson.medbio@gmail.com>
Date: Tue, 3 Mar 2026 20:13:50 +0100
Subject: [PATCH 1/5] start state

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index 237ede12..86b247ab 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 build/
 target/
+test_data/
+

From 6ef161e0fb8f0373d6b3a1c4cfcc1006864bf60d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?P=C3=A4r=20Larsson?= <par.larsson.medbio@gmail.com>
Date: Tue, 3 Mar 2026 20:55:44 +0100
Subject: [PATCH 2/5] avx2 for sw

---
 .gitattributes         |   2 +
 .gitignore             |   3 +
 CMakeLists.txt         |   5 +-
 cpp/ext/ssw/ssw.c      | 117 +++++++++--
 cpp/ext/ssw/ssw_avx2.c | 443 +++++++++++++++++++++++++++++++++++++++++
 cpp/ext/ssw/ssw_avx2.h |  31 +++
 pixi.lock              | 384 +++++++++++++++++++++++++++++++++++
 pixi.toml              |  15 ++
 8 files changed, 982 insertions(+), 18 deletions(-)
 create mode 100644 .gitattributes
 create mode 100644 cpp/ext/ssw/ssw_avx2.c
 create mode 100644 cpp/ext/ssw/ssw_avx2.h
 create mode 100644 pixi.lock
 create mode 100644 pixi.toml

diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 00000000..997504b4
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,2 @@
+# SCM syntax highlighting & preventing 3-way merges
+pixi.lock merge=binary linguist-language=YAML linguist-generated=true -diff
diff --git a/.gitignore b/.gitignore
index 86b247ab..6a38abc6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,6 @@ build/
 target/
 test_data/
 
+# pixi environments
+.pixi/*
+!.pixi/config.toml
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9ac96c3c..6f0d62f4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -114,7 +114,10 @@ endif()
 target_include_directories(salib PUBLIC cpp/ cpp/ext/ ${PROJECT_BINARY_DIR})
 target_link_libraries(salib PUBLIC ZLIB::ZLIB Threads::Threads zstr::zstr ISAL)
 IF(ENABLE_AVX)
-  target_compile_options(salib PUBLIC "-mavx2")
+  target_sources(salib PRIVATE cpp/ext/ssw/ssw_avx2.c)
+  set_source_files_properties(cpp/ext/ssw/ssw_avx2.c
+    PROPERTIES COMPILE_FLAGS "-mavx2")
+  target_compile_definitions(salib PUBLIC SSW_AVX2)
 ENDIF()
 if (TRACE)
    target_compile_definitions(salib PUBLIC "TRACE")
diff --git a/cpp/ext/ssw/ssw.c b/cpp/ext/ssw/ssw.c
index e8599c6f..82325040 100644
--- a/cpp/ext/ssw/ssw.c
+++ b/cpp/ext/ssw/ssw.c
@@ -80,6 +80,10 @@
 #endif
 
 
+#ifdef SSW_AVX2
+#include "ssw_avx2.h"
+#endif
+
 #ifdef __GNUC__
 #define LIKELY(x) __builtin_expect((x),1)
 #define UNLIKELY(x) __builtin_expect((x),0)
@@ -115,6 +119,10 @@ typedef struct {
 struct _profile{
 	__m128i* profile_byte;	// 0: none
 	__m128i* profile_word;	// 0: none
+#ifdef SSW_AVX2
+	void* profile_byte_avx2;	// __m256i*, 0: none
+	void* profile_word_avx2;	// __m256i*, 0: none
+#endif
 	const int8_t* read;
 	const int8_t* mat;
 	int32_t readLen;
@@ -790,6 +798,10 @@ s_profile* ssw_init (const int8_t* read, const int32_t readLen, const int8_t* ma
 	s_profile* p = (s_profile*)calloc(1, sizeof(struct _profile));
 	p->profile_byte = 0;
 	p->profile_word = 0;
+#ifdef SSW_AVX2
+	p->profile_byte_avx2 = 0;
+	p->profile_word_avx2 = 0;
+#endif
 	p->bias = 0;
 
 	if (score_size == 0 || score_size == 2) {
@@ -799,9 +811,17 @@ s_profile* ssw_init (const int8_t* read, const int32_t readLen, const int8_t* ma
 		bias = abs(bias);
 
 		p->bias = bias;
+#ifdef SSW_AVX2
+		p->profile_byte_avx2 = qP_avx2_byte(read, mat, readLen, n, bias);
+#endif
 		p->profile_byte = qP_byte (read, mat, readLen, n, bias);
 	}
-	if (score_size == 1 || score_size == 2) p->profile_word = qP_word (read, mat, readLen, n);
+	if (score_size == 1 || score_size == 2) {
+#ifdef SSW_AVX2
+		p->profile_word_avx2 = qP_avx2_word(read, mat, readLen, n);
+#endif
+		p->profile_word = qP_word (read, mat, readLen, n);
+	}
 	p->read = read;
 	p->mat = mat;
 	p->readLen = readLen;
@@ -812,6 +832,10 @@ s_profile* ssw_init (const int8_t* read, const int32_t readLen, const int8_t* ma
 void init_destroy (s_profile* p) {
 	free(p->profile_byte);
 	free(p->profile_word);
+#ifdef SSW_AVX2
+	free(p->profile_byte_avx2);
+	free(p->profile_word_avx2);
+#endif
 	free(p);
 }
 
@@ -841,6 +865,40 @@ s_align* ssw_align (const s_profile* prof,
 	}
 
 	// Find the alignment scores and ending positions
+#ifdef SSW_AVX2
+	if (prof->profile_byte_avx2) {
+		alignment_end_avx2* bests_avx2 = sw_avx2_byte(ref, 0, refLen, readLen, weight_gapO, weight_gapE, (const __m256i*)prof->profile_byte_avx2, -1, prof->bias, maskLen);
+		if (prof->profile_word_avx2 && bests_avx2[0].score == 255) {
+			free(bests_avx2);
+			bests_avx2 = sw_avx2_word(ref, 0, refLen, readLen, weight_gapO, weight_gapE, (const __m256i*)prof->profile_word_avx2, -1, maskLen);
+			word = 1;
+		} else if (bests_avx2[0].score == 255) {
+			fprintf(stderr, "Please set 2 to the score_size parameter of the function ssw_init, otherwise the alignment results will be incorrect.\n");
+			free(r);
+			return NULL;
+		}
+		/* Copy AVX2 results into the alignment_end format expected below */
+		bests = (alignment_end*)calloc(2, sizeof(alignment_end));
+		bests[0].score = bests_avx2[0].score;
+		bests[0].ref = bests_avx2[0].ref;
+		bests[0].read = bests_avx2[0].read;
+		bests[1].score = bests_avx2[1].score;
+		bests[1].ref = bests_avx2[1].ref;
+		bests[1].read = bests_avx2[1].read;
+		free(bests_avx2);
+	} else if (prof->profile_word_avx2) {
+		alignment_end_avx2* bests_avx2 = sw_avx2_word(ref, 0, refLen, readLen, weight_gapO, weight_gapE, (const __m256i*)prof->profile_word_avx2, -1, maskLen);
+		bests = (alignment_end*)calloc(2, sizeof(alignment_end));
+		bests[0].score = bests_avx2[0].score;
+		bests[0].ref = bests_avx2[0].ref;
+		bests[0].read = bests_avx2[0].read;
+		bests[1].score = bests_avx2[1].score;
+		bests[1].ref = bests_avx2[1].ref;
+		bests[1].read = bests_avx2[1].read;
+		free(bests_avx2);
+		word = 1;
+	} else
+#endif
 	if (prof->profile_byte) {
 		bests = sw_sse2_byte(ref, 0, refLen, readLen, weight_gapO, weight_gapE, prof->profile_byte, -1, prof->bias, maskLen);
 		if (prof->profile_word && bests[0].score == 255) {
@@ -875,23 +933,48 @@ s_align* ssw_align (const s_profile* prof,
 
 	// Find the beginning position of the best alignment.
 	read_reverse = seq_reverse(prof->read, r->read_end1);
-	if (word == 0) {
-		vP = qP_byte(read_reverse, prof->mat, r->read_end1 + 1, prof->n, prof->bias);
-		bests_reverse = sw_sse2_byte(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP, r->score1, prof->bias, maskLen);
-	} else {
-		vP = qP_word(read_reverse, prof->mat, r->read_end1 + 1, prof->n);
-		bests_reverse = sw_sse2_word(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP, r->score1, maskLen);
+#ifdef SSW_AVX2
+	if (prof->profile_byte_avx2 || prof->profile_word_avx2) {
+		alignment_end_avx2* bests_reverse_avx2;
+		if (word == 0) {
+			__m256i* vP_avx2 = qP_avx2_byte(read_reverse, prof->mat, r->read_end1 + 1, prof->n, prof->bias);
+			bests_reverse_avx2 = sw_avx2_byte(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP_avx2, r->score1, prof->bias, maskLen);
+			free(vP_avx2);
+		} else {
+			__m256i* vP_avx2 = qP_avx2_word(read_reverse, prof->mat, r->read_end1 + 1, prof->n);
+			bests_reverse_avx2 = sw_avx2_word(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP_avx2, r->score1, maskLen);
+			free(vP_avx2);
+		}
+		free(read_reverse);
+		r->ref_begin1 = bests_reverse_avx2[0].ref;
+		r->read_begin1 = r->read_end1 - bests_reverse_avx2[0].read;
+
+		if (UNLIKELY(r->score1 > bests_reverse_avx2[0].score)) {
+			fprintf(stderr, "Warning: The alignment path of one pair of sequences may miss a small part. [ssw.c ssw_align]\n");
+			r->flag = 2;
+		}
+		free(bests_reverse_avx2);
+	} else
+#endif
+	{
+		if (word == 0) {
+			vP = qP_byte(read_reverse, prof->mat, r->read_end1 + 1, prof->n, prof->bias);
+			bests_reverse = sw_sse2_byte(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP, r->score1, prof->bias, maskLen);
+		} else {
+			vP = qP_word(read_reverse, prof->mat, r->read_end1 + 1, prof->n);
+			bests_reverse = sw_sse2_word(ref, 1, r->ref_end1 + 1, r->read_end1 + 1, weight_gapO, weight_gapE, vP, r->score1, maskLen);
+		}
+		free(vP);
+		free(read_reverse);
+		r->ref_begin1 = bests_reverse[0].ref;
+		r->read_begin1 = r->read_end1 - bests_reverse[0].read;
+
+		if (UNLIKELY(r->score1 > bests_reverse[0].score)) {
+			fprintf(stderr, "Warning: The alignment path of one pair of sequences may miss a small part. [ssw.c ssw_align]\n");
+			r->flag = 2;
+		}
+		free(bests_reverse);
 	}
-	free(vP);
-	free(read_reverse);
-	r->ref_begin1 = bests_reverse[0].ref;
-	r->read_begin1 = r->read_end1 - bests_reverse[0].read;
-
-    if (UNLIKELY(r->score1 > bests_reverse[0].score)) { // banded_sw result will miss a small part
-		fprintf(stderr, "Warning: The alignment path of one pair of sequences may miss a small part. [ssw.c ssw_align]\n");
-        r->flag = 2;  
-    }
-    free(bests_reverse);
 
 	if ((7&flag) == 0 || ((2&flag) != 0 && r->score1 < filters) || ((4&flag) != 0 && (r->ref_end1 - r->ref_begin1 > filterd || r->read_end1 - r->read_begin1 > filterd))) goto end;
 
diff --git a/cpp/ext/ssw/ssw_avx2.c b/cpp/ext/ssw/ssw_avx2.c
new file mode 100644
index 00000000..c9ed3cdd
--- /dev/null
+++ b/cpp/ext/ssw/ssw_avx2.c
@@ -0,0 +1,443 @@
+/*
+ *  ssw_avx2.c
+ *
+ *  AVX2 (256-bit) implementation of striped Smith-Waterman.
+ *  Ported from the SSE2 implementation in ssw.c.
+ *
+ *  This file must be compiled with -mavx2.
+ */
+
+#include "ssw_avx2.h"
+#include <stdlib.h>
+#include <string.h>
+
+#ifdef __GNUC__
+#define LIKELY(x) __builtin_expect((x),1)
+#define UNLIKELY(x) __builtin_expect((x),0)
+#else
+#define LIKELY(x) (x)
+#define UNLIKELY(x) (x)
+#endif
+
+/* Cross-lane byte shift helpers.
+ * AVX2 _mm256_slli_si256 shifts within each 128-bit lane independently,
+ * so we need explicit cross-lane carry for the striped pattern. */
+
+static inline __m256i avx2_shift_left_1_byte(__m256i v) {
+	/* Shift each 128-bit lane left by 1 byte (zeros inserted at byte 0 and byte 16) */
+	__m256i shifted = _mm256_slli_si256(v, 1);
+	/* Move lower 128-bit lane to upper position, zero the lower lane */
+	__m256i carry = _mm256_permute2x128_si256(v, v, 0x08);
+	/* Isolate byte[15] from (now in upper lane) */
+	carry = _mm256_srli_si256(carry, 15);
+	return _mm256_or_si256(shifted, carry);
+}
+
+static inline __m256i avx2_shift_left_2_bytes(__m256i v) {
+	__m256i shifted = _mm256_slli_si256(v, 2);
+	__m256i carry = _mm256_permute2x128_si256(v, v, 0x08);
+	carry = _mm256_srli_si256(carry, 14);
+	return _mm256_or_si256(shifted, carry);
+}
+
+/* Allocate 32-byte aligned memory, zero-initialized */
+static inline void* aligned_calloc(size_t count, size_t size) {
+	size_t total = count * size;
+	/* aligned_alloc requires size to be a multiple of alignment */
+	size_t aligned_total = (total + 31) & ~((size_t)31);
+	if (aligned_total == 0) aligned_total = 32;
+	void* p = aligned_alloc(32, aligned_total);
+	if (p) memset(p, 0, aligned_total);
+	return p;
+}
+
+/* Generate query profile for byte (8-bit) mode.
+ * 32 elements per __m256i vector. */
+__m256i* qP_avx2_byte(const int8_t* read_num,
+                       const int8_t* mat,
+                       const int32_t readLen,
+                       const int32_t n,
+                       uint8_t bias) {
+
+	int32_t segLen = (readLen + 31) / 32;
+	__m256i* vProfile = (__m256i*)aligned_calloc(n * segLen, sizeof(__m256i));
+	int8_t* t = (int8_t*)vProfile;
+	int32_t nt, i, j, segNum;
+
+	for (nt = 0; LIKELY(nt < n); nt++) {
+		for (i = 0; i < segLen; i++) {
+			j = i;
+			for (segNum = 0; LIKELY(segNum < 32); segNum++) {
+				*t++ = j >= readLen ? bias : mat[nt * n + read_num[j]] + bias;
+				j += segLen;
+			}
+		}
+	}
+	return vProfile;
+}
+
+/* Generate query profile for word (16-bit) mode.
+ * 16 elements per __m256i vector. */
+__m256i* qP_avx2_word(const int8_t* read_num,
+                       const int8_t* mat,
+                       const int32_t readLen,
+                       const int32_t n) {
+
+	int32_t segLen = (readLen + 15) / 16;
+	__m256i* vProfile = (__m256i*)aligned_calloc(n * segLen, sizeof(__m256i));
+	int16_t* t = (int16_t*)vProfile;
+	int32_t nt, i, j, segNum;
+
+	for (nt = 0; LIKELY(nt < n); nt++) {
+		for (i = 0; i < segLen; i++) {
+			j = i;
+			for (segNum = 0; LIKELY(segNum < 16); segNum++) {
+				*t++ = j >= readLen ? 0 : mat[nt * n + read_num[j]];
+				j += segLen;
+			}
+		}
+	}
+	return vProfile;
+}
+
+/* Striped Smith-Waterman, byte (8-bit) mode, AVX2.
+ * Processes 32 query positions in parallel per vector. */
+alignment_end_avx2* sw_avx2_byte(const int8_t* ref,
+                                  int8_t ref_dir,
+                                  int32_t refLen,
+                                  int32_t readLen,
+                                  const uint8_t weight_gapO,
+                                  const uint8_t weight_gapE,
+                                  const __m256i* vProfile,
+                                  uint8_t terminate,
+                                  uint8_t bias,
+                                  int32_t maskLen) {
+
+	/* Horizontal max: reduce 32 bytes in vm to scalar m */
+	#define max32(m, vm) \
+		(vm) = _mm256_max_epu8((vm), _mm256_permute2x128_si256((vm), (vm), 0x01)); \
+		(vm) = _mm256_max_epu8((vm), _mm256_srli_si256((vm), 8)); \
+		(vm) = _mm256_max_epu8((vm), _mm256_srli_si256((vm), 4)); \
+		(vm) = _mm256_max_epu8((vm), _mm256_srli_si256((vm), 2)); \
+		(vm) = _mm256_max_epu8((vm), _mm256_srli_si256((vm), 1)); \
+		(m) = (uint8_t)_mm256_extract_epi8((vm), 0)
+
+	uint8_t max = 0;
+	int32_t end_read = readLen - 1;
+	int32_t end_ref = -1;
+	int32_t segLen = (readLen + 31) / 32;
+
+	uint8_t* maxColumn = (uint8_t*)calloc(refLen, 1);
+	int32_t* end_read_column = (int32_t*)calloc(refLen, sizeof(int32_t));
+
+	__m256i vZero = _mm256_setzero_si256();
+
+	__m256i* pvHStore = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvHLoad  = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvE      = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvHmax   = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+
+	int32_t i, j, k;
+	__m256i vGapO = _mm256_set1_epi8(weight_gapO);
+	__m256i vGapE = _mm256_set1_epi8(weight_gapE);
+	__m256i vBias = _mm256_set1_epi8(bias);
+
+	__m256i vMaxScore = vZero;
+	__m256i vMaxMark = vZero;
+	__m256i vTemp;
+	int32_t edge, begin = 0, end = refLen, step = 1;
+
+	if (ref_dir == 1) {
+		begin = refLen - 1;
+		end = -1;
+		step = -1;
+	}
+	for (i = begin; LIKELY(i != end); i += step) {
+		int32_t cmp;
+		__m256i e, vF = vZero, vMaxColumn = vZero;
+
+		__m256i vH = pvHStore[segLen - 1];
+		vH = avx2_shift_left_1_byte(vH);
+		const __m256i* vP = vProfile + ref[i] * segLen;
+
+		__m256i* pv = pvHLoad;
+		pvHLoad = pvHStore;
+		pvHStore = pv;
+
+		/* Inner loop: process query segments */
+		for (j = 0; LIKELY(j < segLen); ++j) {
+			vH = _mm256_adds_epu8(vH, _mm256_load_si256(vP + j));
+			vH = _mm256_subs_epu8(vH, vBias);
+
+			e = _mm256_load_si256(pvE + j);
+			vH = _mm256_max_epu8(vH, e);
+			vH = _mm256_max_epu8(vH, vF);
+			vMaxColumn = _mm256_max_epu8(vMaxColumn, vH);
+
+			_mm256_store_si256(pvHStore + j, vH);
+
+			vH = _mm256_subs_epu8(vH, vGapO);
+			e = _mm256_subs_epu8(e, vGapE);
+			e = _mm256_max_epu8(e, vH);
+			_mm256_store_si256(pvE + j, e);
+
+			vF = _mm256_subs_epu8(vF, vGapE);
+			vF = _mm256_max_epu8(vF, vH);
+
+			vH = _mm256_load_si256(pvHLoad + j);
+		}
+
+		/* Lazy_F loop */
+		for (k = 0; LIKELY(k < 32); ++k) {
+			vF = avx2_shift_left_1_byte(vF);
+			for (j = 0; LIKELY(j < segLen); ++j) {
+				vH = _mm256_load_si256(pvHStore + j);
+				vH = _mm256_max_epu8(vH, vF);
+				vMaxColumn = _mm256_max_epu8(vMaxColumn, vH);
+				_mm256_store_si256(pvHStore + j, vH);
+				vH = _mm256_subs_epu8(vH, vGapO);
+				vF = _mm256_subs_epu8(vF, vGapE);
+				vTemp = _mm256_subs_epu8(vF, vH);
+				vTemp = _mm256_cmpeq_epi8(vTemp, vZero);
+				if (UNLIKELY(_mm256_movemask_epi8(vTemp) == (int32_t)0xffffffff)) goto end;
+			}
+		}
+
+end:
+		vMaxScore = _mm256_max_epu8(vMaxScore, vMaxColumn);
+		vTemp = _mm256_cmpeq_epi8(vMaxMark, vMaxScore);
+		cmp = _mm256_movemask_epi8(vTemp);
+		if (cmp != (int32_t)0xffffffff) {
+			uint8_t temp;
+			vMaxMark = vMaxScore;
+			max32(temp, vMaxScore);
+			vMaxScore = vMaxMark;
+
+			if (LIKELY(temp > max)) {
+				max = temp;
+				if (max + bias >= 255) break;
+				end_ref = i;
+				for (j = 0; LIKELY(j < segLen); ++j) pvHmax[j] = pvHStore[j];
+			}
+		}
+
+		/* Record the max score of current column. */
+		max32(maxColumn[i], vMaxColumn);
+		if (maxColumn[i] == terminate) break;
+	}
+
+	/* Trace the alignment ending position on read. */
+	uint8_t *t = (uint8_t*)pvHmax;
+	int32_t column_len = segLen * 32;
+	for (i = 0; LIKELY(i < column_len); ++i, ++t) {
+		int32_t temp;
+		if (*t == max) {
+			temp = i / 32 + i % 32 * segLen;
+			if (temp < end_read) end_read = temp;
+		}
+	}
+
+	free(pvHmax);
+	free(pvE);
+	free(pvHLoad);
+	free(pvHStore);
+
+	alignment_end_avx2* bests = (alignment_end_avx2*)calloc(2, sizeof(alignment_end_avx2));
+	bests[0].score = max + bias >= 255 ? 255 : max;
+	bests[0].ref = end_ref;
+	bests[0].read = end_read;
+
+	bests[1].score = 0;
+	bests[1].ref = 0;
+	bests[1].read = 0;
+
+	edge = (end_ref - maskLen) > 0 ? (end_ref - maskLen) : 0;
+	for (i = 0; i < edge; i++) {
+		if (maxColumn[i] > bests[1].score) {
+			bests[1].score = maxColumn[i];
+			bests[1].ref = i;
+		}
+	}
+	edge = (end_ref + maskLen) > refLen ? refLen : (end_ref + maskLen);
+	for (i = edge + 1; i < refLen; i++) {
+		if (maxColumn[i] > bests[1].score) {
+			bests[1].score = maxColumn[i];
+			bests[1].ref = i;
+		}
+	}
+
+	free(maxColumn);
+	free(end_read_column);
+	return bests;
+
+	#undef max32
+}
+
+/* Striped Smith-Waterman, word (16-bit) mode, AVX2.
+ * Processes 16 query positions in parallel per vector. */
+alignment_end_avx2* sw_avx2_word(const int8_t* ref,
+                                  int8_t ref_dir,
+                                  int32_t refLen,
+                                  int32_t readLen,
+                                  const uint8_t weight_gapO,
+                                  const uint8_t weight_gapE,
+                                  const __m256i* vProfile,
+                                  uint16_t terminate,
+                                  int32_t maskLen) {
+
+	/* Horizontal max: reduce 16 words in vm to scalar m */
+	#define max16_avx2(m, vm) \
+		(vm) = _mm256_max_epi16((vm), _mm256_permute2x128_si256((vm), (vm), 0x01)); \
+		(vm) = _mm256_max_epi16((vm), _mm256_srli_si256((vm), 8)); \
+		(vm) = _mm256_max_epi16((vm), _mm256_srli_si256((vm), 4)); \
+		(vm) = _mm256_max_epi16((vm), _mm256_srli_si256((vm), 2)); \
+		(m) = (uint16_t)_mm256_extract_epi16((vm), 0)
+
+	uint16_t max = 0;
+	int32_t end_read = readLen - 1;
+	int32_t end_ref = 0;
+	int32_t segLen = (readLen + 15) / 16;
+
+	uint16_t* maxColumn = (uint16_t*)calloc(refLen, sizeof(uint16_t));
+	int32_t* end_read_column = (int32_t*)calloc(refLen, sizeof(int32_t));
+
+	__m256i vZero = _mm256_setzero_si256();
+
+	__m256i* pvHStore = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvHLoad  = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvE      = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+	__m256i* pvHmax   = (__m256i*)aligned_calloc(segLen, sizeof(__m256i));
+
+	int32_t i, j, k;
+	__m256i vGapO = _mm256_set1_epi16(weight_gapO);
+	__m256i vGapE = _mm256_set1_epi16(weight_gapE);
+
+	__m256i vMaxScore = vZero;
+	__m256i vMaxMark = vZero;
+	__m256i vTemp;
+	int32_t edge, begin = 0, end = refLen, step = 1;
+
+	if (ref_dir == 1) {
+		begin = refLen - 1;
+		end = -1;
+		step = -1;
+	}
+	for (i = begin; LIKELY(i != end); i += step) {
+		int32_t cmp;
+		__m256i e, vF = vZero;
+
+		__m256i vH = pvHStore[segLen - 1];
+		vH = avx2_shift_left_2_bytes(vH);
+
+		__m256i* pv = pvHLoad;
+
+		__m256i vMaxColumn = vZero;
+
+		const __m256i* vP = vProfile + ref[i] * segLen;
+		pvHLoad = pvHStore;
+		pvHStore = pv;
+
+		/* Inner loop: process query segments */
+		for (j = 0; LIKELY(j < segLen); j++) {
+			vH = _mm256_adds_epi16(vH, _mm256_load_si256(vP + j));
+
+			e = _mm256_load_si256(pvE + j);
+			vH = _mm256_max_epi16(vH, e);
+			vH = _mm256_max_epi16(vH, vF);
+			vMaxColumn = _mm256_max_epi16(vMaxColumn, vH);
+
+			_mm256_store_si256(pvHStore + j, vH);
+
+			vH = _mm256_subs_epu16(vH, vGapO);
+			e = _mm256_subs_epu16(e, vGapE);
+			e = _mm256_max_epi16(e, vH);
+			_mm256_store_si256(pvE + j, e);
+
+			vF = _mm256_subs_epu16(vF, vGapE);
+			vF = _mm256_max_epi16(vF, vH);
+
+			vH = _mm256_load_si256(pvHLoad + j);
+		}
+
+		/* Lazy_F loop */
+		for (k = 0; LIKELY(k < 16); ++k) {
+			vF = avx2_shift_left_2_bytes(vF);
+			for (j = 0; LIKELY(j < segLen); ++j) {
+				vH = _mm256_load_si256(pvHStore + j);
+				vH = _mm256_max_epi16(vH, vF);
+				vMaxColumn = _mm256_max_epi16(vMaxColumn, vH);
+				_mm256_store_si256(pvHStore + j, vH);
+				vH = _mm256_subs_epu16(vH, vGapO);
+				vF = _mm256_subs_epu16(vF, vGapE);
+				if (UNLIKELY(!_mm256_movemask_epi8(_mm256_cmpgt_epi16(vF, vH)))) goto end;
+			}
+		}
+
+end:
+		vMaxScore = _mm256_max_epi16(vMaxScore, vMaxColumn);
+		vTemp = _mm256_cmpeq_epi16(vMaxMark, vMaxScore);
+		cmp = _mm256_movemask_epi8(vTemp);
+		if (cmp != (int32_t)0xffffffff) {
+			uint16_t temp;
+			vMaxMark = vMaxScore;
+			max16_avx2(temp, vMaxScore);
+			vMaxScore = vMaxMark;
+
+			if (LIKELY(temp > max)) {
+				max = temp;
+				end_ref = i;
+				for (j = 0; LIKELY(j < segLen); ++j) pvHmax[j] = pvHStore[j];
+			}
+		}
+
+		/* Record the max score of current column. */
+		max16_avx2(maxColumn[i], vMaxColumn);
+		if (maxColumn[i] == terminate) break;
+	}
+
+	/* Trace the alignment ending position on read. */
+	uint16_t *t = (uint16_t*)pvHmax;
+	int32_t column_len = segLen * 16;
+	for (i = 0; LIKELY(i < column_len); ++i, ++t) {
+		int32_t temp;
+		if (*t == max) {
+			temp = i / 16 + i % 16 * segLen;
+			if (temp < end_read) end_read = temp;
+		}
+	}
+
+	free(pvHmax);
+	free(pvE);
+	free(pvHLoad);
+	free(pvHStore);
+
+	alignment_end_avx2* bests = (alignment_end_avx2*)calloc(2, sizeof(alignment_end_avx2));
+	bests[0].score = max;
+	bests[0].ref = end_ref;
+	bests[0].read = end_read;
+
+	bests[1].score = 0;
+	bests[1].ref = 0;
+	bests[1].read = 0;
+
+	edge = (end_ref - maskLen) > 0 ? (end_ref - maskLen) : 0;
+	for (i = 0; i < edge; i++) {
+		if (maxColumn[i] > bests[1].score) {
+			bests[1].score = maxColumn[i];
+			bests[1].ref = i;
+		}
+	}
+	edge = (end_ref + maskLen) > refLen ? refLen : (end_ref + maskLen);
+	for (i = edge; i < refLen; i++) {
+		if (maxColumn[i] > bests[1].score) {
+			bests[1].score = maxColumn[i];
+			bests[1].ref = i;
+		}
+	}
+
+	free(maxColumn);
+	free(end_read_column);
+	return bests;
+
+	#undef max16_avx2
+}
diff --git a/cpp/ext/ssw/ssw_avx2.h b/cpp/ext/ssw/ssw_avx2.h
new file mode 100644
index 00000000..883308b8
--- /dev/null
+++ b/cpp/ext/ssw/ssw_avx2.h
@@ -0,0 +1,31 @@
+#ifndef SSW_AVX2_H
+#define SSW_AVX2_H
+
+#include <stdint.h>
+#include <immintrin.h>
+
+typedef struct {
+	uint16_t score;
+	int32_t ref;
+	int32_t read;
+} alignment_end_avx2;
+
+__m256i* qP_avx2_byte(const int8_t* read_num, const int8_t* mat,
+                       int32_t readLen, int32_t n, uint8_t bias);
+
+__m256i* qP_avx2_word(const int8_t* read_num, const int8_t* mat,
+                       int32_t readLen, int32_t n);
+
+alignment_end_avx2* sw_avx2_byte(const int8_t* ref, int8_t ref_dir,
+                                  int32_t refLen, int32_t readLen,
+                                  uint8_t weight_gapO, uint8_t weight_gapE,
+                                  const __m256i* vProfile, uint8_t terminate,
+                                  uint8_t bias, int32_t maskLen);
+
+alignment_end_avx2* sw_avx2_word(const int8_t* ref, int8_t ref_dir,
+                                  int32_t refLen, int32_t readLen,
+                                  uint8_t weight_gapO, uint8_t weight_gapE,
+                                  const __m256i* vProfile, uint16_t terminate,
+                                  int32_t maskLen);
+
+#endif
diff --git a/pixi.lock b/pixi.lock
new file mode 100644
index 00000000..8d578ee1
--- /dev/null
+++ b/pixi.lock
@@ -0,0 +1,384 @@
+version: 6
+environments:
+  default:
+    channels:
+    - url: https://conda.anaconda.org/conda-forge/
+    - url: https://conda.anaconda.org/bioconda/
+    options:
+      pypi-prerelease-mode: if-necessary-or-explicit
+    packages:
+      linux-64:
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda
+      - conda: https://conda.anaconda.org/bioconda/linux-64/bwa-mem2-2.3-he70b90d_0.tar.bz2
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.6-hb03c661_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.2.25-hbd8a1cb_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/cmake-4.2.3-hc85cc9f_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/isa-l-2.31.1-hb9d3cd8_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.3-hb9d3cd8_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.22.2-ha1258a1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.18.0-hcf29cc6_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20250104-pl5321h7949ede_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-hd590300_2.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.7.4-hecca717_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-15.2.0-he0feb66_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-15.2.0-h69a702a_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libgomp-15.2.0-he0feb66_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.2-hb03c661_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.67.0-had1ee68_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.11.1-hcf80075_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-15.2.0-h934c35e_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-15.2.0-hdf11a46_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libuv-1.51.0-hb03c661_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.1-hb9d3cd8_2.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/nasm-2.16.03-h4bc722e_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.5-h2d0b736_3.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.6.1-h35e630c_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rhash-1.4.6-hb9d3cd8_1.conda
+      - conda: https://conda.anaconda.org/bioconda/linux-64/strobealign-0.17.0-h5ca1c30_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/zlib-1.3.1-hb9d3cd8_2.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.7-hb78ec9c_6.conda
+packages:
+- conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda
+  build_number: 20
+  sha256: 1dd3fffd892081df9726d7eb7e0dea6198962ba775bd88842135a4ddb4deb3c9
+  md5: a9f577daf3de00bca7c3c76c0ecbd1de
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgomp >=7.5.0
+  constrains:
+  - openmp_impl <0.0a0
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 28948
+  timestamp: 1770939786096
+- conda: https://conda.anaconda.org/bioconda/linux-64/bwa-mem2-2.3-he70b90d_0.tar.bz2
+  sha256: 7ccf117aa586d663e12181e60f0cf60215fb162515e0b6e98517f5e108613d41
+  md5: c343a08cd26491362cc2f388e763698a
+  depends:
+  - libgcc
+  - libgcc-ng >=12
+  - libstdcxx
+  - libstdcxx-ng >=12
+  - libzlib >=1.3.1,<2.0a0
+  license: MIT
+  license_family: MIT
+  size: 3139099
+  timestamp: 1751493038078
+- conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda
+  sha256: 0b75d45f0bba3e95dc693336fa51f40ea28c980131fec438afb7ce6118ed05f6
+  md5: d2ffd7602c02f2b316fd921d39876885
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=14
+  license: bzip2-1.0.6
+  license_family: BSD
+  size: 260182
+  timestamp: 1771350215188
+- conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.6-hb03c661_0.conda
+  sha256: cc9accf72fa028d31c2a038460787751127317dcfa991f8d1f1babf216bb454e
+  md5: 920bb03579f15389b9e512095ad995b7
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=14
+  license: MIT
+  license_family: MIT
+  size: 207882
+  timestamp: 1765214722852
+- conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.2.25-hbd8a1cb_0.conda
+  sha256: 67cc7101b36421c5913a1687ef1b99f85b5d6868da3abbf6ec1a4181e79782fc
+  md5: 4492fd26db29495f0ba23f146cd5638d
+  depends:
+  - __unix
+  license: ISC
+  size: 147413
+  timestamp: 1772006283803
+- conda: https://conda.anaconda.org/conda-forge/linux-64/cmake-4.2.3-hc85cc9f_1.conda
+  sha256: 5ece78754577b8d9030ec1f09ce1cd481125f27d8d6fcdcfe2c1017661830c61
+  md5: 51d37989c1758b5edfe98518088bf700
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - bzip2 >=1.0.8,<2.0a0
+  - libcurl >=8.18.0,<9.0a0
+  - libexpat >=2.7.4,<3.0a0
+  - libgcc >=14
+  - liblzma >=5.8.2,<6.0a0
+  - libstdcxx >=14
+  - libuv >=1.51.0,<2.0a0
+  - libzlib >=1.3.1,<2.0a0
+  - ncurses >=6.5,<7.0a0
+  - rhash >=1.4.6,<2.0a0
+  - zstd >=1.5.7,<1.6.0a0
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 22330508
+  timestamp: 1771383666798
+- conda: https://conda.anaconda.org/conda-forge/linux-64/isa-l-2.31.1-hb9d3cd8_1.conda
+  sha256: 75b15f01a6b286630c4a98be0d05e286275a5ef3868e23e6d9644e51b73650e1
+  md5: 00f364ec0a7e975ec9d2fc720b19c129
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 157291
+  timestamp: 1736497194571
+- conda: https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.3-hb9d3cd8_0.conda
+  sha256: 0960d06048a7185d3542d850986d807c6e37ca2e644342dd0c72feefcf26c2a4
+  md5: b38117a3c920364aff79f870c984b4a3
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  license: LGPL-2.1-or-later
+  size: 134088
+  timestamp: 1754905959823
+- conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.22.2-ha1258a1_0.conda
+  sha256: 3e307628ca3527448dd1cb14ad7bb9d04d1d28c7d4c5f97ba196ae984571dd25
+  md5: fb53fb07ce46a575c5d004bbc96032c2
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - keyutils >=1.6.3,<2.0a0
+  - libedit >=3.1.20250104,<3.2.0a0
+  - libedit >=3.1.20250104,<4.0a0
+  - libgcc >=14
+  - libstdcxx >=14
+  - openssl >=3.5.5,<4.0a0
+  license: MIT
+  license_family: MIT
+  size: 1386730
+  timestamp: 1769769569681
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.18.0-hcf29cc6_1.conda
+  sha256: c84e8dccb65ad5149c0121e4b54bdc47fa39303fd5f4979b8c44bb51b39a369b
+  md5: 1707cdd636af2ff697b53186572c9f77
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - krb5 >=1.22.2,<1.23.0a0
+  - libgcc >=14
+  - libnghttp2 >=1.67.0,<2.0a0
+  - libssh2 >=1.11.1,<2.0a0
+  - libzlib >=1.3.1,<2.0a0
+  - openssl >=3.5.5,<4.0a0
+  - zstd >=1.5.7,<1.6.0a0
+  license: curl
+  license_family: MIT
+  size: 463621
+  timestamp: 1770892808818
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20250104-pl5321h7949ede_0.conda
+  sha256: d789471216e7aba3c184cd054ed61ce3f6dac6f87a50ec69291b9297f8c18724
+  md5: c277e0a4d549b03ac1e9d6cbbe3d017b
+  depends:
+  - ncurses
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  - ncurses >=6.5,<7.0a0
+  license: BSD-2-Clause
+  license_family: BSD
+  size: 134676
+  timestamp: 1738479519902
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-hd590300_2.conda
+  sha256: 1cd6048169fa0395af74ed5d8f1716e22c19a81a8a36f934c110ca3ad4dd27b4
+  md5: 172bf1cd1ff8629f2b1179945ed45055
+  depends:
+  - libgcc-ng >=12
+  license: BSD-2-Clause
+  license_family: BSD
+  size: 112766
+  timestamp: 1702146165126
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.7.4-hecca717_0.conda
+  sha256: d78f1d3bea8c031d2f032b760f36676d87929b18146351c4464c66b0869df3f5
+  md5: e7f7ce06ec24cfcfb9e36d28cf82ba57
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=14
+  constrains:
+  - expat 2.7.4.*
+  license: MIT
+  license_family: MIT
+  size: 76798
+  timestamp: 1771259418166
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-15.2.0-he0feb66_18.conda
+  sha256: faf7d2017b4d718951e3a59d081eb09759152f93038479b768e3d612688f83f5
+  md5: 0aa00f03f9e39fb9876085dee11a85d4
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - _openmp_mutex >=4.5
+  constrains:
+  - libgcc-ng ==15.2.0=*_18
+  - libgomp 15.2.0 he0feb66_18
+  license: GPL-3.0-only WITH GCC-exception-3.1
+  license_family: GPL
+  size: 1041788
+  timestamp: 1771378212382
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-15.2.0-h69a702a_18.conda
+  sha256: e318a711400f536c81123e753d4c797a821021fb38970cebfb3f454126016893
+  md5: d5e96b1ed75ca01906b3d2469b4ce493
+  depends:
+  - libgcc 15.2.0 he0feb66_18
+  license: GPL-3.0-only WITH GCC-exception-3.1
+  license_family: GPL
+  size: 27526
+  timestamp: 1771378224552
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libgomp-15.2.0-he0feb66_18.conda
+  sha256: 21337ab58e5e0649d869ab168d4e609b033509de22521de1bfed0c031bfc5110
+  md5: 239c5e9546c38a1e884d69effcf4c882
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  license: GPL-3.0-only WITH GCC-exception-3.1
+  license_family: GPL
+  size: 603262
+  timestamp: 1771378117851
+- conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.2-hb03c661_0.conda
+  sha256: 755c55ebab181d678c12e49cced893598f2bab22d582fbbf4d8b83c18be207eb
+  md5: c7c83eecbb72d88b940c249af56c8b17
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=14
+  constrains:
+  - xz 5.8.2.*
+  license: 0BSD
+  size: 113207
+  timestamp: 1768752626120
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.67.0-had1ee68_0.conda
+  sha256: a4a7dab8db4dc81c736e9a9b42bdfd97b087816e029e221380511960ac46c690
+  md5: b499ce4b026493a13774bcf0f4c33849
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - c-ares >=1.34.5,<2.0a0
+  - libev >=4.33,<4.34.0a0
+  - libev >=4.33,<5.0a0
+  - libgcc >=14
+  - libstdcxx >=14
+  - libzlib >=1.3.1,<2.0a0
+  - openssl >=3.5.2,<4.0a0
+  license: MIT
+  license_family: MIT
+  size: 666600
+  timestamp: 1756834976695
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.11.1-hcf80075_0.conda
+  sha256: fa39bfd69228a13e553bd24601332b7cfeb30ca11a3ca50bb028108fe90a7661
+  md5: eecce068c7e4eddeb169591baac20ac4
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  - libzlib >=1.3.1,<2.0a0
+  - openssl >=3.5.0,<4.0a0
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 304790
+  timestamp: 1745608545575
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-15.2.0-h934c35e_18.conda
+  sha256: 78668020064fdaa27e9ab65cd2997e2c837b564ab26ce3bf0e58a2ce1a525c6e
+  md5: 1b08cd684f34175e4514474793d44bcb
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc 15.2.0 he0feb66_18
+  constrains:
+  - libstdcxx-ng ==15.2.0=*_18
+  license: GPL-3.0-only WITH GCC-exception-3.1
+  license_family: GPL
+  size: 5852330
+  timestamp: 1771378262446
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-15.2.0-hdf11a46_18.conda
+  sha256: 3c902ffd673cb3c6ddde624cdb80f870b6c835f8bf28384b0016e7d444dd0145
+  md5: 6235adb93d064ecdf3d44faee6f468de
+  depends:
+  - libstdcxx 15.2.0 h934c35e_18
+  license: GPL-3.0-only WITH GCC-exception-3.1
+  license_family: GPL
+  size: 27575
+  timestamp: 1771378314494
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libuv-1.51.0-hb03c661_1.conda
+  sha256: c180f4124a889ac343fc59d15558e93667d894a966ec6fdb61da1604481be26b
+  md5: 0f03292cc56bf91a077a134ea8747118
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=14
+  license: MIT
+  license_family: MIT
+  size: 895108
+  timestamp: 1753948278280
+- conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.1-hb9d3cd8_2.conda
+  sha256: d4bfe88d7cb447768e31650f06257995601f89076080e76df55e3112d4e47dc4
+  md5: edb0dca6bc32e4f4789199455a1dbeb8
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  constrains:
+  - zlib 1.3.1 *_2
+  license: Zlib
+  license_family: Other
+  size: 60963
+  timestamp: 1727963148474
+- conda: https://conda.anaconda.org/conda-forge/linux-64/nasm-2.16.03-h4bc722e_1.conda
+  sha256: d01bfa655ad08d33dc5830a5166c7b664143df24fab59d41df15f076c58000b6
+  md5: 35f8ab79609d5bc56d6d040f12dacf3a
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc-ng >=12
+  license: BSD-2-Clause
+  license_family: BSD
+  size: 1221519
+  timestamp: 1721652638250
+- conda: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.5-h2d0b736_3.conda
+  sha256: 3fde293232fa3fca98635e1167de6b7c7fda83caf24b9d6c91ec9eefb4f4d586
+  md5: 47e340acb35de30501a76c7c799c41d7
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  license: X11 AND BSD-3-Clause
+  size: 891641
+  timestamp: 1738195959188
+- conda: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.6.1-h35e630c_1.conda
+  sha256: 44c877f8af015332a5d12f5ff0fb20ca32f896526a7d0cdb30c769df1144fb5c
+  md5: f61eb8cd60ff9057122a3d338b99c00f
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - ca-certificates
+  - libgcc >=14
+  license: Apache-2.0
+  license_family: Apache
+  size: 3164551
+  timestamp: 1769555830639
+- conda: https://conda.anaconda.org/conda-forge/linux-64/rhash-1.4.6-hb9d3cd8_1.conda
+  sha256: d5c73079c1dd2c2a313c3bfd81c73dbd066b7eb08d213778c8bff520091ae894
+  md5: c1c9b02933fdb2cfb791d936c20e887e
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  license: MIT
+  license_family: MIT
+  size: 193775
+  timestamp: 1748644872902
+- conda: https://conda.anaconda.org/bioconda/linux-64/strobealign-0.17.0-h5ca1c30_0.conda
+  sha256: ca7a962fc160859b9327ed7ab1b8c5b8ae878bc1c72e529948ae74dc2be60109
+  md5: 4d25928800b46b289f1d47b226938a94
+  depends:
+  - isa-l >=2.31.1,<3.0a0
+  - libgcc >=13
+  - libstdcxx >=13
+  - libzlib >=1.3.1,<2.0a0
+  - zlib
+  license: MIT
+  size: 280227
+  timestamp: 1766053299578
+- conda: https://conda.anaconda.org/conda-forge/linux-64/zlib-1.3.1-hb9d3cd8_2.conda
+  sha256: 5d7c0e5f0005f74112a34a7425179f4eb6e73c92f5d109e6af4ddeca407c92ab
+  md5: c9f075ab2f33b3bbee9e62d4ad0a6cd8
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libgcc >=13
+  - libzlib 1.3.1 hb9d3cd8_2
+  license: Zlib
+  license_family: Other
+  size: 92286
+  timestamp: 1727963153079
+- conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.7-hb78ec9c_6.conda
+  sha256: 68f0206ca6e98fea941e5717cec780ed2873ffabc0e1ed34428c061e2c6268c7
+  md5: 4a13eeac0b5c8e5b8ab496e6c4ddd829
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - libzlib >=1.3.1,<2.0a0
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 601375
+  timestamp: 1764777111296
diff --git a/pixi.toml b/pixi.toml
new file mode 100644
index 00000000..76bb61ce
--- /dev/null
+++ b/pixi.toml
@@ -0,0 +1,15 @@
+[workspace]
+authors = ["Pär Larsson <par.larsson.medbio@gmail.com>"]
+channels = ["conda-forge", "bioconda"]
+name = "strobealign"
+platforms = ["linux-64"]
+version = "0.1.0"
+
+[tasks]
+
+[dependencies]
+strobealign = ">=0.17.0,<0.18"
+bwa-mem2 = ">=2.3,<3"
+nasm = ">=2.16.3,<3"
+isa-l = ">=2.31.1,<3"
+cmake = ">=4.2.3,<5"

From 744a861ca210191c278d0c878531aebda1458431 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?P=C3=A4r=20Larsson?= <par.larsson.medbio@gmail.com>
Date: Tue, 3 Mar 2026 20:59:48 +0100
Subject: [PATCH 3/5] Prefetching in find_all_hits()

---
 cpp/hits.cpp  |  6 ++++++
 cpp/index.hpp | 14 ++++++++++++++
 2 files changed, 20 insertions(+)

diff --git a/cpp/hits.cpp b/cpp/hits.cpp
index bf83ea44..662da45b 100644
--- a/cpp/hits.cpp
+++ b/cpp/hits.cpp
@@ -66,6 +66,9 @@ std::tuple<std::vector<Hit>, HitsDetails, bool> find_all_hits(
     std::vector<Hit> hits;
     if (mcs_strategy != McsStrategy::FirstStrobe) {
         for (size_t i = 0; i < query_randstrobes.size(); ++i) {
+            if (i + 1 < query_randstrobes.size()) {
+                index.prefetch(query_randstrobes[i + 1].hash);
+            }
             const auto &q = query_randstrobes[i];
             size_t position = index.find_full(q.hash);
             if (position != index.end()) {
@@ -102,6 +105,9 @@ std::tuple<std::vector<Hit>, HitsDetails, bool> find_all_hits(
     ) {
         // Only partial lookups
         for (size_t i = 0; i < query_randstrobes.size(); ++i) {
+            if (i + 1 < query_randstrobes.size()) {
+                index.prefetch(query_randstrobes[i + 1].hash);
+            }
             const auto &q = query_randstrobes[i];
             size_t partial_pos = index.find_partial(q.hash);
             if (partial_pos != index.end()) {
diff --git a/cpp/index.hpp b/cpp/index.hpp
index 8ad33907..92059313 100644
--- a/cpp/index.hpp
+++ b/cpp/index.hpp
@@ -234,6 +234,20 @@ struct StrobemerIndex {
         return (pos - randstrobes.begin() - 1) - position + 1;
     }
 
+    /*
+     * Prefetch the bucket and first randstrobe entries for a given hash key.
+     * Call this for the next query while processing the current one to hide
+     * memory latency.
+     */
+    void prefetch(randstrobe_hash_t key) const {
+        const unsigned int top_N = key >> (64 - bits);
+        __builtin_prefetch(&randstrobe_start_indices[top_N], 0, 0);
+        bucket_index_t pos = randstrobe_start_indices[top_N];
+        if (pos < randstrobes.size()) {
+            __builtin_prefetch(&randstrobes[pos], 0, 0);
+        }
+    }
+
     size_t end() const {
         return -1;
     }

From 2808a5f89b552c4b81167e1fdd2054ec799891b7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?P=C3=A4r=20Larsson?= <par.larsson.medbio@gmail.com>
Date: Tue, 3 Mar 2026 21:06:13 +0100
Subject: [PATCH 4/5] hardware popcnt instructions

---
 CMakeLists.txt      | 2 +-
 cpp/randstrobes.cpp | 7 ++-----
 2 files changed, 3 insertions(+), 6 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6f0d62f4..ee7a4b4e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -29,7 +29,7 @@ set(ISAL "system" CACHE STRING "Where to find ISA-L: system (uses pkg-config), d
 set_property(CACHE ISAL PROPERTY STRINGS "system" "download")
 
 message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
-add_compile_options(-Wall -Wextra -Werror=maybe-uninitialized)
+add_compile_options(-Wall -Wextra -Werror=maybe-uninitialized -mpopcnt)
 
 add_subdirectory(cpp/ext/zstr)
 
diff --git a/cpp/randstrobes.cpp b/cpp/randstrobes.cpp
index b171916e..b46ffddd 100644
--- a/cpp/randstrobes.cpp
+++ b/cpp/randstrobes.cpp
@@ -1,5 +1,4 @@
 #include <deque>
-#include <bitset>
 #include <algorithm>
 #include <cassert>
 #include <array>
@@ -176,8 +175,7 @@ Randstrobe RandstrobeIterator::get(unsigned int strobe1_index) const {
         assert(i < syncmers.size());
 
         // Method 3' skew sample more for prob exact matching
-        std::bitset<64> b = (strobe1.hash ^ syncmers[i].hash) & parameters.q;
-        uint64_t res = b.count();
+        uint64_t res = __builtin_popcountll((strobe1.hash ^ syncmers[i].hash) & parameters.q);
 
         if (res < min_val) {
             min_val = res;
@@ -207,8 +205,7 @@ Randstrobe RandstrobeGenerator::next() {
     for (auto i = parameters.w_min; i < syncmers.size() && syncmers[i].position <= max_position; i++) {
         assert(i <= parameters.w_max);
         // Method 3' skew sample more for prob exact matching
-        std::bitset<64> b = (strobe1.hash ^ syncmers[i].hash) & parameters.q;
-        uint64_t res = b.count();
+        uint64_t res = __builtin_popcountll((strobe1.hash ^ syncmers[i].hash) & parameters.q);
 
         if (res < min_val) {
             min_val = res;

From 7c222d9cde1d75e29be349c2fb533db224f03faf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?P=C3=A4r=20Larsson?= <par.larsson.medbio@gmail.com>
Date: Tue, 3 Mar 2026 21:52:15 +0100
Subject: [PATCH 5/5] Circular buffer

---
 cpp/aln.cpp         | 12 ++++++------
 cpp/randstrobes.cpp |  4 ++++
 cpp/randstrobes.hpp | 22 +++++++++++++++++++++-
 3 files changed, 31 insertions(+), 7 deletions(-)

diff --git a/cpp/aln.cpp b/cpp/aln.cpp
index 53220dba..ce6a764d 100644
--- a/cpp/aln.cpp
+++ b/cpp/aln.cpp
@@ -1019,12 +1019,12 @@ void shuffle_top_nams(std::vector<Nam>& nams, std::minstd_rand& random_engine) {
  * in common with the reference sequence
  */
 bool has_shared_substring(const std::string& read_seq, const std::string& ref_seq, int k) {
-    int sub_size = 2 * k / 3;
-    int step_size = k / 3;
-    std::string submer;
-    for (size_t i = 0; i + sub_size < read_seq.size(); i += step_size) {
-        submer = read_seq.substr(i, sub_size);
-        if (ref_seq.find(submer) != std::string::npos) {
+    size_t sub_size = 2 * k / 3;
+    size_t step_size = k / 3;
+    std::string_view ref_view(ref_seq);
+    std::string_view read_view(read_seq);
+    for (size_t i = 0; i + sub_size <= read_view.size(); i += step_size) {
+        if (ref_view.find(read_view.substr(i, sub_size)) != std::string_view::npos) {
             return true;
         }
     }
diff --git a/cpp/randstrobes.cpp b/cpp/randstrobes.cpp
index b46ffddd..8661d768 100644
--- a/cpp/randstrobes.cpp
+++ b/cpp/randstrobes.cpp
@@ -232,6 +232,10 @@ std::array<std::vector<QueryRandstrobe>, 2> randstrobes_query(const std::string_
         return randstrobes;
     }
 
+    // Pre-allocate: roughly one randstrobe per syncmer
+    randstrobes[0].reserve(syncmers.size());
+    randstrobes[1].reserve(syncmers.size());
+
     // Generate randstrobes for the forward sequence
     RandstrobeIterator randstrobe_fwd_iter{syncmers, parameters.randstrobe};
     while (randstrobe_fwd_iter.has_next()) {
diff --git a/cpp/randstrobes.hpp b/cpp/randstrobes.hpp
index a77e15e8..f79065fe 100644
--- a/cpp/randstrobes.hpp
+++ b/cpp/randstrobes.hpp
@@ -136,6 +136,26 @@ class RandstrobeIterator {
 
 std::ostream& operator<<(std::ostream& os, const Syncmer& syncmer);
 
+/* Fixed-capacity circular buffer for s-mer hashes in SyncmerIterator.
+ * Replaces std::deque for better cache locality (all data inline, no heap).
+ * Capacity 32 covers worst case: k < 32, so k - s + 1 <= 31. */
+class SmerBuffer {
+public:
+    SmerBuffer() = default;
+    void push_back(uint64_t val) { m_data[(m_start + m_size++) & MASK] = val; }
+    uint64_t front() const { return m_data[m_start]; }
+    void pop_front() { m_start = (m_start + 1) & MASK; m_size--; }
+    uint64_t operator[](size_t i) const { return m_data[(m_start + i) & MASK]; }
+    size_t size() const { return m_size; }
+    void clear() { m_start = 0; m_size = 0; }
+private:
+    static constexpr size_t CAP = 32;  // power of 2; must be >= k - s + 1 (max 31)
+    static constexpr size_t MASK = CAP - 1;
+    uint64_t m_data[CAP] = {};
+    size_t m_start = 0;
+    size_t m_size = 0;
+};
+
 class SyncmerIterator {
 public:
     SyncmerIterator(const std::string_view seq, SyncmerParameters parameters)
@@ -151,7 +171,7 @@ class SyncmerIterator {
     const uint64_t smask = (1ULL << 2*parameters.s) - 1;
     const uint64_t kshift = (parameters.k - 1) * 2;
     const uint64_t sshift = (parameters.s - 1) * 2;
-    std::deque<uint64_t> qs;  // s-mer hashes
+    SmerBuffer qs;  // s-mer hashes (circular buffer, capacity 8)
     uint64_t qs_min_val = UINT64_MAX;
     size_t l = 0;
     uint64_t xk[2] = {0, 0};