feat: simd + testing + formatting

5 files changed, 358 insertions, 46 deletions

diff --git a/src/maths/complex.c b/src/maths/complex.c
index 3115fdd..44e255f 100644
--- a/src/maths/complex.c
+++ b/src/maths/complex.c
@@ -4,8 +4,7 @@
 
 struct PsiComplex psi_new_complex(double real, double imaginary)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		real,
 		imaginary,
 	};
@@ -13,8 +12,7 @@ struct PsiComplex psi_new_complex(double real, double imaginary)
 
 struct PsiComplex psi_new_complex_from_real(double real)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		real,
 		0.0,
 	};
@@ -22,8 +20,7 @@ struct PsiComplex psi_new_complex_from_real(double real)
 
 struct PsiComplex psi_conjugate_complex(struct PsiComplex z)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		z.real,
 		-z.imaginary,
 	};
@@ -31,8 +28,7 @@ struct PsiComplex psi_conjugate_complex(struct PsiComplex z)
 
 struct PsiComplex psi_neg_complex(struct PsiComplex z)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		-z.real,
 		-z.imaginary,
 	};
@@ -59,8 +55,7 @@ struct PsiComplex psi_sqrt_complex(struct PsiComplex z)
 	double half_theta = psi_phase_complex(z) / 2.0;
 	double sqrt_r = sqrt(r);
 
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		sqrt_r * cos(half_theta),
 		sqrt_r * sin(half_theta),
 	};
@@ -68,8 +63,7 @@ struct PsiComplex psi_sqrt_complex(struct PsiComplex z)
 
 struct PsiComplex psi_add_complex(struct PsiComplex a, struct PsiComplex b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real + b.real,
 		a.imaginary + b.imaginary,
 	};
@@ -77,8 +71,7 @@ struct PsiComplex psi_add_complex(struct PsiComplex a, struct PsiComplex b)
 
 struct PsiComplex psi_sub_complex(struct PsiComplex a, struct PsiComplex b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real - b.real,
 		a.imaginary - b.imaginary,
 	};
@@ -87,8 +80,7 @@ struct PsiComplex psi_sub_complex(struct PsiComplex a, struct PsiComplex b)
 struct PsiComplex psi_mul_complex(struct PsiComplex a, struct PsiComplex b)
 {
 	// (a + bi)(c + di) = (ac - bd) + (ad + bc)i
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real * b.real - a.imaginary * b.imaginary,
 		a.real * b.imaginary + a.imaginary * b.real,
 	};
@@ -98,8 +90,7 @@ struct PsiComplex psi_div_complex(struct PsiComplex a, struct PsiComplex b)
 {
 	// (a + bi) / (c + di) = ((ac + bd) + (bc - ad)i) / (c² + d²)
 	double denom = b.real * b.real + b.imaginary * b.imaginary;
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		(a.real * b.real + a.imaginary * b.imaginary) / denom,
 		(a.imaginary * b.real - a.real * b.imaginary) / denom,
 	};
@@ -107,8 +98,7 @@ struct PsiComplex psi_div_complex(struct PsiComplex a, struct PsiComplex b)
 
 struct PsiComplex psi_add_complex_real(struct PsiComplex a, double b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real + b,
 		a.imaginary,
 	};
@@ -116,8 +106,7 @@ struct PsiComplex psi_add_complex_real(struct PsiComplex a, double b)
 
 struct PsiComplex psi_sub_complex_real(struct PsiComplex a, double b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real - b,
 		a.imaginary,
 	};
@@ -125,8 +114,7 @@ struct PsiComplex psi_sub_complex_real(struct PsiComplex a, double b)
 
 struct PsiComplex psi_mul_complex_real(struct PsiComplex a, double b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real * b,
 		a.imaginary * b,
 	};
