State vector norm path, more gates, calculation of phase

qasm/run_qasm_on_mtbdd.c

@@ -23,6 +23,7 @@
 #include <qsylvan_simulator_mtbdd.h>
 #include <qsylvan_gates_mtbdd_mpc.h>
 
+#include <sylvan_mpc.h>
 #include "simple_parser.h" // TODO: rename in qsylvan_qasm_parser.h
 
 /**
@@ -133,10 +134,10 @@ static struct argp argp = { options, parse_opt, "<qasm_file>", 0, 0, 0, 0 };
 
 typedef struct stats_s {
 
-    uint64_t applied_gates;     // Number of gates used ?
-    uint64_t final_nodes;       // node count after measurements at end of circuit if any ?
-    uint64_t max_nodes;         // ?
-    uint64_t shots;             // ?
+    uint64_t applied_gates;     // Number of gates used
+    uint64_t final_nodes;       // node count before the measurements
+    uint64_t max_nodes;         // maximum of nodes >= final number of nodes
+    uint64_t shots;             // Number of measurements
     double simulation_time;     // Time of simulation
     double norm;                // Norm L2 of the final state (should be 1.00000), TODO: make mpc type of this ?
     MTBDD final_state;          // State vector MTBDD after simulation
@@ -159,10 +160,9 @@ void fprint_stats(FILE *stream, quantum_circuit_t* circuit)
         fprintf(stream, "  \"state_vector\": [\n");
         for (int k = 0; k < (1<<(circuit->qreg_size)); k++) {
             bool *x = int_to_bitarray(k, circuit->qreg_size, !(circuit->reversed_qubit_order));
-            complex_t c = qmdd_get_amplitude(stats.final_state, x, circuit->qreg_size);
+            MTBDD leaf = mtbdd_getvalue_of_path(stats.final_state, x); // Perhaps reverse qubit sequence on circuit->qreg_size, see gmdd_get_amplitude
             fprintf(stream, "    [\n");
-            fprintf(stream, "      %.16lf,\n", c.r);
-            fprintf(stream, "      %.16lf\n", c.i);
+            mpc_out_str(stream, MPC_BASE_OF_FLOAT, 3, (mpc_ptr)mtbdd_getvalue(leaf), MPC_ROUNDING);
             if (k == (1<<(circuit->qreg_size))-1)
                 fprintf(stream, "    ]\n");
             else
@@ -181,8 +181,8 @@ void fprint_stats(FILE *stream, quantum_circuit_t* circuit)
     fprintf(stream, "    \"seed\": %d,\n", rseed);
     fprintf(stream, "    \"shots\": %" PRIu64 ",\n", stats.shots);
     fprintf(stream, "    \"simulation_time\": %lf,\n", stats.simulation_time);
-    fprintf(stream, "    \"precision\": %d,\n", precision); // Remove
-    fprintf(stream, "    \"rounding\": %d,\n", rounding); // Remove
+    fprintf(stream, "    \"precision\": %d,\n", precision);
+    fprintf(stream, "    \"rounding\": %d,\n", rounding);
     fprintf(stream, "    \"workers\": %d\n", workers);
     fprintf(stream, "  }\n");
     fprintf(stream, "}\n");
@@ -259,13 +259,13 @@ MTBDD apply_gate(MTBDD state, quantum_op_t* gate, int n) // zie master branch
         return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
 
-    else if (strcmp(gate->name, "sx") == 0) { // ???
-        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, S_dd);
+    else if (strcmp(gate->name, "sx") == 0) {
+        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, sqrt_X_dd);
         return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
 
-    else if (strcmp(gate->name, "sxdg") == 0) { // ???
-        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, S_dag_dd);
+    else if (strcmp(gate->name, "sxdg") == 0) {
+        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, sqrt_X_dag_dd);
         return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
 
@@ -274,17 +274,27 @@ MTBDD apply_gate(MTBDD state, quantum_op_t* gate, int n) // zie master branch
         MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, Rx_dd);
         return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
-
-/*    
+
     else if (strcmp(gate->name, "ry") == 0) {
-        return qmdd_gate(state, GATEID_Ry(gate->angle[0]), gate->targets[0]);
+        MTBDD Ry_dd = mtbdd_Ry(gate->angle[0]);
+        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, Ry_dd);
+        return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
+
     else if (strcmp(gate->name, "rz") == 0) {
-        return qmdd_gate(state, GATEID_Rz(gate->angle[0]), gate->targets[0]);
+        MTBDD Rz_dd = mtbdd_Rz(gate->angle[0]);
+        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, Rz_dd);
+        return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
     }
+
     else if (strcmp(gate->name, "p") == 0) {
+        MTBDD P_dd = mtbdd_Phase(gate->angle[0]);
+        MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, P_dd);
+        return mtbdd_matvec_mult(M_dd, state, (1 << n), 0);
         return qmdd_gate(state, GATEID_Phase(gate->angle[0]), gate->targets[0]);
     }
+
+/*    
     else if (strcmp(gate->name, "u2") == 0) {
         fl_t pi_over_2 = flt_acos(0.0);
         return qmdd_gate(state, GATEID_U(pi_over_2, gate->angle[0], gate->angle[1]), gate->targets[0]);

src/qsylvan_gates_mtbdd_mpc.c

@@ -559,20 +559,33 @@ mtbdd_U(double theta, double phi, double lambda)
     mpfr_set_d(mpfr_lambda, lambda, MPC_ROUNDING);
     mpfr_set_d(mpfr_gam, phi + lambda, MPC_ROUNDING);
 
-    mpfr_div_ui(mpfr_theta_2, mpfr_theta, 2, MPC_ROUNDING);  //  theta / 2
+    mpfr_div_ui(mpfr_theta_2, mpfr_theta, 2, MPC_ROUNDING);           // theta / 2
+
+    mpfr_cos(mpfr_cos_theta_2, mpfr_theta_2, MPC_ROUNDING);           // cos(theta/2)
+    mpfr_sin(mpfr_sin_theta_2, mpfr_theta_2, MPC_ROUNDING);           // sin(theta/2)
+    mpfr_neg(mpfr_min_sin_theta_2, mpfr_sin_theta_2, MPC_ROUNDING);   // -sin(theta/2)
+
+    // phi cos / sin
+    mpfr_cos(mpfr_cos_phi, mpfr_phi, MPC_ROUNDING);                   // cos(phi)
+    mpfr_sin(mpfr_sin_phi, mpfr_phi, MPC_ROUNDING);                   // sin(phi)
+
+    // lambda cos / sin
+    mpfr_cos(mpfr_cos_lambda, mpfr_lambda, MPC_ROUNDING);             // cos(lambda)
+    mpfr_sin(mpfr_sin_lambda, mpfr_lambda, MPC_ROUNDING);             // sin(lambda)
+
+    // lambda cos / sin
+    mpfr_cos(mpfr_cos_gam, mpfr_gam, MPC_ROUNDING);                   // cos(gamma)
+    mpfr_sin(mpfr_sin_gam, mpfr_gam, MPC_ROUNDING);                   // sin(gamma)
 
-    mpfr_cos(mpfr_cos_theta_2, mpfr_theta_2, MPC_ROUNDING);  // cos(theta/2)
-    mpfr_sin(mpfr_sin_theta_2, mpfr_theta_2, MPC_ROUNDING);  // sin(theta/2)
-    mpfr_neg(mpfr_min_sin_theta_2, mpfr_sin_theta_2, MPC_ROUNDING); // -sin(theta/2)
 
-    mpc_set_fr_fr(mpc_exp_lambda, mpfr_cos_lambda, mpfr_sin_lambda, MPC_ROUNDING); // cos(lambda) + i sin(lambda)
-    mpc_mul_fr(mpc_exp_lambda_mul_min_sin_theta_2, mpc_exp_lambda, mpfr_min_sin_theta_2, MPC_ROUNDING);
+    mpc_set_fr_fr(mpc_exp_phi, mpfr_cos_phi, mpfr_sin_phi, MPC_ROUNDING);                 // cos(phi) + i sin(phi)
+    mpc_mul_fr(mpc_exp_phi_mul_sin_theta_2, mpc_exp_phi, mpfr_sin_theta_2, MPC_ROUNDING); // (cos(phi) + i sin(phi)) x sin(theta/2)
 
-    mpc_set_fr_fr(mpc_exp_phi, mpfr_cos_phi, mpfr_sin_phi, MPC_ROUNDING); // cos(phi) + i sin(phi)
-    mpc_mul_fr(mpc_exp_phi_mul_sin_theta_2, mpc_exp_phi, mpfr_sin_theta_2, MPC_ROUNDING);
+    mpc_set_fr_fr(mpc_exp_lambda, mpfr_cos_lambda, mpfr_sin_lambda, MPC_ROUNDING);                      // cos(lambda) + i sin(lambda)
+    mpc_mul_fr(mpc_exp_lambda_mul_min_sin_theta_2, mpc_exp_lambda, mpfr_min_sin_theta_2, MPC_ROUNDING); // (cos(lambda) + i sin(lambda)) x -sin(theta/2)
 
-    mpc_set_fr_fr(mpc_exp_gam, mpfr_cos_gam, mpfr_sin_gam, MPC_ROUNDING); // cos(gamma) + i sin(gamma)
-    mpc_mul_fr(mpc_exp_gam_mul_cos_theta_2, mpc_exp_gam, mpfr_cos_theta_2, MPC_ROUNDING);
+    mpc_set_fr_fr(mpc_exp_gam, mpfr_cos_gam, mpfr_sin_gam, MPC_ROUNDING);                 // cos(gamma) + i sin(gamma)
+    mpc_mul_fr(mpc_exp_gam_mul_cos_theta_2, mpc_exp_gam, mpfr_cos_theta_2, MPC_ROUNDING); // (cos(gamma) + i sin(gamma)) x cos(theta/2)
 
     // Convert to MTBDD
 

src/sylvan_mtbdd.c

@@ -5009,3 +5009,30 @@ MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_depth_of_a) // leaf_depth_of_
     return result;
 }
 
+/**
+ * 
+ * 
+ */
+MTBDD // index to leaf
+mtbdd_getvalue_of_path(MTBDD a, bool* path)
+{
+    MTBDD low, high;
+    for (;;) {
+
+        // if the current edge is pointing to the terminal, we're done.
+        if (mtbdd_isleaf(a)) 
+            return a;
+
+        // now we need to choose low or high edge of next node
+        mtbddnode_t node = MTBDD_GETNODE(a);
+        BDDVAR var = mtbddnode_getvariable(node);
+        low = mtbddnode_getlow(node);
+        high = mtbddnode_getlow(node);
+
+        // Condition low/high choice on basis state vector[var]
+        a = (path[var] == 0) ? low : high;
+    }
+
+    return MTBDD_ZERO;
+}
+

src/sylvan_mtbdd.h

@@ -48,6 +48,8 @@
 extern "C" {
 #endif /* __cplusplus */
 
+#include <stdbool.h>
+
 /**
  * A MTBDD node has a 64-bit unsigned integer type. The low 40 bits are an 
  * index into an unique table. The highest 1 bit is the complement edge, indicating negation.
@@ -1344,6 +1346,13 @@ MTBDD mtbdd_mat_tensor_prod(MTBDD M1, MTBDD M2, int n);
  */
 MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_depth_of_a);
 
+/**
+ * 
+ * 
+ */
+MTBDD // index to leaf
+mtbdd_getvalue_of_path(MTBDD a, bool* path);
+
 
 #ifdef __cplusplus
 }