feat: simd + testing + formatting

1 files changed, 115 insertions, 86 deletions

diff --git a/src/core/noise.c b/src/core/noise.c
index 9aaf1d9..5f141c8 100644
--- a/src/core/noise.c
+++ b/src/core/noise.c
@@ -5,30 +5,30 @@
 #include <stdlib.h>
 #include <string.h>
 
-struct PsiKrausOperator psi_new_kraus_operator(const char *name, struct PsiMatrix matrix)
+struct PsiKrausOperator psi_new_kraus_operator(const char* name, struct PsiMatrix matrix)
 {
-	return (struct PsiKrausOperator)
-	{
+	return (struct PsiKrausOperator){
 		matrix,
 		name,
 	};
 }
 
-void psi_free_kraus_operator(struct PsiKrausOperator *op)
+void psi_free_kraus_operator(struct PsiKrausOperator* op)
 {
 	psi_free_matrix(&op->matrix);
 }
 
-struct PsiNoiseChannel psi_new_noise_channel(const char *name, const struct PsiKrausOperator *operators, size_t count, size_t num_qubits)
+struct PsiNoiseChannel psi_new_noise_channel(const char* name,
+                                             const struct PsiKrausOperator* operators, size_t count,
+                                             size_t num_qubits)
 {
-	struct PsiKrausOperator *owned = malloc(count * sizeof(struct PsiKrausOperator));
+	struct PsiKrausOperator* owned = malloc(count * sizeof(struct PsiKrausOperator));
 	assert(owned != NULL || count == 0);
 
 	if (count > 0)
 		memcpy(owned, operators, count * sizeof(struct PsiKrausOperator));
 
-	return (struct PsiNoiseChannel)
-	{
+	return (struct PsiNoiseChannel){
 		name,
 		owned,
 		count,
@@ -36,7 +36,7 @@ struct PsiNoiseChannel psi_new_noise_channel(const char *name, const struct PsiK
 	};
 }
 
-void psi_free_noise_channel(struct PsiNoiseChannel *channel)
+void psi_free_noise_channel(struct PsiNoiseChannel* channel)
 {
 	for (size_t i = 0; i < channel->operator_count; i++)
 		psi_free_matrix(&channel->operators[i].matrix);
@@ -52,18 +52,22 @@ struct PsiNoiseChannel psi_depolarising_channel(double p)
 	double sqrt_p3 = sqrt(p / 3.0);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0", psi_matrix(2, 2,
-			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
-		psi_new_kraus_operator("K1(X)", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p3, 0.0),
-			psi_new_complex(sqrt_p3, 0.0), psi_new_complex(0.0, 0.0))),
-		psi_new_kraus_operator("K2(Y)", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, -sqrt_p3),
-			psi_new_complex(0.0, sqrt_p3), psi_new_complex(0.0, 0.0))),
-		psi_new_kraus_operator("K3(Z)", psi_matrix(2, 2,
-			psi_new_complex(sqrt_p3, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(-sqrt_p3, 0.0))),
+		psi_new_kraus_operator("K0",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_1_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator(
+				"K1(X)",
+				psi_matrix(2, 2, psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p3, 0.0),
+				           psi_new_complex(sqrt_p3, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator(
+				"K2(Y)",
+				psi_matrix(2, 2, psi_new_complex(0.0, 0.0), psi_new_complex(0.0, -sqrt_p3),
+				           psi_new_complex(0.0, sqrt_p3), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K3(Z)",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_p3, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(-sqrt_p3, 0.0))),
 	};
 
 	return psi_new_noise_channel("Depolarising", ops, 4, 1);
@@ -75,12 +79,14 @@ struct PsiNoiseChannel psi_amplitude_damping_channel(double gamma)
 	double sqrt_1_gamma = sqrt(1.0 - gamma);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0", psi_matrix(2, 2,
-			psi_new_complex(1.0, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_gamma, 0.0))),
-		psi_new_kraus_operator("K1", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_gamma, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K0",
+		                       psi_matrix(2, 2, psi_new_complex(1.0, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1",
+		                       psi_matrix(2, 2, psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_gamma, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
 	};
 
 	return psi_new_noise_channel("AmplitudeDamping", ops, 2, 1);
@@ -92,12 +98,14 @@ struct PsiNoiseChannel psi_phase_damping_channel(double gamma)
 	double sqrt_1_gamma = sqrt(1.0 - gamma);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0", psi_matrix(2, 2,
-			psi_new_complex(1.0, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_gamma, 0.0))),
-		psi_new_kraus_operator("K1", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_gamma, 0.0))),
+		psi_new_kraus_operator("K0",
+		                       psi_matrix(2, 2, psi_new_complex(1.0, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1",
+		                       psi_matrix(2, 2, psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_gamma, 0.0))),
 	};
 
 	return psi_new_noise_channel("PhaseDamping", ops, 2, 1);
@@ -109,12 +117,14 @@ struct PsiNoiseChannel psi_bit_flip_channel(double p)
 	double sqrt_p = sqrt(p);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
-			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
-		psi_new_kraus_operator("K1(X)", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p, 0.0),
-			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K0(I)",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_1_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(X)",
+		                       psi_matrix(2, 2, psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_p, 0.0),
+		                                  psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0))),
 	};
 
 	return psi_new_noise_channel("BitFlip", ops, 2, 1);