@@ -134,8 +122,7 @@ struct PsiComplex psi_mul_complex_real(struct PsiComplex a, double b)
 
 struct PsiComplex psi_div_complex_real(struct PsiComplex a, double b)
 {
-	return (struct PsiComplex)
-	{
+	return (struct PsiComplex){
 		a.real / b,
 		a.imaginary / b,
 	};
diff --git a/src/maths/format.c b/src/maths/format.c
index ad80ff0..f3db83d 100644
--- a/src/maths/format.c
+++ b/src/maths/format.c
@@ -15,10 +15,10 @@ static bool approx_eq(double a, double b)
 	return fabs(a - b) < EPSILON;
 }
 
-static bool real_symbolic(double v, char *out, size_t cap)
+static bool real_symbolic(double v, char* out, size_t cap)
 {
 	double abs_v = fabs(v);
-	const char *sign = v < 0.0 ? "-" : "";
+	const char* sign = v < 0.0 ? "-" : "";
 
 	if (approx_eq(abs_v, 0.0))
 	{
@@ -26,7 +26,7 @@ static bool real_symbolic(double v, char *out, size_t cap)
 		return true;
 	}
 
-	const char *sym = NULL;
+	const char* sym = NULL;
 	if (approx_eq(abs_v, 1.0))
 		sym = "1";
 	else if (approx_eq(abs_v, 0.5))
@@ -58,7 +58,7 @@ static bool real_symbolic(double v, char *out, size_t cap)
 	return true;
 }
 
-char *psi_format_amplitude(struct PsiComplex c, char *out, size_t cap)
+char* psi_format_amplitude(struct PsiComplex c, char* out, size_t cap)
 {
 	double re = c.real;
 	double im = c.imaginary;
@@ -108,7 +108,7 @@ char *psi_format_amplitude(struct PsiComplex c, char *out, size_t cap)
 	}
 	else
 	{
-		const char *sign = im > 0.0 ? "+" : "-";
+		const char* sign = im > 0.0 ? "+" : "-";
 		char sym[32];
 		if (real_symbolic(fabs(im), sym, sizeof sym))
 			snprintf(im_str, sizeof im_str, "%s%si", sign, sym);
@@ -120,7 +120,7 @@ char *psi_format_amplitude(struct PsiComplex c, char *out, size_t cap)
 	return out;
 }
 
-char *psi_format_probability(double p, char *out, size_t cap)
+char* psi_format_probability(double p, char* out, size_t cap)
 {
 	if (approx_eq(p, 0.0))
 		snprintf(out, cap, "%s", "0");
diff --git a/src/maths/matrix.c b/src/maths/matrix.c
index 44b953c..2415395 100644
--- a/src/maths/matrix.c
+++ b/src/maths/matrix.c
@@ -6,18 +6,17 @@
 
 struct PsiMatrix psi_new_matrix(size_t rows, size_t cols)
 {
-	struct PsiComplex *data = calloc(rows * cols, sizeof(struct PsiComplex));
+	struct PsiComplex* data = calloc(rows * cols, sizeof(struct PsiComplex));
 	assert(data != NULL || rows * cols == 0);
 
-	return (struct PsiMatrix)
-	{
+	return (struct PsiMatrix){
 		data,
 		rows,
 		cols,
 	};
 }
 
-struct PsiMatrix psi_new_matrix_from(const struct PsiComplex *data, size_t rows, size_t cols)
+struct PsiMatrix psi_new_matrix_from(const struct PsiComplex* data, size_t rows, size_t cols)
 {
 	struct PsiMatrix m = psi_new_matrix(rows, cols);
 	memcpy(m.data, data, rows * cols * sizeof(struct PsiComplex));
@@ -39,7 +38,7 @@ struct PsiMatrix psi_identity_matrix(size_t size)
 	return m;
 }
 
-void psi_free_matrix(struct PsiMatrix *m)
+void psi_free_matrix(struct PsiMatrix* m)
 {
 	free(m->data);
 	m->data = NULL;
@@ -53,7 +52,7 @@ struct PsiComplex psi_get_matrix(struct PsiMatrix m, size_t row, size_t col)
 	return m.data[row * m.cols + col];
 }
 
-void psi_set_matrix(struct PsiMatrix *m, size_t row, size_t col, struct PsiComplex value)
+void psi_set_matrix(struct PsiMatrix* m, size_t row, size_t col, struct PsiComplex value)
 {
 	assert(row < m->rows && col < m->cols);
 	m->data[row * m->cols + col] = value;
@@ -69,7 +68,8 @@ struct PsiMatrix psi_dot_matrix(struct PsiMatrix a, struct PsiMatrix b)
 		{
 			struct PsiComplex sum = psi_new_complex(0.0, 0.0);
 			for (size_t k = 0; k < a.cols; k++)
-				sum = psi_add_complex(sum, psi_mul_complex(a.data[i * a.cols + k], b.data[k * b.cols + j]));
+				sum = psi_add_complex(
+						sum, psi_mul_complex(a.data[i * a.cols + k], b.data[k * b.cols + j]));
 
 			result.data[i * result.cols + j] = sum;
 		}
@@ -89,7 +89,8 @@ struct PsiMatrix psi_kronecker_matrix(struct PsiMatrix a, struct PsiMatrix b)
 				{
 					size_t row = i * b.rows + k;
 					size_t col = j * b.cols + l;
-					result.data[row * result.cols + col] = psi_mul_complex(a_val, b.data[k * b.cols + l]);
+					result.data[row * result.cols + col] =
+							psi_mul_complex(a_val, b.data[k * b.cols + l]);
 				}
 		}
 
diff --git a/src/maths/simd.c b/src/maths/simd.c
new file mode 100644
index 0000000..9e83dcc
--- /dev/null
+++ b/src/maths/simd.c
@@ -0,0 +1,324 @@
+#include "maths/simd.h"
+
+#include <assert.h>
+#include <stdlib.h>
+
+#if defined(__x86_64__) || defined(__i386__)
+#include <immintrin.h>
+#elif defined(__aarch64__)
+#include <arm_neon.h>
+#endif
+
+enum PsiSimdCapability psi_detect_simd(void)
+{
+#if defined(__x86_64__) || defined(__i386__)
+	if (__builtin_cpu_supports("avx512f") && __builtin_cpu_supports("avx512dq"))
+		return PSI_SIMD_AVX512;
+	if (__builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma"))
+		return PSI_SIMD_AVX2;
+
+	return PSI_SIMD_NONE;
+#elif defined(__aarch64__)
+	return PSI_SIMD_NEON;
+#else
+	return PSI_SIMD_NONE;
+#endif
+}
+
+const char* psi_simd_name(enum PsiSimdCapability cap)
+{
+	switch (cap)
+	{
+		case PSI_SIMD_NONE: return "Scalar";
+		case PSI_SIMD_AVX2: return "AVX2+FMA";
+		case PSI_SIMD_AVX512: return "AVX-512";
+		case PSI_SIMD_NEON: return "NEON";
+	}
+
+	return "Scalar";
+}
+
+static void apply_pair(struct PsiComplex* state, size_t i, size_t j, struct PsiComplex g00,
+                       struct PsiComplex g01, struct PsiComplex g10, struct PsiComplex g11)
+{
+	struct PsiComplex s0 = state[i];
+	struct PsiComplex s1 = state[j];
+
+	state[i] = psi_new_complex(s0.real * g00.real - s0.imaginary * g00.imaginary +
+	                                   s1.real * g01.real - s1.imaginary * g01.imaginary,
+	                           s0.real * g00.imaginary + s0.imaginary * g00.real +
+	                                   s1.real * g01.imaginary + s1.imaginary * g01.real);
+	state[j] = psi_new_complex(s0.real * g10.real - s0.imaginary * g10.imaginary +
+	                                   s1.real * g11.real - s1.imaginary * g11.imaginary,
+	                           s0.real * g10.imaginary + s0.imaginary * g10.real +
+	                                   s1.real * g11.imaginary + s1.imaginary * g11.real);
+}
+
+static void apply_scalar(struct PsiComplex* state, const struct PsiComplex gate[2][2],
+                         size_t target, size_t num_qubits)
+{
+	size_t target_bit = num_qubits - 1 - target;
+	size_t step = (size_t)1 << target_bit;
+	size_t dim = (size_t)1 << num_qubits;
+
+	for (size_t i = 0; i < dim; i++)
+	{
+		if ((i >> target_bit) & 1)
+			continue;
+
+		apply_pair(state, i, i | step, gate[0][0], gate[0][1], gate[1][0], gate[1][1]);
+	}
+}
+
+static size_t (*build_pairs(size_t dim, size_t target_bit, size_t step, size_t* out_count))[2]
+{
+	size_t (*pairs)[2] = malloc((dim / 2) * sizeof(*pairs));
+	assert(pairs != NULL || dim == 0);
+
+	size_t n = 0;
+	for (size_t i = 0; i < dim; i++)
+		if (((i >> target_bit) & 1) == 0)
+		{
+			pairs[n][0] = i;
+			pairs[n][1] = i | step;
+			n++;
+		}
+
+	*out_count = n;
+	return pairs;
+}
+
+#if defined(__aarch64__)
+static void apply_neon(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target,
+                       size_t num_qubits)
+{
+	size_t target_bit = num_qubits - 1 - target;
+	size_t step = (size_t)1 << target_bit;
+	size_t dim = (size_t)1 << num_qubits;
+
+	struct PsiComplex g00 = gate[0][0];
+	struct PsiComplex g01 = gate[0][1];
+	struct PsiComplex g10 = gate[1][0];
+	struct PsiComplex g11 = gate[1][1];
+
+	size_t np;
+	size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np);
+	size_t chunks = np / 2;
+
+	float64x2_t g00_re = vdupq_n_f64(g00.real);
+	float64x2_t g00_im = vdupq_n_f64(g00.imaginary);
+	float64x2_t g01_re = vdupq_n_f64(g01.real);
+	float64x2_t g01_im = vdupq_n_f64(g01.imaginary);
+	float64x2_t g10_re = vdupq_n_f64(g10.real);
+	float64x2_t g10_im = vdupq_n_f64(g10.imaginary);
+	float64x2_t g11_re = vdupq_n_f64(g11.real);
+	float64x2_t g11_im = vdupq_n_f64(g11.imaginary);
+
+	for (size_t chunk = 0; chunk < chunks; chunk++)
+	{
+		size_t i0 = pairs[chunk * 2][0];
+		size_t j0 = pairs[chunk * 2][1];
+		size_t i1 = pairs[chunk * 2 + 1][0];
+		size_t j1 = pairs[chunk * 2 + 1][1];
+
+		double s0r[2] = { state[i0].real, state[i1].real };
+		double s0i[2] = { state[i0].imaginary, state[i1].imaginary };
+		double s1r[2] = { state[j0].real, state[j1].real };
+		double s1i[2] = { state[j0].imaginary, state[j1].imaginary };
+
+		float64x2_t s0_re = vld1q_f64(s0r);
+		float64x2_t s0_im = vld1q_f64(s0i);
+		float64x2_t s1_re = vld1q_f64(s1r);
+		float64x2_t s1_im = vld1q_f64(s1i);
+
+		float64x2_t new0_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im),
+		                                vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im));
+		float64x2_t new0_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re),
+		                                vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re));
+		float64x2_t new1_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im),
+		                                vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im));
+		float64x2_t new1_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re),
+		                                vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re));
+
+		state[i0] = psi_new_complex(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0));
+		state[j0] = psi_new_complex(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0));
+		state[i1] = psi_new_complex(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1));
+		state[j1] = psi_new_complex(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1));
+	}
+
+	for (size_t p = chunks * 2; p < np; p++)
+		apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11);
+
+	free(pairs);
+}
+#endif
+
+#if defined(__x86_64__) || defined(__i386__)
+__attribute__((target("avx2,fma"))) static void apply_avx2(struct PsiComplex* state,
+                                                           const struct PsiComplex gate[2][2],
+                                                           size_t target, size_t num_qubits)
+{
+	size_t target_bit = num_qubits - 1 - target;
+	size_t step = (size_t)1 << target_bit;
+	size_t dim = (size_t)1 << num_qubits;
+
+	struct PsiComplex g00 = gate[0][0];
+	struct PsiComplex g01 = gate[0][1];
+	struct PsiComplex g10 = gate[1][0];
+	struct PsiComplex g11 = gate[1][1];
+
+	size_t np;
+	size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np);
+	size_t chunks = np / 2;
+
+	for (size_t chunk = 0; chunk < chunks; chunk++)
+	{
+		size_t i0 = pairs[chunk * 2][0];
+		size_t j0 = pairs[chunk * 2][1];
+		size_t i1 = pairs[chunk * 2 + 1][0];
+		size_t j1 = pairs[chunk * 2 + 1][1];
+
+		__m256d s0_re =
+				_mm256_set_pd(state[j1].real, state[i1].real, state[j0].real, state[i0].real);
+		__m256d s0_im = _mm256_set_pd(state[j1].imaginary, state[i1].imaginary, state[j0].imaginary,
+		                              state[i0].imaginary);
+
+		__m256d g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real);
+		__m256d g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary);
+		__m256d g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real);
+		__m256d g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary);
+
+		__m256d prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0));
+		__m256d prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0));
+		__m256d prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1));
+		__m256d prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1));
+
+		double res0_re[4];
+		double res0_im[4];
+		double res1_re[4];
+		double res1_im[4];
+
+		_mm256_storeu_pd(res0_re, prod0_re);
+		_mm256_storeu_pd(res0_im, prod0_im);
+		_mm256_storeu_pd(res1_re, prod1_re);
+		_mm256_storeu_pd(res1_im, prod1_im);
+
+		state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]);
+		state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]);
+		state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]);
+		state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]);
+	}
+
+	for (size_t p = chunks * 2; p < np; p++)
+		apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11);
+
+	free(pairs);
+}
+
+__attribute__((target("avx512f,avx512dq"))) static void
+apply_avx512(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target,
+             size_t num_qubits)
+{
+	size_t target_bit = num_qubits - 1 - target;
+	size_t step = (size_t)1 << target_bit;
+	size_t dim = (size_t)1 << num_qubits;
+
+	struct PsiComplex g00 = gate[0][0];
+	struct PsiComplex g01 = gate[0][1];
+	struct PsiComplex g10 = gate[1][0];
+	struct PsiComplex g11 = gate[1][1];
+
+	size_t np;
+	size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np);
+	size_t chunks = np / 4;
+
+	for (size_t chunk = 0; chunk < chunks; chunk++)
+	{
+		size_t base = chunk * 4;
+		size_t i0 = pairs[base][0], j0 = pairs[base][1];
+		size_t i1 = pairs[base + 1][0], j1 = pairs[base + 1][1];
+		size_t i2 = pairs[base + 2][0], j2 = pairs[base + 2][1];
+		size_t i3 = pairs[base + 3][0], j3 = pairs[base + 3][1];
+
+		__m512d s0_re =
+				_mm512_set_pd(state[j3].real, state[i3].real, state[j2].real, state[i2].real,
+				              state[j1].real, state[i1].real, state[j0].real, state[i0].real);
+		__m512d s0_im = _mm512_set_pd(state[j3].imaginary, state[i3].imaginary, state[j2].imaginary,
+		                              state[i2].imaginary, state[j1].imaginary, state[i1].imaginary,
+		                              state[j0].imaginary, state[i0].imaginary);
+
+		__m512d g_re_0 = _mm512_set_pd(g01.real, g00.real, g01.real, g00.real, g01.real, g00.real,
+		                               g01.real, g00.real);
+		__m512d g_im_0 = _mm512_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary,
+		                               g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary);
+		__m512d g_re_1 = _mm512_set_pd(g11.real, g10.real, g11.real, g10.real, g11.real, g10.real,
+		                               g11.real, g10.real);
+		__m512d g_im_1 = _mm512_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary,
+		                               g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary);
+
+		__m512d prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0));
+		__m512d prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0));
+		__m512d prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1));
+		__m512d prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1));
+
+		double res0_re[8];
+		double res0_im[8];
+		double res1_re[8];
+		double res1_im[8];
+
+		_mm512_storeu_pd(res0_re, prod0_re);
+		_mm512_storeu_pd(res0_im, prod0_im);
+		_mm512_storeu_pd(res1_re, prod1_re);
+		_mm512_storeu_pd(res1_im, prod1_im);
+
+		state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]);
+		state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]);
+		state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]);
+		state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]);
+		state[i2] = psi_new_complex(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]);
+		state[j2] = psi_new_complex(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]);
+		state[i3] = psi_new_complex(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]);
+		state[j3] = psi_new_complex(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]);
+	}
+
+	for (size_t p = chunks * 4; p < np; p++)
+		apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11);
+
+	free(pairs);
+}
+#endif
+
+void psi_apply_single_qubit_gate_simd(struct PsiComplex* state, const struct PsiComplex gate[2][2],
+                                      size_t target, size_t num_qubits)
+{
+	enum PsiSimdCapability cap = psi_detect_simd();
+
+#if defined(__x86_64__) || defined(__i386__)
+	if (cap == PSI_SIMD_AVX512)
+	{
+		apply_avx512(state, gate, target, num_qubits);
+		return;
+	}
+	if (cap == PSI_SIMD_AVX2)
+	{
+		apply_avx2(state, gate, target, num_qubits);
+		return;
+	}
+#elif defined(__aarch64__)
+	if (cap == PSI_SIMD_NEON)
+	{
+		apply_neon(state, gate, target, num_qubits);
+		return;
+	}
+#endif
+
+	(void)cap;
+	apply_scalar(state, gate, target, num_qubits);
+}
+
+void psi_apply_single_qubit_gate_simd_parallel(struct PsiComplex* state,
+                                               const struct PsiComplex gate[2][2], size_t target,
+                                               size_t num_qubits)
+{
+	apply_scalar(state, gate, target, num_qubits);
+}
diff --git a/src/maths/vector.c b/src/maths/vector.c
index b9713ae..1c5c4f9 100644
--- a/src/maths/vector.c
+++ b/src/maths/vector.c
@@ -14,18 +14,18 @@ static enum PsiVectorKind flip_kind(enum PsiVectorKind kind)
 
 struct PsiVector psi_new_vector(size_t size, enum PsiVectorKind kind)
 {
-	struct PsiComplex *data = calloc(size, sizeof(struct PsiComplex));
+	struct PsiComplex* data = calloc(size, sizeof(struct PsiComplex));
 	assert(data != NULL || size == 0);
 
-	return (struct PsiVector)
-	{
+	return (struct PsiVector){
 		data,
 		size,
 		kind,
 	};
 }
 
-struct PsiVector psi_new_vector_from(const struct PsiComplex *data, size_t size, enum PsiVectorKind kind)
+struct PsiVector psi_new_vector_from(const struct PsiComplex* data, size_t size,
+                                     enum PsiVectorKind kind)
 {
 	struct PsiVector v = psi_new_vector(size, kind);
 	memcpy(v.data, data, size * sizeof(struct PsiComplex));
@@ -38,7 +38,7 @@ struct PsiVector psi_clone_vector(struct PsiVector v)
 	return psi_new_vector_from(v.data, v.size, v.kind);
 }
 
-void psi_free_vector(struct PsiVector *v)
+void psi_free_vector(struct PsiVector* v)
 {
 	free(v->data);
 	v->data = NULL;
@@ -51,7 +51,7 @@ struct PsiComplex psi_get_vector(struct PsiVector v, size_t index)
 	return v.data[index];
 }
 
-void psi_set_vector(struct PsiVector *v, size_t index, struct PsiComplex value)
+void psi_set_vector(struct PsiVector* v, size_t index, struct PsiComplex value)
 {
 	assert(index < v->size);
 	v->data[index] = value;