feat: add matrix, vector, and circuit display + fmt

1 files changed, 266 insertions, 266 deletions

diff --git a/src/core/noise.c b/src/core/noise.c
index 5f141c8..0f73c48 100644
--- a/src/core/noise.c
+++ b/src/core/noise.c
@@ -7,415 +7,415 @@
 
 struct PsiKrausOperator psi_new_kraus_operator(const char* name, struct PsiMatrix matrix)
 {
-	return (struct PsiKrausOperator){
-		matrix,
-		name,
-	};
+    return (struct PsiKrausOperator){
+        matrix,
+        name,
+    };
 }
 
 void psi_free_kraus_operator(struct PsiKrausOperator* op)
 {
-	psi_free_matrix(&op->matrix);
+    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 PsiKrausOperator* owned = malloc(count * sizeof(struct PsiKrausOperator));
-	assert(owned != NULL || count == 0);
+    struct PsiKrausOperator* owned = malloc(count * sizeof(struct PsiKrausOperator));
+    assert(owned != NULL || count == 0);
 
-	if (count > 0)
-		memcpy(owned, operators, count * sizeof(struct PsiKrausOperator));
+    if (count > 0)
+        memcpy(owned, operators, count * sizeof(struct PsiKrausOperator));
 
-	return (struct PsiNoiseChannel){
-		name,
-		owned,
-		count,
-		num_qubits,
-	};
+    return (struct PsiNoiseChannel){
+        name,
+        owned,
+        count,
+        num_qubits,
+    };
 }
 
 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);
+    for (size_t i = 0; i < channel->operator_count; i++)
+        psi_free_matrix(&channel->operators[i].matrix);
 
-	free(channel->operators);
-	channel->operators = NULL;
-	channel->operator_count = 0;
+    free(channel->operators);
+    channel->operators = NULL;
+    channel->operator_count = 0;
 }
 
 struct PsiNoiseChannel psi_depolarising_channel(double p)
 {
-	double sqrt_1_p = sqrt(1.0 - p);
-	double sqrt_p3 = sqrt(p / 3.0);
+    double sqrt_1_p = sqrt(1.0 - 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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("Depolarising", ops, 4, 1);
+    return psi_new_noise_channel("Depolarising", ops, 4, 1);
 }
 
 struct PsiNoiseChannel psi_amplitude_damping_channel(double gamma)
 {
-	double sqrt_gamma = sqrt(gamma);
-	double sqrt_1_gamma = sqrt(1.0 - gamma);
+    double sqrt_gamma = sqrt(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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("AmplitudeDamping", ops, 2, 1);
+    return psi_new_noise_channel("AmplitudeDamping", ops, 2, 1);
 }
 
 struct PsiNoiseChannel psi_phase_damping_channel(double gamma)
 {
-	double sqrt_gamma = sqrt(gamma);
-	double sqrt_1_gamma = sqrt(1.0 - gamma);
+    double sqrt_gamma = sqrt(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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("PhaseDamping", ops, 2, 1);
+    return psi_new_noise_channel("PhaseDamping", ops, 2, 1);
 }
 
 struct PsiNoiseChannel psi_bit_flip_channel(double p)
 {
-	double sqrt_1_p = sqrt(1.0 - p);
-	double sqrt_p = sqrt(p);
+    double sqrt_1_p = sqrt(1.0 - 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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("BitFlip", ops, 2, 1);
+    return psi_new_noise_channel("BitFlip", ops, 2, 1);
 }
 
 struct PsiNoiseChannel psi_phase_flip_channel(double p)
 {
-	double sqrt_1_p = sqrt(1.0 - p);
-	double sqrt_p = sqrt(p);
+    double sqrt_1_p = sqrt(1.0 - 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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("PhaseFlip", ops, 2, 1);
+    return psi_new_noise_channel("PhaseFlip", ops, 2, 1);
 }
 
 struct PsiNoiseChannel psi_bit_phase_flip_channel(double p)
 {
-	double sqrt_1_p = sqrt(1.0 - p);
-	double sqrt_p = sqrt(p);
+    double sqrt_1_p = sqrt(1.0 - 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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("BitPhaseFlip", ops, 2, 1);
+    return psi_new_noise_channel("BitPhaseFlip", ops, 2, 1);
 }
 
 struct PsiNoiseChannel psi_generalised_amplitude_damping_channel(double p, double gamma)
 {
-	double sqrt_p = sqrt(p);
-	double sqrt_1_p = sqrt(1.0 - p);
-	double sqrt_gamma = sqrt(gamma);
-	double sqrt_1_gamma = sqrt(1.0 - gamma);
+    double sqrt_p = sqrt(p);
+    double sqrt_1_p = sqrt(1.0 - p);
+    double sqrt_gamma = sqrt(gamma);
+    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))),
-	};
+    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))),
+    };
 