@@ -126,12 +136,14 @@ struct PsiNoiseChannel psi_phase_flip_channel(double p)
 	double sqrt_p = sqrt(p);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
-			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
-		psi_new_kraus_operator("K1(Z)", psi_matrix(2, 2,
-			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(-sqrt_p, 0.0))),
+		psi_new_kraus_operator("K0(I)",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_1_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(Z)",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(-sqrt_p, 0.0))),
 	};
 
 	return psi_new_noise_channel("PhaseFlip", ops, 2, 1);
@@ -143,12 +155,14 @@ struct PsiNoiseChannel psi_bit_phase_flip_channel(double p)
 	double sqrt_p = sqrt(p);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
-			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
-		psi_new_kraus_operator("K1(Y)", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, -sqrt_p),
-			psi_new_complex(0.0, sqrt_p), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K0(I)",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_1_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(Y)",
+		                       psi_matrix(2, 2, psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(0.0, -sqrt_p),
+		                                  psi_new_complex(0.0, sqrt_p), psi_new_complex(0.0, 0.0))),
 	};
 
 	return psi_new_noise_channel("BitPhaseFlip", ops, 2, 1);
@@ -162,18 +176,22 @@ struct PsiNoiseChannel psi_generalised_amplitude_damping_channel(double p, doubl
 	double sqrt_1_gamma = sqrt(1.0 - gamma);
 
 	struct PsiKrausOperator ops[] = {
-		psi_new_kraus_operator("K0", psi_matrix(2, 2,
-			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p * sqrt_1_gamma, 0.0))),
-		psi_new_kraus_operator("K1", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p * sqrt_gamma, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
-		psi_new_kraus_operator("K2", psi_matrix(2, 2,
-			psi_new_complex(sqrt_1_p * sqrt_1_gamma, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
-		psi_new_kraus_operator("K3", psi_matrix(2, 2,
-			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
-			psi_new_complex(sqrt_1_p * sqrt_gamma, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K0",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_p, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_p * sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1",
+		                       psi_matrix(2, 2, psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_p * sqrt_gamma, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K2",
+		                       psi_matrix(2, 2, psi_new_complex(sqrt_1_p * sqrt_1_gamma, 0.0),
+		                                  psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+		                                  psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator(
+				"K3",
+				psi_matrix(2, 2, psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+				           psi_new_complex(sqrt_1_p * sqrt_gamma, 0.0), psi_new_complex(0.0, 0.0))),
 	};
 
 	return psi_new_noise_channel("GeneralisedAmplitudeDamping", ops, 4, 1);
@@ -182,41 +200,40 @@ struct PsiNoiseChannel psi_generalised_amplitude_damping_channel(double p, doubl
 struct PsiDensityMatrix psi_new_density_matrix(size_t num_qubits)
 {
 	size_t dim = (size_t)1 << num_qubits;
-	struct PsiComplex *data = calloc(dim * dim, sizeof(struct PsiComplex));
+	struct PsiComplex* data = calloc(dim * dim, sizeof(struct PsiComplex));
 	assert(data != NULL);
 	data[0] = psi_new_complex(1.0, 0.0);
 
-	return (struct PsiDensityMatrix)
-	{
+	return (struct PsiDensityMatrix){
 		data,
 		dim,
 		num_qubits,
 	};
 }
 
-struct PsiDensityMatrix psi_new_density_matrix_from_state(const struct PsiComplex *state, size_t len)
+struct PsiDensityMatrix psi_new_density_matrix_from_state(const struct PsiComplex* state,
+                                                          size_t len)
 {
 	size_t dim = len;
 	size_t num_qubits = 0;
 	while (((size_t)1 << num_qubits) < dim)
 		num_qubits++;
 
-	struct PsiComplex *data = malloc(dim * dim * sizeof(struct PsiComplex));
+	struct PsiComplex* data = malloc(dim * dim * sizeof(struct PsiComplex));
 	assert(data != NULL);
 
 	for (size_t i = 0; i < dim; i++)
 		for (size_t j = 0; j < dim; j++)
 			data[i * dim + j] = psi_mul_complex(state[i], psi_conjugate_complex(state[j]));
 
-	return (struct PsiDensityMatrix)
-	{
+	return (struct PsiDensityMatrix){
 		data,
 		dim,
 		num_qubits,
 	};
 }
 
-void psi_free_density_matrix(struct PsiDensityMatrix *dm)
+void psi_free_density_matrix(struct PsiDensityMatrix* dm)
 {
 	free(dm->data);
 	dm->data = NULL;
@@ -230,7 +247,8 @@ struct PsiComplex psi_get_density_matrix(struct PsiDensityMatrix dm, size_t row,
 	return dm.data[row * dm.dim + col];
 }
 
-void psi_set_density_matrix(struct PsiDensityMatrix *dm, size_t row, size_t col, struct PsiComplex value)
+void psi_set_density_matrix(struct PsiDensityMatrix* dm, size_t row, size_t col,
+                            struct PsiComplex value)
 {
 	assert(row < dm->dim && col < dm->dim);
 	dm->data[row * dm->dim + col] = value;
@@ -250,7 +268,8 @@ double psi_purity_density_matrix(struct PsiDensityMatrix dm)
 	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
 	for (size_t i = 0; i < dm.dim; i++)
 		for (size_t j = 0; j < dm.dim; j++)
-			sum = psi_add_complex(sum, psi_mul_complex(dm.data[i * dm.dim + j], dm.data[j * dm.dim + i]));
+			sum = psi_add_complex(
+					sum, psi_mul_complex(dm.data[i * dm.dim + j], dm.data[j * dm.dim + i]));
 
 	return sum.real;
 }
@@ -260,19 +279,20 @@ bool psi_is_pure_density_matrix(struct PsiDensityMatrix dm, double tolerance)
 	return fabs(psi_purity_density_matrix(dm) - 1.0) < tolerance;
 }
 
-void psi_density_matrix_probabilities(struct PsiDensityMatrix dm, double *out)
+void psi_density_matrix_probabilities(struct PsiDensityMatrix dm, double* out)
 {
 	for (size_t i = 0; i < dm.dim; i++)
 		out[i] = dm.data[i * dm.dim + i].real;
 }
 
-void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMatrix gate, const size_t *targets, size_t target_count)
+void psi_apply_unitary_density_matrix(struct PsiDensityMatrix* dm, struct PsiMatrix gate,
+                                      const size_t* targets, size_t target_count)
 {
 	size_t g = target_count;
 	size_t gate_dim = (size_t)1 << g;
 	size_t dim = dm->dim;
 
-	size_t *target_bits = malloc(g * sizeof(size_t));
+	size_t* target_bits = malloc(g * sizeof(size_t));
 	assert(target_bits != NULL || g == 0);
 	for (size_t t = 0; t < g; t++)
 		target_bits[t] = dm->num_qubits - 1 - targets[t];
@@ -281,7 +301,7 @@ void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMat
 	for (size_t t = 0; t < g; t++)
 		non_target_mask &= ~((size_t)1 << target_bits[t]);
 
-	struct PsiComplex *new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+	struct PsiComplex* new_data = calloc(dim * dim, sizeof(struct PsiComplex));
 	assert(new_data != NULL);
 
 	for (size_t i = 0; i < dim; i++)
@@ -314,10 +334,12 @@ void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMat
 					}
 
 					struct PsiComplex u_ik = gate.data[tgt_i * gate_dim + k];
-					struct PsiComplex u_jl_dag = psi_conjugate_complex(gate.data[tgt_j * gate_dim + l]);
+					struct PsiComplex u_jl_dag =
+							psi_conjugate_complex(gate.data[tgt_j * gate_dim + l]);
 					struct PsiComplex rho_kl = dm->data[src_i * dim + src_j];
 
-					sum = psi_add_complex(sum, psi_mul_complex(psi_mul_complex(u_ik, rho_kl), u_jl_dag));
+					sum = psi_add_complex(sum,
+					                      psi_mul_complex(psi_mul_complex(u_ik, rho_kl), u_jl_dag));
 				}
 
 			new_data[i * dim + j] = sum;
@@ -328,14 +350,15 @@ void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMat
 	dm->data = new_data;
 }
 
-void psi_apply_noise_channel(struct PsiDensityMatrix *dm, struct PsiNoiseChannel channel, size_t target)
+void psi_apply_noise_channel(struct PsiDensityMatrix* dm, struct PsiNoiseChannel channel,
+                             size_t target)
 {
 	assert(channel.num_qubits == 1);
 
 	size_t dim = dm->dim;
 	size_t target_bit = dm->num_qubits - 1 - target;
 
-	struct PsiComplex *new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+	struct PsiComplex* new_data = calloc(dim * dim, sizeof(struct PsiComplex));
 	assert(new_data != NULL);
 
 	for (size_t op = 0; op < channel.operator_count; op++)
@@ -355,10 +378,12 @@ void psi_apply_noise_channel(struct PsiDensityMatrix *dm, struct PsiNoiseChannel
 						size_t src_j = (j & ~((size_t)1 << target_bit)) | (kj << target_bit);
 
 						struct PsiComplex k_elem = k.data[i_target * 2 + ki];
-						struct PsiComplex k_dag_elem = psi_conjugate_complex(k.data[j_target * 2 + kj]);
+						struct PsiComplex k_dag_elem =
+								psi_conjugate_complex(k.data[j_target * 2 + kj]);
 						struct PsiComplex rho_elem = dm->data[src_i * dim + src_j];
 
-						struct PsiComplex term = psi_mul_complex(psi_mul_complex(k_elem, rho_elem), k_dag_elem);
+						struct PsiComplex term =
+								psi_mul_complex(psi_mul_complex(k_elem, rho_elem), k_dag_elem);
 						new_data[i * dim + j] = psi_add_complex(new_data[i * dim + j], term);
 					}
 			}
@@ -368,7 +393,8 @@ void psi_apply_noise_channel(struct PsiDensityMatrix *dm, struct PsiNoiseChannel
 	dm->data = new_data;
 }
 
-double psi_measure_probability_density_matrix(struct PsiDensityMatrix dm, size_t qubit, size_t outcome)
+double psi_measure_probability_density_matrix(struct PsiDensityMatrix dm, size_t qubit,
+                                              size_t outcome)
 {
 	size_t target_bit = dm.num_qubits - 1 - qubit;
 	double prob = 0.0;
@@ -380,13 +406,16 @@ double psi_measure_probability_density_matrix(struct PsiDensityMatrix dm, size_t
 	return prob;
 }
 
-double psi_fidelity_density_matrix(struct PsiDensityMatrix dm, const struct PsiComplex *state)
+double psi_fidelity_density_matrix(struct PsiDensityMatrix dm, const struct PsiComplex* state)
 {
 	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
 
 	for (size_t i = 0; i < dm.dim; i++)
 		for (size_t j = 0; j < dm.dim; j++)
-			sum = psi_add_complex(sum, psi_mul_complex(psi_mul_complex(psi_conjugate_complex(state[i]), dm.data[i * dm.dim + j]), state[j]));
+			sum = psi_add_complex(sum,
+			                      psi_mul_complex(psi_mul_complex(psi_conjugate_complex(state[i]),
+			                                                      dm.data[i * dm.dim + j]),
+			                                      state[j]));
 
 	return sum.real;
 }