aboutsummaryrefslogtreecommitdiff
path: root/src/maths
diff options
context:
space:
mode:
Diffstat (limited to 'src/maths')
-rw-r--r--src/maths/complex.c39
-rw-r--r--src/maths/format.c12
-rw-r--r--src/maths/matrix.c17
-rw-r--r--src/maths/simd.c324
-rw-r--r--src/maths/vector.c12
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;