-	return psi_new_noise_channel("GeneralisedAmplitudeDamping", ops, 4, 1);
+    return psi_new_noise_channel("GeneralisedAmplitudeDamping", ops, 4, 1);
 }
 
 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));
-	assert(data != NULL);
-	data[0] = psi_new_complex(1.0, 0.0);
+    size_t dim = (size_t)1 << num_qubits;
+    struct PsiComplex* data = calloc(dim * dim, sizeof(struct PsiComplex));
+    assert(data != NULL);
+    data[0] = psi_new_complex(1.0, 0.0);
 
-	return (struct PsiDensityMatrix){
-		data,
-		dim,
-		num_qubits,
-	};
+    return (struct PsiDensityMatrix){
+        data,
+        dim,
+        num_qubits,
+    };
 }
 
 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++;
+    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));
-	assert(data != NULL);
+    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]));
+    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){
-		data,
-		dim,
-		num_qubits,
-	};
+    return (struct PsiDensityMatrix){
+        data,
+        dim,
+        num_qubits,
+    };
 }
 
 void psi_free_density_matrix(struct PsiDensityMatrix* dm)
 {
-	free(dm->data);
-	dm->data = NULL;
-	dm->dim = 0;
-	dm->num_qubits = 0;
+    free(dm->data);
+    dm->data = NULL;
+    dm->dim = 0;
+    dm->num_qubits = 0;
 }
 
 struct PsiComplex psi_get_density_matrix(struct PsiDensityMatrix dm, size_t row, size_t col)
 {
-	assert(row < dm.dim && col < dm.dim);
-	return dm.data[row * dm.dim + col];
+    assert(row < dm.dim && col < dm.dim);
+    return dm.data[row * dm.dim + col];
 }
 
 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;
+    assert(row < dm->dim && col < dm->dim);
+    dm->data[row * dm->dim + col] = value;
 }
 
 struct PsiComplex psi_trace_density_matrix(struct PsiDensityMatrix dm)
 {
-	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
-	for (size_t i = 0; i < dm.dim; i++)
-		sum = psi_add_complex(sum, dm.data[i * dm.dim + i]);
+    struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+    for (size_t i = 0; i < dm.dim; i++)
+        sum = psi_add_complex(sum, dm.data[i * dm.dim + i]);
 
-	return sum;
+    return sum;
 }
 
 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]));
+    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]));
 
-	return sum.real;
+    return sum.real;
 }
 
 bool psi_is_pure_density_matrix(struct PsiDensityMatrix dm, double tolerance)
 {
-	return fabs(psi_purity_density_matrix(dm) - 1.0) < tolerance;
+    return fabs(psi_purity_density_matrix(dm) - 1.0) < tolerance;
 }
 
 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;
+    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)
 {
-	size_t g = target_count;
-	size_t gate_dim = (size_t)1 << g;
-	size_t dim = dm->dim;
+    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));
-	assert(target_bits != NULL || g == 0);
-	for (size_t t = 0; t < g; t++)
-		target_bits[t] = dm->num_qubits - 1 - targets[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];
 
-	size_t non_target_mask = dim - 1;
-	for (size_t t = 0; t < g; t++)
-		non_target_mask &= ~((size_t)1 << target_bits[t]);
+    size_t non_target_mask = dim - 1;
+    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));
-	assert(new_data != NULL);
+    struct PsiComplex* new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+    assert(new_data != NULL);
 
-	for (size_t i = 0; i < dim; i++)
-		for (size_t j = 0; j < dim; j++)
-		{
-			struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+    for (size_t i = 0; i < dim; i++)
+        for (size_t j = 0; j < dim; j++)
+        {
+            struct PsiComplex sum = psi_new_complex(0.0, 0.0);
 
-			for (size_t k = 0; k < gate_dim; k++)
-				for (size_t l = 0; l < gate_dim; l++)
-				{
-					size_t src_i = i & non_target_mask;
-					size_t src_j = j & non_target_mask;
+            for (size_t k = 0; k < gate_dim; k++)
+                for (size_t l = 0; l < gate_dim; l++)
+                {
+                    size_t src_i = i & non_target_mask;
+                    size_t src_j = j & non_target_mask;
 
-					for (size_t idx = 0; idx < g; idx++)
-					{
-						if ((k >> (g - 1 - idx)) & 1)
-							src_i |= (size_t)1 << target_bits[idx];
-						if ((l >> (g - 1 - idx)) & 1)
-							src_j |= (size_t)1 << target_bits[idx];
-					}
+                    for (size_t idx = 0; idx < g; idx++)
+                    {
+                        if ((k >> (g - 1 - idx)) & 1)
+                            src_i |= (size_t)1 << target_bits[idx];
+                        if ((l >> (g - 1 - idx)) & 1)
+                            src_j |= (size_t)1 << target_bits[idx];
+                    }
 
-					size_t tgt_i = 0;
-					size_t tgt_j = 0;
-					for (size_t idx = 0; idx < g; idx++)
-					{
-						if ((i >> target_bits[idx]) & 1)
-							tgt_i |= (size_t)1 << (g - 1 - idx);
-						if ((j >> target_bits[idx]) & 1)
-							tgt_j |= (size_t)1 << (g - 1 - idx);
-					}
+                    size_t tgt_i = 0;
+                    size_t tgt_j = 0;
+                    for (size_t idx = 0; idx < g; idx++)
+                    {
+                        if ((i >> target_bits[idx]) & 1)
+                            tgt_i |= (size_t)1 << (g - 1 - idx);
+                        if ((j >> target_bits[idx]) & 1)
+                            tgt_j |= (size_t)1 << (g - 1 - idx);
+                    }
 
-					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 rho_kl = dm->data[src_i * dim + src_j];
+                    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 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;
-		}
+            new_data[i * dim + j] = sum;
+        }
 
-	free(target_bits);
-	free(dm->data);
-	dm->data = new_data;
+    free(target_bits);
+    free(dm->data);
+    dm->data = new_data;
 }
 
 void psi_apply_noise_channel(struct PsiDensityMatrix* dm, struct PsiNoiseChannel channel,
                              size_t target)
 {
-	assert(channel.num_qubits == 1);
+    assert(channel.num_qubits == 1);
 
-	size_t dim = dm->dim;
-	size_t target_bit = dm->num_qubits - 1 - target;
+    size_t dim = dm->dim;
+    size_t target_bit = dm->num_qubits - 1 - target;
 
-	struct PsiComplex* new_data = calloc(dim * dim, sizeof(struct PsiComplex));
-	assert(new_data != NULL);
+    struct PsiComplex* new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+    assert(new_data != NULL);
 
-	for (size_t op = 0; op < channel.operator_count; op++)
-	{
-		struct PsiMatrix k = channel.operators[op].matrix;
+    for (size_t op = 0; op < channel.operator_count; op++)
+    {
+        struct PsiMatrix k = channel.operators[op].matrix;
 
-		for (size_t i = 0; i < dim; i++)
-			for (size_t j = 0; j < dim; j++)
-			{
-				size_t i_target = (i >> target_bit) & 1;
-				size_t j_target = (j >> target_bit) & 1;
+        for (size_t i = 0; i < dim; i++)
+            for (size_t j = 0; j < dim; j++)
+            {
+                size_t i_target = (i >> target_bit) & 1;
+                size_t j_target = (j >> target_bit) & 1;
 
-				for (size_t ki = 0; ki < 2; ki++)
-					for (size_t kj = 0; kj < 2; kj++)
-					{
-						size_t src_i = (i & ~((size_t)1 << target_bit)) | (ki << target_bit);
-						size_t src_j = (j & ~((size_t)1 << target_bit)) | (kj << target_bit);
+                for (size_t ki = 0; ki < 2; ki++)
+                    for (size_t kj = 0; kj < 2; kj++)
+                    {
+                        size_t src_i = (i & ~((size_t)1 << target_bit)) | (ki << target_bit);
+                        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 rho_elem = dm->data[src_i * dim + src_j];
+                        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 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);
-						new_data[i * dim + j] = psi_add_complex(new_data[i * dim + j], term);
-					}
-			}
-	}
+                        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);
+                    }
+            }
+    }
 
-	free(dm->data);
-	dm->data = new_data;
+    free(dm->data);
+    dm->data = new_data;
 }
 
 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;
+    size_t target_bit = dm.num_qubits - 1 - qubit;
+    double prob = 0.0;
 
-	for (size_t i = 0; i < dm.dim; i++)
-		if (((i >> target_bit) & 1) == outcome)
-			prob += dm.data[i * dm.dim + i].real;
+    for (size_t i = 0; i < dm.dim; i++)
+        if (((i >> target_bit) & 1) == outcome)
+            prob += dm.data[i * dm.dim + i].real;
 
-	return prob;
+    return prob;
 }
 
 double psi_fidelity_density_matrix(struct PsiDensityMatrix dm, const struct PsiComplex* state)
 {
-	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+    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]));
+    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]));
 
-	return sum.real;
+    return sum.real;
 }