diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 1ac660a..17080f5 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -35,4 +35,4 @@ add_example(test_algs test_algs.c) target_sources(test_algs PRIVATE ${ALGORITHM_EXAMPLES}) add_executable(circuit_equivalence circuit_equivalence.c) -target_link_libraries(circuit_equivalence qsylvan qsylvan_qasm_parser) +target_link_libraries(circuit_equivalence qsylvan qsylvan_qasm_parser_qmdd) diff --git a/qasm/CMakeLists.txt b/qasm/CMakeLists.txt index 341adae..07ec4bf 100644 --- a/qasm/CMakeLists.txt +++ b/qasm/CMakeLists.txt @@ -9,14 +9,42 @@ # allow the parser to be included in other places as a library # (without compiling eval_expr.cpp again) -add_library(qsylvan_qasm_parser) -target_sources(qsylvan_qasm_parser + +#add_library(qsylvan_qasm_parser) +#target_sources(qsylvan_qasm_parser +# PRIVATE +# simple_parser.cpp +# parse_math/eval_expr.cpp +# PUBLIC +# simple_parser.h) +# +#add_executable(sim_qasm sim_qasm.c) +#target_link_libraries(sim_qasm qsylvan qsylvan_qasm_parser) + + +# QASM simulator with qmddd + +add_library(qsylvan_qasm_parser_qmdd) +target_sources(qsylvan_qasm_parser_qmdd PRIVATE simple_parser.cpp parse_math/eval_expr.cpp PUBLIC simple_parser.h) -add_executable(sim_qasm sim_qasm.c) -target_link_libraries(sim_qasm qsylvan qsylvan_qasm_parser) +add_executable(run_qasm_on_qmdd run_qasm_on_qmdd.c) +target_link_libraries(run_qasm_on_qmdd qsylvan qsylvan_qasm_parser_qmdd) + + +# QASM simulator with mtbdd + +add_library(qsylvan_qasm_parser_mtbdd) +target_sources(qsylvan_qasm_parser_mtbdd + PRIVATE + simple_parser.cpp + parse_math/eval_expr.cpp + PUBLIC + simple_parser.h) +add_executable(run_qasm_on_mtbdd run_qasm_on_mtbdd.c) +target_link_libraries(run_qasm_on_mtbdd qsylvan qsylvan_qasm_parser_mtbdd) diff --git a/qasm/QASM_to_circuit.h b/qasm/QASM_to_circuit.h index 7456484..bcb8495 100644 --- a/qasm/QASM_to_circuit.h +++ b/qasm/QASM_to_circuit.h @@ -182,7 +182,7 @@ BDDVAR partition(BDDVAR* qubits, int low, int high); * PARAMETERS: * - qubits: the circuit in which to store the barrier column * - low: the first index of the array - * - high: the last index of the array + * - high: the last index of the array */ void sort(BDDVAR* qubits, int low, int high); diff --git a/qasm/circuits/qft_n4.qasm b/qasm/circuits/qft_n4.qasm index 0b0fe08..2067ed9 100644 --- a/qasm/circuits/qft_n4.qasm +++ b/qasm/circuits/qft_n4.qasm @@ -1,18 +1,18 @@ -// quantum Fourier transform -OPENQASM 2.0; -include "qelib1.inc"; -qreg q[4]; -creg c[4]; -x q[0]; -x q[2]; -barrier q; -h q[0]; -cu1(pi/2) q[1],q[0]; -h q[1]; -cu1(pi/4) q[2],q[0]; -cu1(pi/2) q[2],q[1]; -h q[2]; -cu1(pi/8) q[3],q[0]; -cu1(pi/4) q[3],q[1]; -cu1(pi/2) q[3],q[2]; -h q[3]; +// quantum Fourier transform +OPENQASM 2.0; +include "qelib1.inc"; +qreg q[4]; +creg c[4]; +x q[0]; +x q[2]; +barrier q; +h q[0]; +cu1(pi/2) q[1],q[0]; +h q[1]; +cu1(pi/4) q[2],q[0]; +cu1(pi/2) q[2],q[1]; +h q[2]; +cu1(pi/8) q[3],q[0]; +cu1(pi/4) q[3],q[1]; +cu1(pi/2) q[3],q[2]; +h q[3]; diff --git a/qasm/circuits/vqe_n4.qasm b/qasm/circuits/vqe_n4.qasm index 3e5350a..b922d4e 100644 --- a/qasm/circuits/vqe_n4.qasm +++ b/qasm/circuits/vqe_n4.qasm @@ -1,98 +1,98 @@ -OPENQASM 2.0; -include "qelib1.inc"; -qreg q[4]; -creg meas[4]; -sx q[0]; -rz(5.0300511584448) q[0]; -sx q[0]; -rz(3*pi) q[0]; -rz(2.61519568487349) q[0]; -sx q[1]; -rz(5.36606826220224) q[1]; -sx q[1]; -rz(3*pi) q[1]; -rz(0.667082990176662) q[1]; -cx q[0],q[1]; -sx q[0]; -rz(4.09739784898316) q[0]; -sx q[0]; -rz(3*pi) q[0]; -rz(1.92219255913748) q[0]; -sx q[2]; -rz(3.20626074964735) q[2]; -sx q[2]; -rz(3*pi) q[2]; -rz(0.571219981217032) q[2]; -cx q[1],q[2]; -sx q[1]; -rz(4.79016360412963) q[1]; -sx q[1]; -rz(3*pi) q[1]; -rz(0.438232887845333) q[1]; -cx q[0],q[1]; -sx q[0]; -rz(4.57437876552885) q[0]; -sx q[0]; -rz(3*pi) q[0]; -rz(1.86112525741656) q[0]; -sx q[3]; -rz(6.18865431978628) q[3]; -sx q[3]; -rz(3*pi) q[3]; -rz(0.576182260790784) q[3]; -cx q[2],q[3]; -sx q[2]; -rz(4.49858795091057) q[2]; -sx q[2]; -rz(3*pi) q[2]; -rz(0.917799481623813) q[2]; -cx q[1],q[2]; -sx q[1]; -rz(5.60829568567739) q[1]; -sx q[1]; -rz(3*pi) q[1]; -rz(0.145928275357359) q[1]; -cx q[0],q[1]; -sx q[0]; -rz(3.34595826021666) q[0]; -sx q[0]; -rz(3*pi) q[0]; -rz(0.956972379417358) q[0]; -sx q[3]; -rz(4.05651598094723) q[3]; -sx q[3]; -rz(3*pi) q[3]; -rz(1.15095967544708) q[3]; -cx q[2],q[3]; -sx q[2]; -rz(3.76888634073298) q[2]; -sx q[2]; -rz(3*pi) q[2]; -rz(1.90865844345986) q[2]; -cx q[1],q[2]; -sx q[1]; -rz(6.12260448652247) q[1]; -sx q[1]; -rz(3*pi) q[1]; -rz(0.30684599582304) q[1]; -sx q[3]; -rz(4.75710778753287) q[3]; -sx q[3]; -rz(3*pi) q[3]; -rz(0.535717334235832) q[3]; -cx q[2],q[3]; -sx q[2]; -rz(6.17521515476781) q[2]; -sx q[2]; -rz(3*pi) q[2]; -rz(2.1495814494341) q[2]; -sx q[3]; -rz(5.68124782361394) q[3]; -sx q[3]; -rz(3*pi) q[3]; -rz(1.38277984079156) q[3]; -barrier q[0],q[1],q[2],q[3]; -measure q[0] -> meas[0]; -measure q[1] -> meas[1]; -measure q[2] -> meas[2]; -measure q[3] -> meas[3]; +OPENQASM 2.0; +include "qelib1.inc"; +qreg q[4]; +creg meas[4]; +sx q[0]; +rz(5.0300511584448) q[0]; +sx q[0]; +rz(3*pi) q[0]; +rz(2.61519568487349) q[0]; +sx q[1]; +rz(5.36606826220224) q[1]; +sx q[1]; +rz(3*pi) q[1]; +rz(0.667082990176662) q[1]; +cx q[0],q[1]; +sx q[0]; +rz(4.09739784898316) q[0]; +sx q[0]; +rz(3*pi) q[0]; +rz(1.92219255913748) q[0]; +sx q[2]; +rz(3.20626074964735) q[2]; +sx q[2]; +rz(3*pi) q[2]; +rz(0.571219981217032) q[2]; +cx q[1],q[2]; +sx q[1]; +rz(4.79016360412963) q[1]; +sx q[1]; +rz(3*pi) q[1]; +rz(0.438232887845333) q[1]; +cx q[0],q[1]; +sx q[0]; +rz(4.57437876552885) q[0]; +sx q[0]; +rz(3*pi) q[0]; +rz(1.86112525741656) q[0]; +sx q[3]; +rz(6.18865431978628) q[3]; +sx q[3]; +rz(3*pi) q[3]; +rz(0.576182260790784) q[3]; +cx q[2],q[3]; +sx q[2]; +rz(4.49858795091057) q[2]; +sx q[2]; +rz(3*pi) q[2]; +rz(0.917799481623813) q[2]; +cx q[1],q[2]; +sx q[1]; +rz(5.60829568567739) q[1]; +sx q[1]; +rz(3*pi) q[1]; +rz(0.145928275357359) q[1]; +cx q[0],q[1]; +sx q[0]; +rz(3.34595826021666) q[0]; +sx q[0]; +rz(3*pi) q[0]; +rz(0.956972379417358) q[0]; +sx q[3]; +rz(4.05651598094723) q[3]; +sx q[3]; +rz(3*pi) q[3]; +rz(1.15095967544708) q[3]; +cx q[2],q[3]; +sx q[2]; +rz(3.76888634073298) q[2]; +sx q[2]; +rz(3*pi) q[2]; +rz(1.90865844345986) q[2]; +cx q[1],q[2]; +sx q[1]; +rz(6.12260448652247) q[1]; +sx q[1]; +rz(3*pi) q[1]; +rz(0.30684599582304) q[1]; +sx q[3]; +rz(4.75710778753287) q[3]; +sx q[3]; +rz(3*pi) q[3]; +rz(0.535717334235832) q[3]; +cx q[2],q[3]; +sx q[2]; +rz(6.17521515476781) q[2]; +sx q[2]; +rz(3*pi) q[2]; +rz(2.1495814494341) q[2]; +sx q[3]; +rz(5.68124782361394) q[3]; +sx q[3]; +rz(3*pi) q[3]; +rz(1.38277984079156) q[3]; +barrier q[0],q[1],q[2],q[3]; +measure q[0] -> meas[0]; +measure q[1] -> meas[1]; +measure q[2] -> meas[2]; +measure q[3] -> meas[3]; diff --git a/qasm/run_qasm_on_mtbdd.c b/qasm/run_qasm_on_mtbdd.c new file mode 100644 index 0000000..4fdce2f --- /dev/null +++ b/qasm/run_qasm_on_mtbdd.c @@ -0,0 +1,503 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +#include +#include +#include + +#include +#include +#include + +#include +#include "simple_parser.h" // TODO: rename in qsylvan_qasm_parser.h + +/** + * + * Arguments variables of command line interface (configured via argp). + * + */ +static int workers = 1; // Number of threads running on separate CPU core +static int rseed = 0; +static int precision = 128; +static int rounding = 0; + +static bool count_nodes = false; +static bool output_vector = false; + +static size_t min_tablesize = 1LL<<25; +static size_t max_tablesize = 1LL<<25; +static size_t min_cachesize = 1LL<<16; +static size_t max_cachesize = 1LL<<16; + +static char* qasm_inputfile = NULL; +static char* json_outputfile = NULL; + +/** + * + * Command Line Interface argument help list. + * + */ +static struct argp_option options[] = +{ + {"workers", 'w', "", 0, "Number of workers/threads (default=1)", 0}, + {"rseed", 'r', "", 0, "Set random seed as integer", 0}, + {"precision", 'p', "", 0, "Precision of mantissa multiprecision complex float in bits (default=128)", 0}, + {"rounding", 'o', "<...>", 0, "Rounding strategy", 0}, + {"json", 'j', "", 0, "Write statistics in json format to ", 0}, + {"count-nodes", 'c', 0, 0, "Track maximum number of nodes", 0}, + {"state-vector", 'v', 0, 0, "Show the complete state vector after simulation", 0}, + {0, 0, 0, 0, 0, 0} +}; + +/** + * + * Command Line Interface argument processing. + * + */ + +static error_t +parse_opt(int key, char *arg, struct argp_state *state) +{ + switch (key) { + + case 'w': + workers = atoi(arg); + break; + + case 'r': + rseed = atoi(arg); + break; + + case 'p': // Precision + precision = atoi(arg); // ASCII to integer, TODO: validate value + break; + + case 'o': // Rounding + rounding = atoi(arg); // ASCII to float, TODO: validate value + break; + + case 'j': + json_outputfile = arg; + break; + + case 'c': + count_nodes = true; + break; + + case 'v': + output_vector = true; + break; + + case ARGP_KEY_ARG: + if (state->arg_num >= 1) argp_usage(state); + qasm_inputfile = arg; + break; + + case ARGP_KEY_END: + if (state->arg_num < 1) argp_usage(state); + break; + + default: + return ARGP_ERR_UNKNOWN; + } + return 0; +} +static struct argp argp = { options, parse_opt, "", 0, 0, 0, 0 }; + + +/** + * + * Statistics of the simulation. + * + * Store statistical info in struct below. + * + * NOTE: If there are only measurements at the end of the circuit, 'final_nodes' + * and 'norm' will contain the node count and the norm of the state MTBDD before + * the measurements. + * + */ + +typedef struct stats_s { + + 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 + +} stats_t; + +stats_t stats; + +/** + * Print the statistics after the simulation the given file. + */ +void fprint_stats(FILE *stream, quantum_circuit_t* circuit) +{ + fprintf(stream, "{\n"); + fprintf(stream, " \"measurement_results\": {\n"); + fprintf(stream, " \""); fprint_creg(stream, circuit); fprintf(stream, "\": 1\n"); + fprintf(stream, " },\n"); + if (output_vector) + { + 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)); + 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"); + 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 + fprintf(stream, " ],\n"); + free(x); + } + fprintf(stream, " ],\n"); + } + fprintf(stream, " \"statistics\": {\n"); + fprintf(stream, " \"applied_gates\": %" PRIu64 ",\n", stats.applied_gates); + fprintf(stream, " \"benchmark\": \"%s\",\n", circuit->name); + fprintf(stream, " \"final_nodes\": %" PRIu64 ",\n", stats.final_nodes); + fprintf(stream, " \"max_nodes\": %" PRIu64 ",\n", stats.max_nodes); + fprintf(stream, " \"n_qubits\": %d,\n", circuit->qreg_size); + fprintf(stream, " \"norm\": %.5e,\n", stats.norm); + 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); + fprintf(stream, " \"rounding\": %d,\n", rounding); + fprintf(stream, " \"workers\": %d\n", workers); + fprintf(stream, " }\n"); + fprintf(stream, "}\n"); +} + +static double +wctime() +{ + struct timeval tv; + gettimeofday(&tv, NULL); + return (tv.tv_sec + 1E-6 * tv.tv_usec); +} + +/** + * Here we match the QASM name of a gate with the MTBDD of the corresponding gate. + * + * If the gate is supported the state vector |fi> will be applied to the matrix of the gate: + * + * M = (I (x) ... (x) I (x) G (x) I (x) ... (x) I) |fi> + * + * with I the 2x2 identity matrix, G the 2x2 choosen gate and |fi> the state vector. + * + * The function returns the corresponding MTBDD of M. + * + * The type of the matrix elements are complex numbers based on the GNU multiprecision complex library mpc. + * + */ +MTBDD apply_gate(MTBDD state, quantum_op_t* gate, int n) // zie master branch +{ + stats.applied_gates++; + + if (strcmp(gate->name, "id") == 0) { + stats.applied_gates--; + return state; + } + + else if (strcmp(gate->name, "x") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, X_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "y") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, Y_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "z") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, Z_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "h") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, H_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "s") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, S_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "sdg") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, S_dag_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "t") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, T_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "tdg") == 0) { + MTBDD M_dd = mtbdd_create_single_gate_for_qubits_mpc(n, gate->targets[0], I_dd, T_dag_dd); + 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, 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, sqrt_X_dag_dd); + return mtbdd_matvec_mult(M_dd, state, (1 << n), 0); + } + + else if (strcmp(gate->name, "rx") == 0) { + MTBDD Rx_dd = mtbdd_Rx(gate->angle[0]); + 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) { + 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) { + 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]); + } + else if (strcmp(gate->name, "u") == 0) { + return qmdd_gate(state, GATEID_U(gate->angle[0], gate->angle[1], gate->angle[2]), gate->targets[0]); + } + else if (strcmp(gate->name, "cx") == 0) { + return qmdd_cgate(state, GATEID_X, gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "cy") == 0) { + return qmdd_cgate(state, GATEID_Y, gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "cz") == 0) { + return qmdd_cgate(state, GATEID_Z, gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "ch") == 0) { + return qmdd_cgate(state, GATEID_H, gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "csx") == 0) { + return qmdd_cgate(state, GATEID_sqrtX, gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "crx") == 0) { + return qmdd_cgate(state, GATEID_Rx(gate->angle[0]), gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "cry") == 0) { + return qmdd_cgate(state, GATEID_Ry(gate->angle[0]), gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "crz") == 0) { + return qmdd_cgate(state, GATEID_Rz(gate->angle[0]), gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "cp") == 0) { + return qmdd_cgate(state, GATEID_Phase(gate->angle[0]), gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "cu") == 0) { + return qmdd_cgate(state, GATEID_U(gate->angle[0], gate->angle[1], gate->angle[2]), gate->ctrls[0], gate->targets[0]); + } + else if (strcmp(gate->name, "ccx") == 0) { + return qmdd_cgate2(state, GATEID_X, gate->ctrls[0], gate->ctrls[1], gate->targets[0]); + } + else if (strcmp(gate->name, "c3x") == 0) { + return qmdd_cgate3(state, GATEID_X, gate->ctrls[0], gate->ctrls[1], gate->ctrls[2], gate->targets[0]); + } + else if (strcmp(gate->name, "c3sx") == 0) { + return qmdd_cgate3(state, GATEID_sqrtX, gate->ctrls[0], gate->ctrls[1], gate->ctrls[2], gate->targets[0]); + } + else if (strcmp(gate->name, "swap") == 0) { + // no native SWAP gates in Q-Sylvan + stats.applied_gates += 4; + return qmdd_circuit_swap(state, gate->targets[0], gate->targets[1]); + } + else if (strcmp(gate->name, "cswap") == 0) { + // no native CSWAP gates in Q-Sylvan + stats.applied_gates += 4; + // CCNOT + state = qmdd_cgate2(state, GATEID_X, gate->ctrls[0], gate->targets[0], gate->targets[1]); + // upside down CCNOT (equivalent) + state = qmdd_cgate(state, GATEID_H, gate->ctrls[0], gate->targets[0]); + state = qmdd_cgate2(state, GATEID_Z, gate->ctrls[0], gate->targets[0], gate->targets[1]); + state = qmdd_cgate(state, GATEID_H, gate->ctrls[0], gate->targets[0]); + // CCNOT + state = qmdd_cgate2(state, GATEID_X, gate->ctrls[0], gate->targets[0], gate->targets[1]); + return state; + } + else if (strcmp(gate->name, "rccx") == 0) { + // no native RCCX (simplified Toffoli) gates in Q-Sylvan + stats.applied_gates += 3; + state = qmdd_cgate2(state, GATEID_X, gate->ctrls[0], gate->ctrls[1], gate->targets[0]); + state = qmdd_gate(state, GATEID_X, gate->ctrls[1]); + state = qmdd_cgate2(state, GATEID_Z, gate->ctrls[0], gate->ctrls[1], gate->targets[0]); + state = qmdd_gate(state, GATEID_X, gate->ctrls[1]); + return state; + } + else if (strcmp(gate->name, "rzz") == 0 ) { + // no native RZZ gates in Q-Sylvan + stats.applied_gates += 2; + state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1]); + state = qmdd_gate(state, GATEID_Phase(gate->angle[0]), gate->targets[1]); + state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1]); + return state; + } + else if (strcmp(gate->name, "rxx") == 0) { + // no native RXX gates in Q-Sylvan + fl_t pi = flt_acos(0.0) * 2; + stats.applied_gates += 6; + state = qmdd_gate(state, GATEID_U(pi/2.0, gate->angle[0], 0), gate->targets[0]); + state = qmdd_gate(state, GATEID_H, gate->targets[1]); + state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1]); + state = qmdd_gate(state, GATEID_Phase(-(gate->angle[0])), gate->targets[1]); + state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1]); + state = qmdd_gate(state, GATEID_H, gate->targets[1]); + state = qmdd_gate(state, GATEID_U(pi/2.0, -pi, pi-gate->angle[0]), gate->targets[0]); + return state; + } +*/ + + else { + fprintf(stderr, "Gate '%s' currently unsupported\n", gate->name); + return state; + } +} + +/** + * + * Calculate the final state after observation (measurement). + * + */ +/* TODO: convert to MTBDD +MTBDD measure(QMDD state, quantum_op_t *meas, quantum_circuit_t* circuit) +{ + double p; + int m; + printf("measure qubit %d, store result in creg[%d]\n", meas->targets[0], meas->meas_dest); + +// qmdd_measure_qubit(state, meas->targets[0], circuit->qreg_size, &m, &p); + + circuit->creg[meas->meas_dest] = m; + return state; +} +*/ + +/** + * + * Simulate the circuit as parsed. + * + */ +void simulate_circuit(quantum_circuit_t* circuit) +{ + double t_start = wctime(); + MTBDD state = mtbdd_create_all_zero_state_mpc(circuit->qreg_size); + quantum_op_t *op = circuit->operations; + while (op != NULL) { + if (op->type == op_gate) { + state = apply_gate(state, op, circuit->qreg_size); // , n); + } + else if (op->type == op_measurement) { + + assert(false && "Measurements not yet supported"); + fprintf(stderr, "Measurements not yet supported\n"); + exit(1); +/* + if (circuit->has_intermediate_measurements) { + state = measure(state, op, circuit); + } + else { + double p; + // don't set state = post measurement state + qmdd_measure_all(state, circuit->qreg_size, circuit->creg, &p); // TODO: convert to mtbdd + if (circuit->reversed_qubit_order) { + reverse_bit_array(circuit->creg, circuit->qreg_size); + } + break; + } +*/ + } + if (count_nodes) { +uint64_t count = 0; // TODO: make mtbdd_countnodes(state); + if (count > stats.max_nodes) stats.max_nodes = count; + } + op = op->next; + } + stats.simulation_time = wctime() - t_start; + stats.final_state = state; + stats.shots = 1; +//stats.final_nodes = mtbdd_countnodes(state); +//stats.norm = qmdd_get_norm(state, circuit->qreg_size); // TODO: convert to mtbdd (calculate L2 norm vector), Sebastiaan +} + +/** + * + * Command Line Interface, parse arguments, parse QASM file, start simulation. + * + */ +int main(int argc, char *argv[]) +{ + argp_parse(&argp, argc, argv, 0, 0, 0); + + quantum_circuit_t* circuit = parse_qasm_file(qasm_inputfile); + optimize_qubit_order(circuit, false); + + if (rseed == 0) rseed = time(NULL); + srand(rseed); + + // Standard Lace initialization + lace_start(workers, 0); + + // Simple Sylvan initialization + sylvan_set_sizes(min_tablesize, max_tablesize, min_cachesize, max_cachesize); + sylvan_init_package(); + sylvan_init_mtbdd(); + + simulate_circuit(circuit); + + fprint_stats(stdout, circuit); + if (json_outputfile != NULL) { + FILE *fp = fopen(json_outputfile, "w"); + fprint_stats(fp, circuit); + fclose(fp); + } + + sylvan_quit(); + lace_stop(); + free_quantum_circuit(circuit); + + return 0; +} diff --git a/qasm/sim_qasm.c b/qasm/run_qasm_on_qmdd.c similarity index 95% rename from qasm/sim_qasm.c rename to qasm/run_qasm_on_qmdd.c index 2a5d834..a925504 100644 --- a/qasm/sim_qasm.c +++ b/qasm/run_qasm_on_qmdd.c @@ -1,3 +1,20 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + #include #include #include @@ -5,7 +22,6 @@ #include "qsylvan.h" #include "simple_parser.h" - /*********************************************/ static int workers = 1; @@ -40,6 +56,7 @@ static struct argp_option options[] = {"disable-inv-caching", 12, 0, 0, "Disable storing inverse of MUL and DIV in cache.", 0}, {0, 0, 0, 0, 0, 0} }; + static error_t parse_opt(int key, char *arg, struct argp_state *state) { @@ -160,13 +177,17 @@ wctime() return (tv.tv_sec + 1E-6 * tv.tv_usec); } - +/** + * Here we match the name of a gate in QASM to + * the GATEID + */ QMDD apply_gate(QMDD state, quantum_op_t* gate, BDDVAR nqubits) { // TODO: move this relation between parsed quantum_op and internal gate // somewhere else? stats.applied_gates++; - if (strcmp(gate->name, "id") == 0) { + + if (strcmp(gate->name, "id") == 0) { stats.applied_gates--; return state; } @@ -274,6 +295,7 @@ QMDD apply_gate(QMDD state, quantum_op_t* gate, BDDVAR nqubits) state = qmdd_cgate(state, GATEID_H, gate->ctrls[0], gate->targets[0], nqubits); // CCNOT state = qmdd_cgate2(state, GATEID_X, gate->ctrls[0], gate->targets[0], gate->targets[1], nqubits); + return state; } else if (strcmp(gate->name, "rccx") == 0) { @@ -291,6 +313,7 @@ QMDD apply_gate(QMDD state, quantum_op_t* gate, BDDVAR nqubits) state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1], nqubits); state = qmdd_gate(state, GATEID_Phase(gate->angle[0]), gate->targets[1]); state = qmdd_cgate(state, GATEID_X, gate->targets[0], gate->targets[1], nqubits); + return state; } else if (strcmp(gate->name, "rxx") == 0) { @@ -332,6 +355,7 @@ void simulate_circuit(quantum_circuit_t* circuit) while (op != NULL) { if (op->type == op_gate) { state = apply_gate(state, op, circuit->qreg_size); + } else if (op->type == op_measurement) { if (circuit->has_intermediate_measurements) { diff --git a/qasm/simple_parser.cpp b/qasm/simple_parser.cpp index b20f799..b8af1e3 100644 --- a/qasm/simple_parser.cpp +++ b/qasm/simple_parser.cpp @@ -1,3 +1,20 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + #include #include #include @@ -198,7 +215,7 @@ class QASMParser { case ins_type::classical_cond: parse_error("Classical conditioning currently unsupported"); default: - parse_error("Unsupported instruction type"); + parse_error("Unsupported instruction type (" + line + ")"); } } diff --git a/qasm/simple_parser.h b/qasm/simple_parser.h index 085e991..1387357 100644 --- a/qasm/simple_parser.h +++ b/qasm/simple_parser.h @@ -1,55 +1,91 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + #ifdef __cplusplus extern "C" { #endif + +/** + * The structure quantum_circuit_t contains the QASM quantum circuit as defined in the QASM file. + * + * It is composed of the structure quantum_op_t which is a linked-list of gates or gate-control + * quantum subcircuits of the quantum circuit. + * + * Every element is put on the heap (malloc'ed) and should be released after use with free_quantum_ops(). + * + */ typedef enum operation_type { - op_gate, - op_measurement, - op_blank + + op_gate, // gate element + op_measurement, // measurement element + op_blank // empty element + } operation_type_t; typedef struct quantum_op_s { - operation_type_t type; - char name[32]; - double angle[3]; - int targets[2]; - int ctrls[3]; - int meas_dest; + + operation_type_t type; // Element can be a gate or a measurement or blank + char name[32]; // QASM name of this gate + double angle[3]; // Arguments for Rotation gates(theta) or U(theta, phi, lambda), could be empty = 0,0,0, double -> mpfr, after library parser + int targets[2]; // Index of qubit which is connected to the first or second non-identity gate + int ctrls[3]; // Index of qubit which is connected to the first, second or third control o + int meas_dest; // Measurement dest struct quantum_op_s* next; + } quantum_op_t; typedef struct quantum_circuit_s { - char name[200]; - int qreg_size; - int creg_size; - bool *creg; - bool has_intermediate_measurements; - quantum_op_t *operations; - bool reversed_qubit_order; + + char name[200]; // Circuit name + int qreg_size; // qubit register size + int creg_size; // classical register size + bool *creg; // bool is used as bit here, classical register ( size = sizeof(bool)*creg_size ) + bool has_intermediate_measurements; // If true measurements not only at the end + quantum_op_t *operations; // First element of linked-list with all gates or gate-coontrol pairs + bool reversed_qubit_order; // Order of indices of the qubits, if true numbering is from bottom to top + } quantum_circuit_t; /** - * quantum_op_t is a linked-list where every element is malloc'ed. - * free_quantum_ops() should be called on this after use. + * Parser of the QASM file, returns the quantum circuit in above structures. */ quantum_circuit_t* parse_qasm_file(char *filepath); - /** * Invert qubit order if that yields less controls below target qubits. -*/ + */ void optimize_qubit_order(quantum_circuit_t *circuit, bool allow_swaps); + /** - * free()'s all quantum_ops from and including 'first'. -*/ + * Free all quantum elements found in quantum_op_s including 'first'. What is 'first'? + */ void free_quantum_circuit(quantum_circuit_t* circuit); +/** + * Print operations of quantum circuit, quantum elements or control registers. + */ void print_quantum_op(quantum_op_t *op); void print_quantum_circuit(quantum_circuit_t* circuit); void fprint_creg(FILE *stream, quantum_circuit_t* circuit); + #ifdef __cplusplus } #endif \ No newline at end of file diff --git a/qasm/test/test_sim_qasm.py b/qasm/test/test_sim_qasm.py index 3cf32f8..1802aec 100644 --- a/qasm/test/test_sim_qasm.py +++ b/qasm/test/test_sim_qasm.py @@ -9,7 +9,7 @@ TOLERANCE = 1e-6 -SIM_QASM = './build/qasm/sim_qasm' +SIM_QASM = './build/qasm/run_qasm_on_qmdd' QASM_DIR = 'qasm/circuits/' diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b15fcc8..c54dc46 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,7 +5,9 @@ target_sources(qsylvan PRIVATE sha2.c qsylvan_gates.c + qsylvan_gates_mtbdd_mpc.c qsylvan_simulator.c + qsylvan_simulator_mtbdd.c sylvan_aadd.c sylvan_bdd.c sylvan_cache.c @@ -27,7 +29,9 @@ target_sources(qsylvan sylvan.h qsylvan.h qsylvan_gates.h + qsylvan_gates_mtbdd_mpc.h qsylvan_simulator.h + qsylvan_simulator_mtbdd.h sylvan_aadd.h sylvan_aadd_int.h sylvan_bdd.h diff --git a/src/qsylvan_gates.h b/src/qsylvan_gates.h index f3dc681..f8f0244 100644 --- a/src/qsylvan_gates.h +++ b/src/qsylvan_gates.h @@ -43,20 +43,25 @@ void qmdd_gates_init(); // we would like a (for example) pi/16 gate to always have the same unique ID // throughout the entire run of the circuit. // TODO: maybe just use GATEID_dynamic for this + void qmdd_phase_gates_init(int n); + static inline uint32_t GATEID_Rk(int k) { return k + n_predef_gates; }; + // Another 255 parameterized phase gates, but this time with negative angles. static inline uint32_t GATEID_Rk_dag(int k){ return k + (n_predef_gates+256); }; // Reserve 'num_dynamic_gates' gate IDs (comibned) for of the following: // Rx, Ry, Rz. These gate IDs are re-used when a user requires more than // 'num_dynamic_gates' custom gates. + /** * Rotation around x-axis with angle theta. * NOTE: This GATEID is re-used for parametrized gates. The returned ID only * corresponds to the Rx(theta) gate until the next parametrized gate is created. */ uint32_t GATEID_Rx(fl_t theta); + /** * Rotation around y-axis with angle theta. * NOTE: This GATEID is re-used for parametrized gates. The returned ID only @@ -64,18 +69,22 @@ uint32_t GATEID_Rx(fl_t theta); * a custom gate id. */ uint32_t GATEID_Ry(fl_t theta); + /** * Rotation around z-axis with angle theta. * NOTE: This GATEID is re-used for parametrized gates. The returned ID only * corresponds to the Rz(theta) gate until the next parametrized gate is created. */ + uint32_t GATEID_Rz(fl_t theta); /** * Rotation around z-axis with angle theta (but different global phase than Rz) * NOTE: This GATEID is re-used for parametrized gates. The returned ID only * corresponds to the P(theta) gate until the next parametrized gate is created. */ + uint32_t GATEID_Phase(fl_t theta); + /** * Generic single-qubit rotation gate with 3 Euler angles. * NOTE: This GATEID is re-used for parametrized gates. The returned ID only @@ -84,6 +93,4 @@ uint32_t GATEID_Phase(fl_t theta); */ uint32_t GATEID_U(fl_t theta, fl_t phi, fl_t lambda); - - #endif diff --git a/src/qsylvan_gates_mtbdd_mpc.c b/src/qsylvan_gates_mtbdd_mpc.c new file mode 100644 index 0000000..3ad08f5 --- /dev/null +++ b/src/qsylvan_gates_mtbdd_mpc.c @@ -0,0 +1,717 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +#include +//#include +//#include +#include + +/** + * Definition of Quantum mpc values to fill the fixed and dynamic gates. + * + * They are defined to be used globally to improve computation performance for the dynamic gates. + * + * After gate initialization and use, call gate exit function to clear these mpc values. + * +*/ + +// +// For the fixed gates +// + +// static MTBDD I_dd = MTBDD_ZERO; // moved to header file +// static MTBDD X_dd = MTBDD_ZERO; +// static MTBDD Y_dd = MTBDD_ZERO; +// static MTBDD Z_dd = MTBDD_ZERO; + +// static MTBDD H_dd = MTBDD_ZERO; +// static MTBDD S_dd = MTBDD_ZERO; +// static MTBDD S_dag_dd = MTBDD_ZERO; +// static MTBDD T_dd = MTBDD_ZERO; +// static MTBDD T_dag_dd = MTBDD_ZERO; + +// static MTBDD sqrt_X_dd = MTBDD_ZERO; +// static MTBDD sqrt_X_dag_dd = MTBDD_ZERO; +// static MTBDD sqrt_Y_dd = MTBDD_ZERO; +// static MTBDD sqrt_Y_dag_dd = MTBDD_ZERO; + +static mpc_ptr mpc_pi; +static mpc_ptr mpc_zero; // 0.0 + i 0.0 +static mpc_ptr mpc_re_one; // 1.0 + i 0.0 +static mpc_ptr mpc_re_one_min; // 0.0 - i 1.0 +static mpc_ptr mpc_im_one; // 0.0 + i 1.0 +static mpc_ptr mpc_im_one_min; // 0.0 - i 0.0 +static mpc_ptr mpc_half_half; // 0.5 + i 0.5 +static mpc_ptr mpc_half_half_min; // 0.5 - i 0.5 +static mpc_ptr mpc_half_min_half; // -0.5 + i 0.5 +static mpc_ptr mpc_half_min_half_min; // -0.5 - i 0.5 + +static mpc_ptr mpc_sqrt_2; // sqrt(2) +static mpc_ptr mpc_res_sqrt_2; // 1/sqrt(2) +static mpc_ptr mpc_res_sqrt_2_min; // -1/sqrt(2) +static mpc_ptr mpc_res_sqrt_2_res_sqrt_2; // cos(pi/4) + i sin(pi/4) = 1/sqrt(2) + i/sqrt(2) +static mpc_ptr mpc_res_sqrt_2_res_sqrt_2_min; // cos(pi/4) - i sin(pi/4) = 1/sqrt(2) - i/sqrt(2) + +// +// For the dynamic initialized gates with theta Rx, Ry, Rz +// + +#define MAX_QUBITS 128 + +static MTBDD R_dd[MAX_QUBITS], R_dag_dd[MAX_QUBITS]; + +static mpfr_t mpfr_pi; +static mpfr_t mpfr_zero; +static mpfr_t mpfr_theta; +static mpfr_t mpfr_theta_2; +static mpfr_t mpfr_min_theta_2; +static mpfr_t mpfr_cos_theta_2; +static mpfr_t mpfr_sin_theta_2; + +static mpc_t mpc_cos_theta_2; // cos(theta/2) + i 0.0 +static mpc_t mpc_zero_sin_min_theta_2; // 0.0 - i sin(theta/2) + +static mpc_t mpc_sin_min_theta_2; // -sin(theta/2) + i 0.0 +static mpc_t mpc_sin_theta_2; // sin(theta/2) + i 0.0 + +static mpfr_t mpfr_cos_min_theta_2; +static mpfr_t mpfr_sin_min_theta_2; +static mpc_t mpc_exp_min_theta_2; +static mpc_t mpc_exp_theta_2; + +static mpfr_t mpfr_cos_theta; // cos(theta) +static mpfr_t mpfr_sin_theta; // sin(theta/2) + +static mpc_t mpc_exp_theta; // cos(theta) + i sin(theta) + +static mpfr_t mpfr_phi; +static mpfr_t mpfr_lambda; +static mpfr_t mpfr_gam; + +static mpfr_t mpfr_min_sin_theta_2, mpfr_sin_theta_2; // -sin(theta/2) + +static mpc_t mpc_exp_lambda; +static mpfr_t mpfr_cos_lambda, mpfr_sin_lambda; // cos(lambda) + i sin(lambda) +static mpc_t mpc_exp_lambda_mul_min_sin_theta_2, mpc_exp_lambda; +static mpfr_t mpfr_min_sin_theta_2; + +static mpc_t mpc_exp_phi; +static mpfr_t mpfr_cos_phi, mpfr_sin_phi; // cos(phi) + i sin(phi) +static mpc_t mpc_exp_phi_mul_sin_theta_2, mpc_exp_phi; +static mpfr_t mpfr_sin_theta_2; + +static mpc_t mpc_exp_gam; +static mpfr_t mpfr_cos_gam, mpfr_sin_gam; // cos(gamma) + i sin(gamma) +static mpc_t mpc_exp_gam_mul_cos_theta_2, mpc_exp_gam; +static mpfr_t mpfr_cos_theta_2; + + +// All fixed 2x2 gates +//uint64_t G[nr_predef_gates+256+256][2][2]; // G[gateid][row][col] + +/** + * "Constructor" + */ +void +mtbdd_gates_init_mpc() +{ + mtbdd_gate_init_fixed_variables(); + mtbdd_gate_init_dynamic_variables(); + mtbdd_fixed_gates_init_mpc(); + + return; +} + +/** + * "Destructor" + */ +void +mtbdd_gate_exit_mpc() +{ + mpc_clear(mpc_pi); + + mpc_clear(mpc_pi); + mpc_clear(mpc_zero); + mpc_clear(mpc_re_one); + mpc_clear(mpc_re_one_min); + mpc_clear(mpc_im_one); + mpc_clear(mpc_im_one_min); + mpc_clear(mpc_half_half); + mpc_clear(mpc_half_half_min); + mpc_clear(mpc_half_min_half); + mpc_clear(mpc_half_min_half_min); + + mpc_clear(mpc_sqrt_2); + mpc_clear(mpc_res_sqrt_2); + mpc_clear(mpc_res_sqrt_2_min); + mpc_clear(mpc_res_sqrt_2_res_sqrt_2); + mpc_clear(mpc_res_sqrt_2_res_sqrt_2_min); + + mpfr_clear(mpfr_theta); + mpfr_clear(mpfr_theta_2); + mpfr_clear(mpfr_cos_theta_2); + mpfr_clear(mpfr_sin_theta_2); + mpfr_clear(mpfr_cos_min_theta_2); + mpfr_clear(mpfr_sin_min_theta_2); + + mpc_clear(mpc_exp_min_theta_2); + mpc_clear(mpc_exp_theta_2); + + return; +} + +/** + * Initialize mpc variables to be used globally. + */ +void +mtbdd_gate_init_fixed_variables() +{ + mpfr_init2(mpfr_pi, MPC_PRECISION); + mpfr_const_pi(mpfr_pi, MPC_ROUNDING); + + mpc_init2(mpc_pi, MPC_PRECISION); + mpc_set_fr(mpc_pi, mpfr_pi, MPC_ROUNDING); + + mpc_assign(mpc_zero, 0.0, 0.0); // 0.0 + i 0.0 + mpc_assign(mpc_re_one, 1.0, 0.0); // 1.0 + i 0.0 + mpc_assign(mpc_re_one_min, -1.0, 0.0); // ... + mpc_assign(mpc_im_one, 0.0, 1.0); + mpc_assign(mpc_im_one_min, 0.0, -1.0); + mpc_assign(mpc_half_half, 0.5, 0.5); + mpc_assign(mpc_half_half_min, 0.5, -0.5); + mpc_assign(mpc_half_min_half, -0.5, 0.5); + mpc_assign(mpc_half_min_half_min, -0.5, -0.5); + + mpc_sqrt_assign(mpc_sqrt_2, 2.0, 0.0); // sqrt(2.0) + + mpc_init2(mpc_res_sqrt_2, MPC_PRECISION); + mpc_div(mpc_res_sqrt_2, mpc_re_one, mpc_sqrt_2, MPC_ROUNDING); // 1/sqrt(2.0) + + mpc_init2(mpc_res_sqrt_2, MPC_PRECISION); + mpc_div(mpc_res_sqrt_2_min, mpc_re_one_min, mpc_sqrt_2, MPC_ROUNDING); // -1/sqrt(2.0) + + mpc_init2(mpc_res_sqrt_2_res_sqrt_2, MPC_PRECISION); + mpc_add(mpc_res_sqrt_2_res_sqrt_2, mpc_res_sqrt_2, mpc_res_sqrt_2, MPC_ROUNDING); // 1/sqrt(2.0) + i 1/sqrt(2.0) + + mpc_init2(mpc_res_sqrt_2_res_sqrt_2_min, MPC_PRECISION); + mpc_add(mpc_res_sqrt_2_res_sqrt_2_min, mpc_res_sqrt_2, mpc_res_sqrt_2_min, MPC_ROUNDING); // 1/sqrt(2.0) - i 1/sqrt(2.0) + + return; +} + +void +mtbdd_gate_init_dynamic_variables() +{ + mpfr_init2(mpfr_theta, MPC_PRECISION); + mpfr_init2(mpfr_theta_2, MPC_PRECISION); + mpfr_init2(mpfr_cos_theta_2, MPC_PRECISION); + mpfr_init2(mpfr_sin_theta_2, MPC_PRECISION); + mpfr_init2(mpfr_cos_min_theta_2, MPC_PRECISION); + mpfr_init2(mpfr_sin_min_theta_2, MPC_PRECISION); + + mpc_init2(mpc_exp_min_theta_2, MPC_PRECISION); + mpc_init2(mpc_exp_theta_2, MPC_PRECISION); + + return; +} + +/** + * Initialize all fixed gates with mpc type. + */ +void +mtbdd_fixed_gates_init_mpc() +{ + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_re_one; + + I_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_zero; + G_arr[0][1] = mpc_re_one; + G_arr[1][0] = mpc_re_one; + G_arr[1][1] = mpc_zero; + + X_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_zero; + G_arr[0][1] = mpc_im_one_min; + G_arr[1][0] = mpc_im_one; + G_arr[1][1] = mpc_zero; + + Y_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_re_one_min; + + Z_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_res_sqrt_2; + G_arr[0][1] = mpc_res_sqrt_2; + G_arr[1][0] = mpc_res_sqrt_2; + G_arr[1][1] = mpc_res_sqrt_2_min; + + H_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_im_one; + + S_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_im_one_min; + + S_dag_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_res_sqrt_2_res_sqrt_2; + + T_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_res_sqrt_2_res_sqrt_2_min; + + T_dag_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_half_half; + G_arr[0][1] = mpc_half_half_min; + G_arr[1][0] = mpc_half_half_min; + G_arr[1][1] = mpc_half_half; + + sqrt_X_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_half_half_min; + G_arr[0][1] = mpc_half_half; + G_arr[1][0] = mpc_half_half; + G_arr[1][1] = mpc_half_half_min; + + sqrt_X_dag_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_half_half; + G_arr[0][1] = mpc_half_min_half_min; + G_arr[1][0] = mpc_half_half; + G_arr[1][1] = mpc_half_half; + + sqrt_Y_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_half_half_min; + G_arr[0][1] = mpc_half_half_min; + G_arr[1][0] = mpc_half_min_half; + G_arr[1][1] = mpc_half_half_min; + + sqrt_Y_dag_dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return; +} + +/* + k = GATEID_I; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_re_one; + + k = GATEID_X; + G[k][0][0] = mpc_zero; G[k][0][1] = mpc_re_one; + G[k][1][0] = mpc_re_one; G[k][1][1] = mpc_zero; + + k = GATEID_Y; + G[k][0][0] = mpc_zero; G[k][0][1] = mpc_im_one_min; + G[k][1][0] = mpc_im_one; G[k][1][1] = mpc_zero; + + k = GATEID_Z; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_re_one_min; + + k = GATEID_H; + G[k][0][0] = mpc_res_sqrt_2; G[k][0][1] = mpc_res_sqrt_2; + G[k][1][0] = mpc_res_sqrt_2; G[k][1][1] = mpc_res_sqrt_2_min; + + k = GATEID_S; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_im_one; + + k = GATEID_Sdag; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_im_one_min; + + k = GATEID_T; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_res_sqrt_2_res_sqrt_2; + + k = GATEID_Tdag; + G[k][0][0] = mpc_re_one; G[k][0][1] = mpc_zero; + G[k][1][0] = mpc_zero; G[k][1][1] = mpc_res_sqrt_2_res_sqrt_2_min; + + k = GATEID_sqrtX; + G[k][0][0] = mpc_half_half; G[k][0][1] = mpc_half_half_min; + G[k][1][0] = mpc_half_half_min; G[k][1][1] = mpc_half_half; + + k = GATEID_sqrtXdag; + G[k][0][0] = mpc_half_half_min; G[k][0][1] = mpc_half_half; + G[k][1][0] = mpc_half_half; G[k][1][1] = mpc_half_half_min; + + k = GATEID_sqrtY; + G[k][0][0] = mpc_half_half; G[k][0][1] = mpc_half_min_half_min; + G[k][1][0] = mpc_half_half; G[k][1][1] = mpc_half_half; + + k = GATEID_sqrtYdag; + G[k][0][0] = mpc_half_half_min; G[k][0][1] = mpc_half_half_min; + G[k][1][0] = mpc_half_min_half; G[k][1][1] = mpc_half_half_min; +*/ + +/* + qmdd_phase_gates_init(255); + + // init dynamic gate + // (necessary when qmdd_gates_init() is called after gc to re-init all gates) + k = GATEID_dynamic; + gates[k][0] = weight_lookup(&dynamic_gate[0]); + gates[k][1] = weight_lookup(&dynamic_gate[1]); + gates[k][2] = weight_lookup(&dynamic_gate[2]); + gates[k][3] = weight_lookup(&dynamic_gate[3]); + + return; +} +*/ + +/** + * Parameterized rotation gates Rx(theta), Ry(theta), Rz(theta) + */ + +MTBDD +mtbdd_Rx(double theta) // TODO: change in mpfr! +{ + // Calculate complex numbers for gate elements + + mpfr_set_d(mpfr_zero, 0.0, MPC_ROUNDING); + mpfr_set_d(mpfr_theta, theta, MPC_ROUNDING); + mpfr_div_ui(mpfr_theta_2, mpfr_theta, 2, MPC_ROUNDING); // theta / 2 + mpfr_neg(mpfr_min_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_cos(mpfr_cos_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // cos(-theta/2) + mpfr_sin(mpfr_sin_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // sin(-theta/2) + + mpc_set_fr_fr(mpc_cos_theta_2, mpfr_cos_theta_2, mpfr_zero, MPC_ROUNDING); // cos(theta/2) + i 0.0 + mpc_set_fr_fr(mpc_zero_sin_min_theta_2, mpfr_zero, mpfr_sin_min_theta_2, MPC_ROUNDING); // 0.0 - i sin(theta/2) + + // Convert to MTBDD + + MTBDD dd = MTBDD_ZERO; + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_cos_theta_2; + G_arr[0][1] = mpc_zero_sin_min_theta_2; + G_arr[1][0] = mpc_zero_sin_min_theta_2; + G_arr[1][1] = mpc_cos_theta_2; + + dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return dd; +} + +MTBDD +mtbdd_Ry(double theta) // TODO: change in mpfr! +{ + // Calculate complex numbers for gate elements + + mpfr_set_d(mpfr_theta, theta, MPC_ROUNDING); + mpfr_div_ui(mpfr_theta_2, mpfr_theta, 2, MPC_ROUNDING); // theta / 2 + mpfr_neg(mpfr_min_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_cos(mpfr_cos_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // cos(-theta/2) + mpfr_sin(mpfr_sin_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // sin(-theta/2) + + mpc_set_fr_fr(mpc_cos_theta_2, mpfr_cos_theta_2, mpfr_zero, MPC_ROUNDING); // cos(theta/2) + i 0.0 + mpc_set_fr_fr(mpc_sin_min_theta_2, mpfr_sin_min_theta_2, mpfr_zero, MPC_ROUNDING); // -sin(theta/2) + i 0.0 + mpc_set_fr_fr(mpc_sin_theta_2, mpfr_sin_theta_2, mpfr_zero, MPC_ROUNDING); // sin(theta/2) + i 0.0 + + // Convert to MTBDD + + MTBDD dd = MTBDD_ZERO; + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_cos_theta_2; + G_arr[0][1] = mpc_sin_min_theta_2; + G_arr[1][0] = mpc_sin_theta_2; + G_arr[1][1] = mpc_cos_theta_2; + + dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return dd; +} + +MTBDD +mtbdd_Rz(double theta) // TODO: change in mpfr! +{ + // Calculate complex numbers for gate elements + + mpfr_set_d(mpfr_theta, theta, MPC_ROUNDING); + mpfr_div_ui(mpfr_theta_2, mpfr_theta, 2, MPC_ROUNDING); // theta / 2 + mpfr_neg(mpfr_min_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_cos(mpfr_cos_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // cos(-theta/2) + mpfr_sin(mpfr_sin_min_theta_2, mpfr_min_theta_2, MPC_ROUNDING); // sin(-theta/2) + + mpc_set_fr_fr(mpc_exp_theta_2, mpfr_cos_theta_2, mpfr_sin_theta_2, MPC_ROUNDING); // cos(theta/2) + i sin(theta/2) + mpc_set_fr_fr(mpc_exp_min_theta_2, mpfr_cos_min_theta_2, mpfr_sin_min_theta_2, MPC_ROUNDING); // cos(-theta/2) + i sin(-theta/2) + + // Convert to MTBDD + + MTBDD dd = MTBDD_ZERO; + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_exp_min_theta_2; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_exp_theta_2; + + dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return dd; +} + +MTBDD +mtbdd_Phase(double theta) // TODO: change in mpfr! +{ + // Calculate complex numbers for gate elements + + mpfr_set_d(mpfr_theta, theta, MPC_ROUNDING); + + mpfr_cos(mpfr_cos_theta, mpfr_theta, MPC_ROUNDING); // cos(theta) + mpfr_sin(mpfr_sin_theta, mpfr_theta, MPC_ROUNDING); // sin(theta) + + mpc_set_fr_fr(mpc_exp_theta, mpfr_cos_theta, mpfr_sin_theta, MPC_ROUNDING); // cos(theta) + i sin(theta) + + // Convert to MTBDD + + MTBDD dd = MTBDD_ZERO; + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_re_one; // 1.0 + G_arr[0][1] = mpc_zero; // 0.0 + G_arr[1][0] = mpc_zero; // 0.0 + G_arr[1][1] = mpc_exp_theta; // exp(i theta) = cos(theta) + i sin(theta) + + dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return dd; +} + +MTBDD +mtbdd_U(double theta, double phi, double lambda) +{ + // Calculate complex numbers for gate elements + + mpfr_set_d(mpfr_theta, theta, MPC_ROUNDING); + mpfr_set_d(mpfr_phi, phi, MPC_ROUNDING); + 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_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) + + + 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_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); // (cos(gamma) + i sin(gamma)) x cos(theta/2) + + // Convert to MTBDD + + MTBDD dd = MTBDD_ZERO; + int n = 1; + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_cos_theta_2; // cos(theta/2) + G_arr[0][1] = mpc_exp_lambda_mul_min_sin_theta_2; // (cos(lambda) + i sin(lambda)) x -sin(theta/2) + G_arr[1][0] = mpc_exp_phi_mul_sin_theta_2; // (cos(phi) + i sin(phi)) x sin(theta/2) + G_arr[1][1] = mpc_exp_gam_mul_cos_theta_2; // (cos(phi+lambda) + i sin(phi+lambda)) x cos(theta/2) + + dd = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + return dd; +} + +/* +void +qmdd_phase_gates_init(int n) // n is number of qubits? +{ + // add gate R_k to gates table + // (note that R_0 = I, R_1 = Z, R_2 = S, R_4 = T) + uint32_t gate_id; + fl_t angle; + complex_t cartesian; +*/ + +void +mtbdd_phase_gates_init(int n) +{ + mpfr_t mpfr_2_pi; + mpfr_t mpfr_2; + + mpfr_t mpfr_2_power_k; + mpfr_t mpfr_2_pi_div_2_power_k; + mpfr_t mpfr_2_pi_div_2_power_k_min; + + mpfr_t mpfr_cos_theta_min; + mpfr_t mpfr_sin_theta_min; + + mpc_t mpc_exp_theta_min; + + mpfr_init2(mpfr_2_pi, MPC_PRECISION); + mpfr_init2(mpfr_2, MPC_PRECISION); + mpfr_init2(mpfr_2_power_k, MPC_PRECISION); + mpfr_init2(mpfr_2_pi_div_2_power_k, MPC_PRECISION); + + mpfr_mul_si(mpfr_2_pi, mpfr_pi, 2, MPC_ROUNDING); // 2 * pi + mpfr_set_si(mpfr_2, 2, MPC_ROUNDING); + + // MTBDD R_dd[n], R_dag_dd[n]; + + for (int k=0; k<=n; k++) { + + mpfr_pow_si(mpfr_2_power_k, mpfr_2, (long)k, MPC_ROUNDING); // 2 ^ k + mpfr_div(mpfr_2_pi_div_2_power_k, mpfr_2_pi, mpfr_2_power_k, MPC_ROUNDING); // 2 * pi / (2 ^ k) + mpfr_neg(mpfr_2_pi_div_2_power_k_min, mpfr_2_pi_div_2_power_k, MPC_ROUNDING); // -2 * pi / (2 ^ k) + + mpfr_cos(mpfr_cos_theta, mpfr_2_pi_div_2_power_k, MPC_ROUNDING); // cos(theta) + mpfr_sin(mpfr_sin_theta, mpfr_2_pi_div_2_power_k, MPC_ROUNDING); // sin(theta) + + mpc_set_fr_fr(mpc_exp_theta, mpfr_cos_theta, mpfr_sin_theta, MPC_ROUNDING); // cos(theta) + i sin(theta) + + mpfr_cos(mpfr_cos_theta_min, mpfr_2_pi_div_2_power_k_min, MPC_ROUNDING); // cos(-theta) + mpfr_sin(mpfr_sin_theta_min, mpfr_2_pi_div_2_power_k_min, MPC_ROUNDING); // sin(-theta) + + mpc_set_fr_fr(mpc_exp_theta_min, mpfr_cos_theta_min, mpfr_sin_theta_min, MPC_ROUNDING); // cos(-theta) + i sin(-theta) + + R_dd[k] = MTBDD_ZERO; + R_dag_dd[k] = MTBDD_ZERO; + + mpc_ptr **G_arr = NULL; + allocate_matrix_array_mpc(&G_arr, n); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_exp_theta; + + R_dd[k] = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + G_arr[0][0] = mpc_re_one; + G_arr[0][1] = mpc_zero; + G_arr[1][0] = mpc_zero; + G_arr[1][1] = mpc_exp_theta_min; + + R_dag_dd[k] = matrix_array_to_mtbdd_mpc(G_arr, n, ALTERNATE_ROW_FIRST_WISE_MODE); + + free_matrix_array_mpc(G_arr, n); + + } + + mpfr_clear(mpfr_2_pi); + mpfr_clear(mpfr_2); + + mpfr_clear(mpfr_2_power_k); + mpfr_clear(mpfr_2_pi_div_2_power_k); + mpfr_clear(mpfr_2_pi_div_2_power_k_min); + + mpfr_clear(mpfr_cos_theta_min); + mpfr_clear(mpfr_sin_theta_min); + + mpc_clear(mpc_exp_theta_min); + +/* + for (int k=0; k<=n; k++) { + + // forward rotation + angle = 2 * Pi / (1< +#include + +/** + * + * Global defined mtbdd's of fixed gates. + * + */ + +static MTBDD I_dd = MTBDD_ZERO; +static MTBDD X_dd = MTBDD_ZERO; +static MTBDD Y_dd = MTBDD_ZERO; +static MTBDD Z_dd = MTBDD_ZERO; + +static MTBDD H_dd = MTBDD_ZERO; +static MTBDD S_dd = MTBDD_ZERO; +static MTBDD S_dag_dd = MTBDD_ZERO; +static MTBDD T_dd = MTBDD_ZERO; +static MTBDD T_dag_dd = MTBDD_ZERO; + +static MTBDD sqrt_X_dd = MTBDD_ZERO; +static MTBDD sqrt_X_dag_dd = MTBDD_ZERO; +static MTBDD sqrt_Y_dd = MTBDD_ZERO; +static MTBDD sqrt_Y_dag_dd = MTBDD_ZERO; + +/** + * + * Functions + * + */ + +// "Constructor" +void +mtbdd_gates_init_mpc(); + +// "Destructor" +void +mtbdd_gate_exit_mpc(); + +// The reason why these are initialized beforehand instead of on-demand is that +// we would like a (for example) pi/16 gate to always have the same unique ID +// throughout the entire run of the circuit. + +void +mtbdd_gate_init_fixed_variables(); + +void +mtbdd_gate_init_dynamic_variables(); + +void +mtbdd_fixed_gates_init_mpc(); + +void +mtbdd_phase_gates_init(int n); + +/** + * Rotation around x-axis with angle theta. + */ +MTBDD +mtbdd_Rx(double theta); + +/** + * Rotation around y-axis with angle theta. + */ +MTBDD +mtbdd_Ry(double theta); + +/** + * Rotation around z-axis with angle theta. + */ +MTBDD +mtbdd_Rz(double theta); + +/** + * Rotation around z-axis with angle theta (but different global phase than Rz) + */ +MTBDD +mtbdd_Phase(double theta); + +/** + * Generic single-qubit rotation gate with 3 Euler angles. + */ +MTBDD +mtbdd_U(double theta, double phi, double lambda); + +#endif diff --git a/src/qsylvan_mtbdd_gates.c b/src/qsylvan_mtbdd_gates.c deleted file mode 100644 index f1dcc0c..0000000 --- a/src/qsylvan_mtbdd_gates.c +++ /dev/null @@ -1,98 +0,0 @@ -/* - * Copyright 2023 System Verification Lab, LIACS, Leiden University - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - * - */ - -#include - -MTBDD gate_to_mtbdd(quantum_gate_t g) -{ - if(g == I) { - mpc_t G[2][2]; - G[0][0] = MTBDD_ONE; G[0][1] = MTBDD_ZERO; - G[1][0] = MTBDD_ZERO; G[1][1] = MTBDD_ONE; - return array_matrix_to_mtbdd(G,2); - } - - // ... - - return MTBDD_ZERO; -} - -/* - Pi = 2.0 * flt_acos(0.0); - - // initialize 2x2 gates (complex values from gates currently stored in - // same table as complex amplitude values) - uint32_t k; - - k = GATEID_I; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = AADD_ONE; - - k = GATEID_X; - gates[k][0] = AADD_ZERO; gates[k][1] = AADD_ONE; - gates[k][2] = AADD_ONE; gates[k][3] = AADD_ZERO; - - k = GATEID_Y; - gates[k][0] = AADD_ZERO; gates[k][1] = weight_lookup(cmake(0.0, -1.0)); - gates[k][2] = weight_lookup(cmake(0.0, 1.0)); gates[k][3] = AADD_ZERO; - - k = GATEID_Z; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = AADD_MIN_ONE; - - k = GATEID_H; - gates[k][0] = gates[k][1] = gates[k][2] = weight_lookup(cmake(1.0/flt_sqrt(2.0),0)); - gates[k][3] = weight_lookup(cmake(-1.0/flt_sqrt(2.0),0)); - - k = GATEID_S; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = weight_lookup(cmake(0.0, 1.0)); - - k = GATEID_Sdag; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = weight_lookup(cmake(0.0, -1.0)); - - k = GATEID_T; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = weight_lookup(cmake(1.0/flt_sqrt(2.0), 1.0/flt_sqrt(2.0))); - - k = GATEID_Tdag; - gates[k][0] = AADD_ONE; gates[k][1] = AADD_ZERO; - gates[k][2] = AADD_ZERO; gates[k][3] = weight_lookup(cmake(1.0/flt_sqrt(2.0), -1.0/flt_sqrt(2.0))); - - k = GATEID_sqrtX; - gates[k][0] = weight_lookup(cmake(0.5, 0.5)); gates[k][1] = weight_lookup(cmake(0.5,-0.5)); - gates[k][2] = weight_lookup(cmake(0.5,-0.5)); gates[k][3] = weight_lookup(cmake(0.5, 0.5)); - - k = GATEID_sqrtXdag; - gates[k][0] = weight_lookup(cmake(0.5,-0.5)); gates[k][1] = weight_lookup(cmake(0.5, 0.5)); - gates[k][2] = weight_lookup(cmake(0.5, 0.5)); gates[k][3] = weight_lookup(cmake(0.5,-0.5)); - - k = GATEID_sqrtY; - gates[k][0] = weight_lookup(cmake(0.5, 0.5)); gates[k][1] = weight_lookup(cmake(-0.5,-0.5)); - gates[k][2] = weight_lookup(cmake(0.5, 0.5)); gates[k][3] = weight_lookup(cmake(0.5, 0.5)); - - k = GATEID_sqrtYdag; - gates[k][0] = weight_lookup(cmake(0.5,-0.5)); gates[k][1] = weight_lookup(cmake(0.5,-0.5)); - gates[k][2] = weight_lookup(cmake(-0.5,0.5)); gates[k][3] = weight_lookup(cmake(0.5,-0.5)); - - qmdd_phase_gates_init(255); - - next_custom_id = 0; - -*/ - diff --git a/src/qsylvan_mtbdd_gates.h b/src/qsylvan_mtbdd_gates.h deleted file mode 100644 index 288395b..0000000 --- a/src/qsylvan_mtbdd_gates.h +++ /dev/null @@ -1,47 +0,0 @@ -/* - * Copyright 2023 System Verification Lab, LIACS, Leiden University - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - * - */ - -#ifndef QSYLVAN_MTBDD_GATES_H -#define QSYLVAN_MTBDD_GATES_H - -#include - -/** - * - * Definition of the gate matrices 2 x 2, with mpc floatingpoint type. - * - */ - -typedef enum quantum_gate { - I, - X, - Y, - Z, - H, - S, - Sdag, - T, - Tdag, - sqrtX, - sqrtXdag, - sqrtY, - sqrtYdag -} quantum_gate_t; - -MTBDD quantum_gate_to_mtbdd(quantum_gate_t g); - -#endif diff --git a/src/qsylvan_mtbdd_simulator.c b/src/qsylvan_mtbdd_simulator.c deleted file mode 100644 index 6e3789d..0000000 --- a/src/qsylvan_mtbdd_simulator.c +++ /dev/null @@ -1,1219 +0,0 @@ -/* - * Copyright 2023 System Verification Lab, LIACS, Leiden University - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - * - */ - -// -// This code is copied in September 2023 from qsylvan_simulator.h/c made by Sebastian Brand -// - -#include - -/* -#include - -static bool testing_mode = 0; // turns on/off (expensive) sanity checks -static int granularity = 1; // operation cache access granularity - - -// *************** ************** - -// Pack 2 BDDVARs (assumed max 8 bits each) -static inline uint32_t -QMDD_PARAM_PACK_16(BDDVAR a, BDDVAR b) -{ - return b<<8 | a; -} - -// Pack 24 bit gateid, 2 possible qubit parameters (e.g. control/target) -static inline uint64_t -GATE_OPID_40(uint32_t gateid, BDDVAR a, BDDVAR b) -{ - uint64_t res = ((uint64_t)b)<<32 | ((uint64_t)a)<<24 | gateid; - return res; -} - -// Pack 24 bit gateid w/ 5 possible qubit parameters (e.g. target/control/range) -static inline uint64_t -GATE_OPID_64(uint32_t gateid, BDDVAR a, BDDVAR b, BDDVAR c, BDDVAR d, BDDVAR e) -{ - uint64_t res = ((uint64_t)e)<<56 | - ((uint64_t)d)<<48 | - ((uint64_t)c)<<40 | - ((uint64_t)b)<<32 | - ((uint64_t)a)<<24 | - gateid; - return res; -} - -// **************************** - -*/ - - - -// ************************************************************ - -/* -void -qsylvan_init_simulator(size_t wgt_tab_size, double wgt_tab_tolerance, int edge_weigth_backend, int norm_strat) -{ - sylvan_init_mtbdd(wgt_tab_size, wgt_tab_tolerance, edge_weigth_backend, norm_strat, &qmdd_gates_init); -} - -void -qsylvan_init_defaults(size_t wgt_tab_size) -{ - qsylvan_init_simulator(wgt_tab_size, -1, COMP_HASHMAP, NORM_LOW); -} -*/ - -// *********************************************************** - - -// **************************************************** - -MTBDD -mtbdd_create_all_zero_state(int n) -{ - if(n == 0) - return MTBDD_ZERO; - - bool qubit[n]; - for (int k=0; k= 0; k--) { - if (qubit[k] == 0) { - low = mtbdd_makenode(k, prev, AADD_ONE); - high = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - } - if (qubit[k] == 1) { - low = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - high = aadd_bundle(AADD_TARGET(prev), AADD_ONE); - } - // add node to unique table - prev = aadd_makenode(k, low, high); - } - - return prev; -*/ -} - -MTBDD -mtbdd_stack_matrix(MTBDD below, BDDVAR k, gate_t gateid) -{ - -/* - // This function effectively does a Kronecker product gate \tensor below - BDDVAR s, t; - QMDD u00, u01, u10, u11, low, high, res; - - // Even + uneven variable are used to encode the 4 values - s = 2*k; - t = s + 1; - - // Matrix U = [u00 u01 - // u10 u11] endoded in a small tree - u00 = aadd_bundle(AADD_TARGET(below), gates[gateid][0]); - u10 = aadd_bundle(AADD_TARGET(below), gates[gateid][2]); - u01 = aadd_bundle(AADD_TARGET(below), gates[gateid][1]); - u11 = aadd_bundle(AADD_TARGET(below), gates[gateid][3]); - low = aadd_makenode(t, u00, u10); - high = aadd_makenode(t, u01, u11); - res = aadd_makenode(s, low, high); - - // Propagate common factor on previous root amp to new root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(below), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; -*/ -} - -QMDD -qmdd_stack_control(QMDD case0, QMDD case1, BDDVAR k) -{ - // Effectively does |0><0| \tensor case0 + |1><1| \tensor case1 - BDDVAR s, t; - QMDD u00, u01, u10, u11, low, high, res; - - s = 2*k; - t = s + 1; - - u00 = case0; - u10 = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - u01 = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - u11 = case1; - low = aadd_makenode(t, u00, u10); - high = aadd_makenode(t, u01, u11); - res = aadd_makenode(s, low, high); - - // Weights of case0/case1 already dealt with by aadd_makenode - return res; -} - -QMDD -qmdd_create_all_identity_matrix(BDDVAR n) -{ - // Start at terminal and build backwards - QMDD prev = aadd_bundle(AADD_TERMINAL, AADD_ONE); - for (int k = n-1; k >= 0; k--) { - prev = qmdd_stack_matrix(prev, k, GATEID_I); - } - return prev; -} - -QMDD -qmdd_create_single_qubit_gate(BDDVAR n, BDDVAR t, gate_id_t gateid) -{ - // Start at terminal and build backwards - QMDD prev = aadd_bundle(AADD_TERMINAL, AADD_ONE); - for (int k = n-1; k >= 0; k--) { - if ((unsigned int)k == t) - prev = qmdd_stack_matrix(prev, k, gateid); - else - prev = qmdd_stack_matrix(prev, k, GATEID_I); - } - return prev; -} - -QMDD -qmdd_create_single_qubit_gates(BDDVAR n, gate_id_t *gateids) -{ - // Start at terminal and build backwards - QMDD prev = aadd_bundle(AADD_TERMINAL, AADD_ONE); - for (int k = n-1; k >= 0; k--) { - prev = qmdd_stack_matrix(prev, k, gateids[k]); - } - return prev; -} - -QMDD -qmdd_create_single_qubit_gates_same(BDDVAR n, gate_id_t gateid) -{ - // Start at terminal and build backwards - QMDD prev = aadd_bundle(AADD_TERMINAL, AADD_ONE); - for (int k = n-1; k >= 0; k--) { - prev = qmdd_stack_matrix(prev, k, gateid); - } - return prev; -} - -QMDD -qmdd_create_cgate(BDDVAR n, BDDVAR c, BDDVAR t, gate_id_t gateid) -{ - // for now, assume t > c - assert(t > c); - // Start at terminal and build backwards - QMDD prev = aadd_bundle(AADD_TERMINAL, AADD_ONE); - QMDD branch0 = AADD_TERMINAL, branch1 = AADD_TERMINAL; - for (int k = n-1; k>= 0; k--) { - if ((unsigned int)k > t || (unsigned int) k < c) { - prev = qmdd_stack_matrix(prev, k, GATEID_I); - } - else if ((unsigned int) k == t) { - branch0 = qmdd_stack_matrix(prev, k, GATEID_I); - branch1 = qmdd_stack_matrix(prev, k, gateid); - } - else if ((unsigned int) k < t && (unsigned int) k > c) { - branch0 = qmdd_stack_matrix(branch0, k, GATEID_I); - branch1 = qmdd_stack_matrix(branch1, k, GATEID_I); - } - else if ((unsigned int) k == c) { - prev = qmdd_stack_control(branch0, branch1, k); - } - else { - assert("all cases should have been covered" && false); - } - } - return prev; -} - -QMDD -qmdd_create_multi_cgate_rec(BDDVAR n, int *c_options, gate_id_t gateid, BDDVAR k) -{ - // (assumes controls above target) - // c_options[k] = -1 -> ignore qubit k (apply I) - // c_options[k] = 0 -> control on q_k = |0> - // c_options[k] = 1 -> control on q_k = |1> - // c_options[k] = 2 -> target qubit (for now assume 1 target) - - // Terminal case - if (k == n) { - return aadd_bundle(AADD_TERMINAL, AADD_ONE); - } - - // TODO: catching GATEID_I avoids exp number of recrusive calls, but this - // entire function might be handled in a better way(?) - if (gateid == GATEID_I) { - QMDD identities = aadd_bundle(AADD_TERMINAL, AADD_ONE); - for (int j = n-1; j >= (int) k; j--) { - identities = qmdd_stack_matrix(identities, j, GATEID_I); - } - return identities; - } - - // Recursively build matrix - BDDVAR next_k = k + 1; - - // -1 : Ignore qubit (apply I) - if (c_options[k] == -1) { - QMDD below = qmdd_create_multi_cgate_rec(n, c_options, gateid, next_k); - return qmdd_stack_matrix(below, k, GATEID_I); - } - // 0 : control on q_k = |0> (apply gateid to low branch) - else if (c_options[k] == 0) { - QMDD case0 = qmdd_create_multi_cgate_rec(n, c_options, gateid, next_k); - QMDD case1 = qmdd_create_multi_cgate_rec(n, c_options, GATEID_I, next_k); - return qmdd_stack_control(case0, case1, k); - } - // 1 : control on q_k = |1> (apply gateid to high branch) - else if (c_options[k] == 1) { - QMDD case0 = qmdd_create_multi_cgate_rec(n, c_options, GATEID_I, next_k); - QMDD case1 = qmdd_create_multi_cgate_rec(n, c_options, gateid, next_k); - return qmdd_stack_control(case0, case1, k); - } - // 2 : target qubit - else if (c_options[k] == 2) { - QMDD below = qmdd_create_multi_cgate_rec(n, c_options, gateid, next_k); - return qmdd_stack_matrix(below, k, gateid); - } - else { - printf("Invalid option %d for qubit %d (options = {-1,0,1,2}\n", c_options[k], k); - exit(1); - } -} - -QMDD -qmdd_create_all_control_phase(BDDVAR n, bool *x) -{ - QMDD identity = aadd_bundle(AADD_TERMINAL, AADD_ONE); - QMDD ccphase = aadd_bundle(AADD_TERMINAL, AADD_ONE); - - // Start with (-1)Z gate on last qubit. Z if control on 1 and -Z if 0. - if (x[n-1] == 1) { - ccphase = qmdd_stack_matrix(ccphase, n-1, GATEID_Z); - } - else if (x[n-1] == 0) { - ccphase = qmdd_stack_matrix(ccphase, n-1, GATEID_Z); - ccphase = aadd_bundle(AADD_TARGET(ccphase), wgt_neg(AADD_WEIGHT(ccphase))); - } - - // Stack remaining controls - for (int k = n-2; k >= 0; k--) { - // "Identity stack" for doing nothing on each qubit's non-control branch - identity = qmdd_stack_matrix(identity, k+1, GATEID_I); - - // Check if this qubit should be controlled on 0 or 1 - if (x[k] == 1) - ccphase = qmdd_stack_control(identity, ccphase, k); - else if (x[k] == 0) - ccphase = qmdd_stack_control(ccphase, identity, k); - } - - return ccphase; -} - -/*****************************************************/ - - - - - -/**************************************************************/ - -static int periodic_gc_nodetable = 0; // trigger for gc of node table -static uint64_t gate_counter = 0; - -static void -qmdd_do_before_gate(QMDD* qmdd) -{ - // check if ctable needs gc - if (aadd_test_gc_wgt_table()) { - aadd_protect(qmdd); - aadd_gc_wgt_table(); - aadd_unprotect(qmdd); - } - - if (periodic_gc_nodetable) { - gate_counter++; - if (gate_counter % periodic_gc_nodetable == 0) { - aadd_protect(qmdd); - sylvan_gc(); - aadd_unprotect(qmdd); - } - } - - // log stuff (if logging is enabled) - qmdd_stats_log(*qmdd); -} - -/* Wrapper for applying a single qubit gate. */ -TASK_IMPL_3(QMDD, qmdd_gate, QMDD, qmdd, gate_id_t, gate, BDDVAR, target) -{ - qmdd_do_before_gate(&qmdd); - return qmdd_gate_rec(qmdd, gate, target); -} - -/* Wrapper for applying a controlled gate with 1 control qubit. */ -TASK_IMPL_4(QMDD, qmdd_cgate, QMDD, qmdd, gate_id_t, gate, BDDVAR, c, BDDVAR, t) -{ - qmdd_do_before_gate(&qmdd); - BDDVAR cs[4] = {c, AADD_INVALID_VAR, AADD_INVALID_VAR, AADD_INVALID_VAR}; - return qmdd_cgate_rec(qmdd, gate, cs, t); -} - -/* Wrapper for applying a controlled gate with 2 control qubits. */ -TASK_IMPL_5(QMDD, qmdd_cgate2, QMDD, qmdd, gate_id_t, gate, BDDVAR, c1, BDDVAR, c2, BDDVAR, t) -{ - qmdd_do_before_gate(&qmdd); - BDDVAR cs[4] = {c1, c2, AADD_INVALID_VAR, AADD_INVALID_VAR}; - return qmdd_cgate_rec(qmdd, gate, cs, t); -} - -/* Wrapper for applying a controlled gate with 3 control qubits. */ -TASK_IMPL_6(QMDD, qmdd_cgate3, QMDD, qmdd, gate_id_t, gate, BDDVAR, c1, BDDVAR, c2, BDDVAR, c3, BDDVAR, t) -{ - qmdd_do_before_gate(&qmdd); - BDDVAR cs[4] = {c1, c2, c3, AADD_INVALID_VAR}; // last pos is a buffer - return qmdd_cgate_rec(qmdd, gate, cs, t); -} - -/* Wrapper for applying a controlled gate where the controls are a range. */ -TASK_IMPL_5(QMDD, qmdd_cgate_range, QMDD, qmdd, gate_id_t, gate, BDDVAR, c_first, BDDVAR, c_last, BDDVAR, t) -{ - qmdd_do_before_gate(&qmdd); - return qmdd_cgate_range_rec(qmdd,gate,c_first,c_last,t); -} - -TASK_IMPL_3(QMDD, qmdd_gate_rec, QMDD, q, gate_id_t, gate, BDDVAR, target) -{ - // Trivial cases - if (AADD_WEIGHT(q) == AADD_ZERO) return q; - - BDDVAR var; - QMDD res, low, high; - aadd_get_topvar(q, target, &var, &low, &high); - assert(var <= target); - - // Check cache - bool cachenow = ((var % granularity) == 0); - if (cachenow) { - if (cache_get3(CACHE_QMDD_GATE, sylvan_false, AADD_TARGET(q), GATE_OPID_40(gate, target, 0), &res)) { - sylvan_stats_count(QMDD_GATE_CACHED); - // Multiply root of res with root of input qmdd - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; - } - } - - if (var == target) { - AMP a_u00 = wgt_mul(AADD_WEIGHT(low), gates[gate][0]); - AMP a_u10 = wgt_mul(AADD_WEIGHT(low), gates[gate][2]); - AMP b_u01 = wgt_mul(AADD_WEIGHT(high), gates[gate][1]); - AMP b_u11 = wgt_mul(AADD_WEIGHT(high), gates[gate][3]); - QMDD low1, low2, high1, high2; - low1 = aadd_bundle(AADD_TARGET(low), a_u00); - low2 = aadd_bundle(AADD_TARGET(high),b_u01); - high1 = aadd_bundle(AADD_TARGET(low), a_u10); - high2 = aadd_bundle(AADD_TARGET(high),b_u11); - aadd_refs_spawn(SPAWN(aadd_plus, high1, high2)); - low = CALL(aadd_plus, low1, low2); - aadd_refs_push(low); - high = aadd_refs_sync(SYNC(aadd_plus)); - aadd_refs_pop(1); - res = aadd_makenode(target, low, high); - } - else { // var < target: not at target qubit yet, recursive calls down - aadd_refs_spawn(SPAWN(qmdd_gate_rec, high, gate, target)); - low = CALL(qmdd_gate_rec, low, gate, target); - aadd_refs_push(low); - high = aadd_refs_sync(SYNC(qmdd_gate_rec)); - aadd_refs_pop(1); - res = aadd_makenode(var, low, high); - } - - // Store not yet "root normalized" result in cache - if (cachenow) { - if (cache_put3(CACHE_QMDD_GATE, sylvan_false, AADD_TARGET(q), GATE_OPID_40(gate, target, 0), res)) - sylvan_stats_count(QMDD_GATE_CACHEDPUT); - } - // Multiply amp res with amp of input qmdd - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; -} - -TASK_IMPL_5(QMDD, qmdd_cgate_rec, QMDD, q, gate_id_t, gate, BDDVAR*, cs, uint32_t, ci, BDDVAR, t) -{ - // Get current control qubit. If no more control qubits, apply gate here - BDDVAR c = cs[ci]; - if (c == AADD_INVALID_VAR || ci > MAX_CONTROLS) { - return CALL(qmdd_gate_rec, q, gate, t); - } - - assert(c < t && "ctrl < target required"); - if (ci > 0) - assert(cs[ci-1] < cs[ci] && "order required for multiple controls"); - - BDDVAR var; - QMDD res, low, high; - aadd_get_topvar(q, c, &var, &low, &high); - assert(var <= c); - - // Check cache - bool cachenow = ((var % granularity) == 0); - if (cachenow) { - if (cache_get3(CACHE_QMDD_CGATE, sylvan_false, AADD_TARGET(q), - GATE_OPID_64(gate, ci, cs[0], cs[1], cs[2], t), - &res)) { - sylvan_stats_count(QMDD_CGATE_CACHED); - // Multiply root amp of res with input root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; - } - } - - // If current node is (one of) the control qubit(s), - // control on q_c = |1> (high edge) - if (var == c) { - ci++; - high = CALL(qmdd_cgate_rec, high, gate, cs, ci, t); - ci--; - } - // Not at control qubit yet, apply to both childeren. - else { - aadd_refs_spawn(SPAWN(qmdd_cgate_rec, high, gate, cs, ci, t)); - low = CALL(qmdd_cgate_rec, low, gate, cs, ci, t); - aadd_refs_push(low); - high = aadd_refs_sync(SYNC(qmdd_cgate_rec)); - aadd_refs_pop(1); - } - res = aadd_makenode(var, low, high); - - // Store not yet "root normalized" result in cache - if (cachenow) { - if (cache_put3(CACHE_QMDD_CGATE, sylvan_false, AADD_TARGET(q), - GATE_OPID_64(gate, ci, cs[0], cs[1], cs[2], t), - res)) { - sylvan_stats_count(QMDD_CGATE_CACHEDPUT); - } - } - // Multiply root amp of res with input root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; -} - -TASK_IMPL_6(QMDD, qmdd_cgate_range_rec, QMDD, q, gate_id_t, gate, BDDVAR, c_first, BDDVAR, c_last, BDDVAR, t, BDDVAR, k) -{ - // Past last control (done with "control part" of controlled gate) - if (k > c_last) { - return CALL(qmdd_gate_rec, q, gate, t); - } - - assert(c_first <= c_last); - assert(c_last < t); - - // Check cache - QMDD res; - bool cachenow = ((k % granularity) == 0); - if (cachenow) { - if (cache_get3(CACHE_QMDD_CGATE_RANGE, sylvan_false, AADD_TARGET(q), - GATE_OPID_64(gate, c_first, c_last, k, t, 0), - &res)) { - sylvan_stats_count(QMDD_CGATE_CACHED); - // Multiply root amp of result with the input root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; - } - } - - // Get top node - BDDVAR var, nextvar; - QMDD low, high; - if (k < c_first) { - // possibly skip to c_first if we are before c_first - aadd_get_topvar(q, c_first, &var, &low, &high); - assert(var <= c_first); - } - else { - // k is a control qubit, so get node with var = k (insert if skipped) - assert(c_first <= k && k <= c_last); - aadd_get_topvar(q, k, &var, &low, &high); - assert(var == k); - } - nextvar = var + 1; - - // Not at first control qubit yet, apply to both children - if (var < c_first) { - aadd_refs_spawn(SPAWN(qmdd_cgate_range_rec, high, gate, c_first, c_last, t, nextvar)); - low = CALL(qmdd_cgate_range_rec, low, gate, c_first, c_last, t, var); - aadd_refs_push(low); - high = aadd_refs_sync(SYNC(qmdd_cgate_range_rec)); - aadd_refs_pop(1); - } - // Current var is a control qubit, control on q_k = |1> (high edge) - else { - high = CALL(qmdd_cgate_range_rec, high, gate, c_first, c_last, t, nextvar); - } - res = aadd_makenode(var, low, high); - - // Store not yet "root normalized" result in cache - if (cachenow) { - if (cache_put3(CACHE_QMDD_CGATE_RANGE, sylvan_false, AADD_TARGET(q), - GATE_OPID_64(gate, c_first, c_last, k, t, 0), - res)) { - sylvan_stats_count(QMDD_CGATE_CACHEDPUT); - } - } - // Multiply root amp of result with the input root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(q), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; -} - -/*************************************************************/ - - - - - -/******************************************/ - -QMDD -qmdd_circuit_swap(QMDD qmdd, BDDVAR qubit1, BDDVAR qubit2) -{ - assert (qubit1 < qubit2); - - QMDD res; - - // CNOT - res = qmdd_cgate(qmdd, GATEID_X, qubit1, qubit2); - // upside down CNOT (equivalent) - res = qmdd_gate(res, GATEID_H, qubit1); - res = qmdd_cgate(res, GATEID_Z, qubit1, qubit2); - res = qmdd_gate(res, GATEID_H, qubit1); - // CNOT - res = qmdd_cgate(res, GATEID_X, qubit1, qubit2); - - return res; -} - -QMDD -qmdd_circuit_reverse_range(QMDD qmdd, BDDVAR first, BDDVAR last) -{ - QMDD res = qmdd; - BDDVAR a, b; - int num_qubits = (last - first) + 1; - for (int j = 0; j < (int)(num_qubits/2); j++) { - a = first + j; - b = last - j; - res = qmdd_circuit_swap(res, a, b); - } - return res; -} - -QMDD -qmdd_circuit_QFT(QMDD qmdd, BDDVAR first, BDDVAR last) -{ - int k; - QMDD res = qmdd; - BDDVAR a, b; - for (a = first; a <= last; a++) { - - // H gate on current qubit - res = qmdd_gate(res, GATEID_H, a); - - // Controlled phase gates on all qubits below - for (b = a+1; b <= last; b++) { - k = (b - a) + 1; - res = qmdd_cgate(res, GATEID_Rk(k), a, b); - } - } - - // Note that we're not swapping the qubit order in this function - - return res; -} - -QMDD -qmdd_circuit_QFT_inv(QMDD qmdd, BDDVAR first, BDDVAR last) -{ - int k; - QMDD res = qmdd; - BDDVAR a, b; - - // Note that we're not swapping the qubit order in this function - - // H gates and phase gates (but now backwards) - for (a = last + 1; a-- > first; ) { // weird for-loop because BDDVARs are unsigned - - // Controlled phase gates (negative angles this time) - for (b = last; b >= (a+1); b--){ - k = (b - a) + 1; - res = qmdd_cgate(res, GATEID_Rk_dag(k), a, b); - } - - // H on current qubit - res = qmdd_gate(res, GATEID_H, a); - } - - return res; -} - -QMDD -qmdd_circuit(QMDD qmdd, circuit_id_t circ_id, BDDVAR t1, BDDVAR t2) -{ - switch (circ_id) { // don't judge me please - case CIRCID_swap : return qmdd_circuit_swap(qmdd, t1, t2); - case CIRCID_reverse_range : return qmdd_circuit_reverse_range(qmdd, t1, t2); - case CIRCID_QFT : return qmdd_circuit_QFT(qmdd, t1, t2); - case CIRCID_QFT_inv : return qmdd_circuit_QFT_inv(qmdd, t1, t2); - default : - assert ("Invalid circuit ID" && false); - return AADD_TERMINAL; - } -} - -TASK_IMPL_6(QMDD, qmdd_ccircuit, QMDD, qmdd, circuit_id_t, circ_id, BDDVAR*, cs, uint32_t, ci, BDDVAR, t1, BDDVAR, t2) -{ - // Cache lookup - QMDD res; - bool cachenow = 1; - if (cachenow) { - if (cache_get3(CACHE_QMDD_SUBCIRC, GATE_OPID_40(0, ci, 0), qmdd, - GATE_OPID_64(circ_id, cs[0], cs[1], cs[2], t1, t2), - &res)) { - return res; - } - } - - // Get current control qubit - BDDVAR c = cs[ci]; - - // If no more control qubits, apply sub-circ here - if (c == AADD_INVALID_VAR || ci > MAX_CONTROLS) { - res = qmdd_circuit(qmdd, circ_id, t1, t2); - // the gates in qmdd_circuit already took care of multiplying the input - // root amp with normalization, so no need to do that here again - } - else { - BDDVAR var; - QMDD low, high; - aadd_get_topvar(qmdd, c, &var, &low, &high); - assert(var <= c); - - if (var == c) { - ci++; // next control qubit - high = CALL(qmdd_ccircuit, high, circ_id, cs, ci, t1, t2); - ci--; - } - else { - // recursive call to both children - aadd_refs_spawn(SPAWN(qmdd_ccircuit, high, circ_id, cs, ci, t1, t2)); - low = CALL(qmdd_ccircuit, low, circ_id, cs, ci, t1, t2); - aadd_refs_push(low); - high = aadd_refs_sync(SYNC(qmdd_ccircuit)); - aadd_refs_pop(1); - } - res = aadd_makenode(var, low, high); - // Multiply root amp of sum with input root amp - AMP new_root_amp = wgt_mul(AADD_WEIGHT(qmdd), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - } - - // Add to cache, return - if (cachenow) { - cache_put3(CACHE_QMDD_SUBCIRC, GATE_OPID_40(0, ci, 0), qmdd, - GATE_OPID_64(circ_id, cs[0], cs[1], cs[2], t1, t2), - res); - } - return res; -} - -QMDD -qmdd_all_control_phase_rec(QMDD qmdd, BDDVAR k, BDDVAR n, bool *x) -{ - assert(k < n); - - bool skipped_k = false; - aaddnode_t node; - if (AADD_TARGET(qmdd) == AADD_TERMINAL) { - skipped_k = true; - } - else { - node = AADD_GETNODE(AADD_TARGET(qmdd)); - BDDVAR var = aaddnode_getvar(node); - if(var > k) { - skipped_k = true; - } - } - - QMDD low, high; - if (skipped_k) { - // insert skipped node - low = aadd_bundle(AADD_TARGET(qmdd), AADD_ONE); - high = aadd_bundle(AADD_TARGET(qmdd), AADD_ONE); - } - else { - // case var == k (var < k shouldn't happen) - aaddnode_getchilderen(node, &low, &high); - } - - // terminal case, apply phase depending on x[k] (control k on 0 or 1) - if (k == (n-1)) { - if (x[k] == 1) { - AMP new_amp = wgt_mul(AADD_WEIGHT(high), AADD_MIN_ONE); - high = aadd_bundle(AADD_TARGET(high), new_amp); - } - else { - AMP new_amp = wgt_mul(AADD_WEIGHT(low), AADD_MIN_ONE); - low = aadd_bundle(AADD_TARGET(low), new_amp); - } - } - // non terminal case, choose low/high depending on x[k] (control k on 0 or 1) - else { - if (x[k] == 1) { - k++; // next level - high = qmdd_all_control_phase_rec(high, k, n, x); - k--; - } - else { - k++; - low = qmdd_all_control_phase_rec(low, k, n, x); - k--; - } - } - - QMDD res = aadd_makenode(k, low, high); - - // multiply by existing edge weight on qmdd - AMP new_root_amp = wgt_mul(AADD_WEIGHT(qmdd), AADD_WEIGHT(res)); - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - return res; -} - -QMDD -qmdd_all_control_phase(QMDD qmdd, BDDVAR n, bool *x) -{ - return qmdd_all_control_phase_rec(qmdd, 0, n, x); -} - -/*****************************************/ - - - - - -/**********************************************/ - -QMDD -qmdd_measure_qubit(QMDD qmdd, BDDVAR k, BDDVAR nvars, int *m, double *p) -{ - if (k == 0) return qmdd_measure_q0(qmdd, nvars, m, p); - qmdd = qmdd_circuit_swap(qmdd, 0, k); - qmdd = qmdd_measure_q0(qmdd, nvars, m, p); - qmdd = qmdd_circuit_swap(qmdd, 0, k); - return qmdd; -} - -QMDD -qmdd_measure_q0(QMDD qmdd, BDDVAR nvars, int *m, double *p) -{ - // get probabilities for q0 = |0> and q0 = |1> - double prob_low, prob_high, prob_root; - - QMDD low, high; - BDDVAR var; - aadd_get_topvar(qmdd, 0, &var, &low, &high); - - if (testing_mode) assert(qmdd_is_unitvector(qmdd, nvars)); - - // TODO: don't use doubles here but allow for mpreal ? - // (e.g. by using AMPs) - prob_low = qmdd_unnormed_prob(low, 1, nvars); - prob_high = qmdd_unnormed_prob(high, 1, nvars); - prob_root = qmdd_amp_to_prob(AADD_WEIGHT(qmdd)); - prob_low *= prob_root; - prob_high *= prob_root; - if (fabs(prob_low + prob_high - 1.0) > 1e-6) { - fprintf(stderr, "WARNING: prob sum = %.14lf\n", prob_low + prob_high); - } - - // flip a coin - float rnd = ((float)rand())/RAND_MAX; - *m = (rnd < prob_low) ? 0 : 1; - *p = prob_low; - - // produce post-measurement state - AMP norm; - if (*m == 0) { - high = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - norm = qmdd_amp_from_prob(prob_low); - } - else { - low = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - norm = qmdd_amp_from_prob(prob_high); - } - - QMDD res = aadd_makenode(0, low, high); - - AMP new_root_amp = wgt_mul(AADD_WEIGHT(qmdd), AADD_WEIGHT(res)); - new_root_amp = wgt_div(new_root_amp, norm); - - res = aadd_bundle(AADD_TARGET(res), new_root_amp); - res = qmdd_remove_global_phase(res); - return res; -} - -QMDD -qmdd_measure_all(QMDD qmdd, BDDVAR n, bool* ms, double *p) -{ - aaddnode_t node; - bool skipped; - BDDVAR var; - double prob_low, prob_high, prob_path = 1.0, prob_roots = 1.0; - - for (BDDVAR k=0; k < n; k++) { - // find relevant node (assuming it should be the next one) - skipped = false; - if (AADD_TARGET(qmdd) == AADD_TERMINAL) { - skipped = true; - } - else { - node = AADD_GETNODE(AADD_TARGET(qmdd)); - var = aaddnode_getvar(node); - assert(var >= k); - if (var > k) skipped = true; - } - QMDD low, high; - if (skipped) { - // if skipped q0 is a don't care, treat separately? - low = aadd_bundle(AADD_TARGET(qmdd), AADD_ONE); - high = aadd_bundle(AADD_TARGET(qmdd), AADD_ONE); - } - else { - aaddnode_getchilderen(node, &low, &high); - } - - prob_low = qmdd_unnormed_prob(low, k+1, n); - prob_high = qmdd_unnormed_prob(high, k+1, n); - prob_roots *= qmdd_amp_to_prob(AADD_WEIGHT(qmdd)); - prob_high = prob_high * prob_roots / prob_path; - prob_low = prob_low * prob_roots / prob_path; - - if (fabs(prob_low + prob_high - 1.0) > 1e-6) { - fprintf(stderr, "WARNING: prob sum = %.14lf\n", prob_low + prob_high); - } - - // flip a coin - float rnd = ((float)rand())/RAND_MAX; - ms[k] = (rnd < prob_low) ? 0 : 1; - - // Get next edge - qmdd = (ms[k] == 0) ? low : high; - prob_path *= (ms[k] == 0) ? prob_low : prob_high; - } - - *p = prob_path; - - //TODO: replace below ith "return qmdd_create_basis_state(n, ms);"" - QMDD low, high, prev = AADD_TERMINAL; - - for (int k = n-1; k >= 0; k--) { - if (ms[k] == 0) { - low = aadd_bundle(AADD_TARGET(prev), AADD_ONE); - high = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - } - else if (ms[k] == 1) { - low = aadd_bundle(AADD_TERMINAL, AADD_ZERO); - high = aadd_bundle(AADD_TARGET(prev), AADD_ONE); - } - // add node to unique table - prev = aadd_makenode(k, low, high); - } - return prev; -} - -// Container for disguising doubles as ints so they can go in Sylvan's cache -// (see also union "hack" in mtbdd_satcount) -typedef union { - double as_double; - uint64_t as_int; -} double_hack_t; - -TASK_IMPL_3(double, qmdd_unnormed_prob, QMDD, qmdd, BDDVAR, topvar, BDDVAR, nvars) -{ - assert(topvar <= nvars); - - if (topvar == nvars) { - assert(AADD_TARGET(qmdd) == AADD_TERMINAL); - return qmdd_amp_to_prob(AADD_WEIGHT(qmdd)); - } - - // Look in cache - bool cachenow = 1; - if (cachenow) { - uint64_t prob_bits; - if (cache_get3(CACHE_QMDD_PROB, sylvan_false, qmdd, QMDD_PARAM_PACK_16(topvar, nvars), &prob_bits)) { - sylvan_stats_count(QMDD_PROB_CACHED); - double_hack_t container = (double_hack_t) prob_bits; - return container.as_double; - } - } - - // Check if the node we want is being skipped - BDDVAR var; - QMDD low, high; - aadd_get_topvar(qmdd, topvar, &var, &low, &high); - - double prob_low, prob_high, prob_root, prob_res; // "prob" = absolute amps squared - BDDVAR nextvar = topvar + 1; - - SPAWN(qmdd_unnormed_prob, high, nextvar, nvars); - prob_low = CALL(qmdd_unnormed_prob, low, nextvar, nvars); - prob_high = SYNC(qmdd_unnormed_prob); - prob_root = qmdd_amp_to_prob(AADD_WEIGHT(qmdd)); - prob_res = prob_root * (prob_low + prob_high); - - // Put in cache and return - if (cachenow) { - double_hack_t container = (double_hack_t) prob_res; - if (cache_put3(CACHE_QMDD_PROB, sylvan_false, qmdd, QMDD_PARAM_PACK_16(topvar, nvars), container.as_int)) - sylvan_stats_count(QMDD_PROB_CACHEDPUT); - } - return prob_res; -} - -complex_t -qmdd_get_amplitude(QMDD q, bool *basis_state) -{ - complex_t res; - weight_value(aadd_getvalue(q, basis_state), &res); - return res; -} - -double -qmdd_amp_to_prob(AMP a) -{ - complex_t c; - weight_value(a, &c); - double abs = flt_sqrt( c.r*c.r + c.i*c.i ); - return (abs*abs); -} - -AMP -qmdd_amp_from_prob(double a) -{ - complex_t c; - c.r = flt_sqrt(a); - c.i = 0; - return weight_lookup(c); -} - -/*********************************************/ - - - - - -/***************************************************************/ - -QMDD -qmdd_remove_global_phase(QMDD qmdd) -{ - // Remove global phase by replacing amp of qmdd with absolute value of amp - AMP abs = wgt_abs(AADD_WEIGHT(qmdd)); - QMDD res = aadd_bundle(AADD_TARGET(qmdd), abs); - return res; -} - -void -qmdd_set_periodic_gc_nodetable(int every_n_gates) -{ - periodic_gc_nodetable = every_n_gates; -} - -/**************************************************************/ - - - - - -/***************************************************************/ - -bool qmdd_stats_logging = false; -uint32_t statslog_granularity = 1; -uint64_t statslog_buffer = 10; // TODO: remove manual buffer flushing -FILE *qmdd_logfile; -uint64_t nodes_peak = 0; -double nodes_avg = 0; -uint64_t logcounter = 0; -uint64_t logtrycounter = 0; - -void -qmdd_stats_start(FILE *out) -{ - if (out == NULL) return; - qmdd_stats_logging = true; - qmdd_logfile = out; - fprintf(qmdd_logfile, "nodes, amps\n"); - nodes_peak = 0; - logcounter = 0; - logtrycounter = 0; -} - -void -qmdd_stats_set_granularity(uint32_t g) -{ - if (g == 0) statslog_granularity = 1; - else statslog_granularity = g; -} - -void -qmdd_stats_log(QMDD qmdd) -{ - if (!qmdd_stats_logging) return; - - // only log every 'statslog_granularity' calls of this function - if (logtrycounter++ % statslog_granularity != 0) return; - - // Insert info - uint64_t num_nodes = aadd_countnodes(qmdd); - uint64_t num_amps = sylvan_edge_weights_count_entries(); - fprintf(qmdd_logfile, "%" PRIu64 ",%" PRIu64 "\n", num_nodes, num_amps); - logcounter++; - - // manually flush every 'statslog_buffer' entries - if (logcounter % statslog_buffer == 0) - fflush(qmdd_logfile); - - // peak nodes - if (num_nodes > nodes_peak) - nodes_peak = num_nodes; - - // (online) avg nodes - double a = 1.0/(double)logcounter; - double b = 1.0 - a; - nodes_avg = a * (double)num_nodes + b * nodes_avg; -} - -uint64_t -qmdd_stats_get_nodes_peak() -{ - return nodes_peak; -} - -double -qmdd_stats_get_nodes_avg() -{ - return nodes_avg; -} - -uint64_t -qmdd_stats_get_logcounter() -{ - return logtrycounter; -} - -void -qmdd_stats_finish() -{ - if (!qmdd_stats_logging) return; - fflush(qmdd_logfile); - qmdd_stats_logging = false; - nodes_peak = 0; - logcounter = 0; - logtrycounter = 0; -} - -/**************************************************************/ - - - - - -/*****************************************************************/ - -void -qmdd_set_testing_mode(bool on) -{ - testing_mode = on; -} - -bool -qmdd_is_close_to_unitvector(QMDD qmdd, BDDVAR n, double tol) -{ - bool WRITE_TO_FILE = false; - bool has_next = true; - AMP a; - bool x[n]; - for(BDDVAR k=0; k*********************************/ diff --git a/src/qsylvan_mtbdd_simulator.h b/src/qsylvan_mtbdd_simulator.h deleted file mode 100644 index c23de57..0000000 --- a/src/qsylvan_mtbdd_simulator.h +++ /dev/null @@ -1,414 +0,0 @@ -/* - * Copyright 2023 System Verification Lab, LIACS, Leiden University - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - * - */ - -// -// This code is copied in September 2023 from qsylvan_simulator.h/c made by Sebastian Brand -// - -#ifndef QSYLVAN_MTBDD_SIMULATOR_H -#define QSYLVAN_MTBDD_SIMULATOR_H - -#include -#include - -//typedef AADD QMDD; // QMDD edge (contains AMP and PTR) -//typedef AADD_WGT AMP; // edge weight index -//typedef AADD_TARG PTR; // node index - -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - -// ************************************************************* - -/** - * TODO: description of params. - * - * This could be a constructor of the QSylvan class "QSimulator()". - * -*/ - -/* -class QSimulator() -{ - size_t weight_table_size = DEFAULT_WEIGHT_TABLE_SIZE; - double weight_table_tolerance = DEFAULT_TABLE_TOLERANCE; - int edge_weight_backend = DEFAULT_WEIGHT_BACKEND; - int norm_strategy = DEFAULT_NORM_STRATEGY; - - QSimulator(size_t weight_table_size, double weight_table_tolerance, int edge_weight_backend, int norm_strategy) - { - qsylvan_init_simulator(weight_table_size, weight_table_tolerance, edge_weight_backend, norm_strategy); - qsylvan_init_defaults(size_t wgt_tab_size); - } - -} -*/ - -void qsylvan_init_simulator(size_t wgt_tab_size, double wgt_tab_tolerance, int edge_weigth_backend, int norm_strat); -void qsylvan_init_defaults(size_t wgt_tab_size); - -/************************************************************* - - -/****************************************************** - -/** - * Creates a QMDD for an n-qubit state |00...0>. - * - * @param n Number of qubits. - * - * @return A QMDD encoding the n-qubit state |00..0>. - */ -QMDD qmdd_create_all_zero_state(BDDVAR n); - -/** - * Creates a QMDD for an n-qubit state |x>. - * - * @param n Number of qubits. - * @param x A bitstring x \in {0,1}^n. - * - * @return A QMDD encoding of the n-qubit state |x>. - */ -QMDD qmdd_create_basis_state(BDDVAR n, bool* x); - -/** - * Creates a QMDD representing the matrix I \tensor I \tensor ... \tensor I. - * - * @param n Number of qubits. - * - * @return A QMDD encoding of I \tensor I \tensor ... \tensor I - */ -QMDD qmdd_create_all_identity_matrix(BDDVAR n); - -/** - * Creates a QMDD matrix which applies gate U to qubit t and I to all others. - * - * @param n Total number of qubits. - * @param t Target qubit for given gate. - * @param gateid Gate ID of predefined single qubit gate U. - * - * @return A QMDD encoding of I_0 \tensor ... U_t ... \tensor I_{n-1} - */ -QMDD qmdd_create_single_qubit_gate(BDDVAR n, BDDVAR t, gate_id_t gateid); - -/** - * Creates a QMDD matrix which applies the given list of n gates to n qubits. - * - * @param n Total number of qubits. - * @param gateids List of length n of gate ID of predefined single qubit gates. - * - * @return A QMDD encoding of U_0 \tensor U_1 \tensor ... \tensor U_{n-1} - */ -QMDD qmdd_create_single_qubit_gates(BDDVAR n, gate_id_t *gateid); - -/** - * Creates a QMDD matrix which applies single qubit gate U to all qubits. - * - * @param n Total number of qubits. - * @param gateid Gate ID of predefined single qubit gate U. - * - * @return A QMDD encoding of U^{\tensor n} - */ -QMDD qmdd_create_single_qubit_gates_same(BDDVAR n, gate_id_t gateid); - -/** - * Creates a QMDD matrix which does gateid(c,t) and I on all other qubits. - * - * @param n Total number of qubits. - * @param c Control qubit (for now we assume/need that c < t). - * @param t Target qubit for given gate. - * @param gateid Gate ID of predefined single qubit gate U. - * - * @return A matrix QMDD encoding of the controlled gate on given qubits. - */ -QMDD qmdd_create_controlled_gate(BDDVAR n, BDDVAR c, BDDVAR t, gate_id_t gateid); - -/** - * Creates a controlled-`gateid` gate which acts on the qubits as specified by - * `c_options`. - * - * @param n Total number of qubits. - * @param c_options Array of length n with option for each qubit k: { - * -1 : ignore qubit k (apply I), - * 0 : control on q_k = |0> - * 1 : control on q_k = |1> - * 2 : target qubit } - * @param gateid Gate ID of predefined single qubit gate U. - * - * @return A matrix QMDD encoding of the multi-controlled gate on given qubits. - */ -#define qmdd_create_multi_cgate(n,c_options,gateid) qmdd_create_multi_cgate_rec(n,c_options,gateid,0) -QMDD qmdd_create_multi_cgate_rec(BDDVAR n, int *c_options, gate_id_t gateid, BDDVAR k); - -/** - * Creates an n-qubit controlled Z gate, controlled on all qubits. - * - * @param n Number of qubits. - * @param x A bitstring x whether to control qubit k on x[k] = 0 or x[k] = 1. - * - * @return A matrix QMDD encoding of an n-qubit controlled Z gate. - */ -QMDD qmdd_create_all_control_phase(BDDVAR n, bool *x); - -/*****************************************************/ - - - - - -/**********************************************/ - -/** - * Computational basis measurement on qubit q_k. - * - * @param qmdd A QMDD encoding of some n qubit state. - * @param k Which qubit to measure. - * @param m Return of measurement outcome (0 or 1). - * @param p Return of measurement probability. - * - * @return QMDD of post-measurement state corresponding to measurement outcome. - */ -QMDD qmdd_measure_qubit(QMDD qqd, BDDVAR k, BDDVAR nvars, int *m, double *p); -QMDD qmdd_measure_q0(QMDD qmdd, BDDVAR nvars, int *m, double *p); - -/** - * Computational basis measurement of all n qubits in the qmdd. - * - * @param qmdd A QMDD encoding an n qubit state |psi>. - * @param n Number of qubits. - * @param ms Array of lenght n where the measurement outcomes are put. - * @param p Return of measurement probability ||^2. - * - * @return QMDD of post-measurement state (computational basis state |ms>). - */ -QMDD qmdd_measure_all(QMDD qmdd, BDDVAR n, bool* ms, double *p); - -/** - * (Recursive) helper function for obtaining probabilities for measurements - */ -#define qmdd_unnormed_prob(qmdd, topvar, nvars) (RUN(qmdd_unnormed_prob,qmdd,topvar,nvars)) -TASK_DECL_3(double, qmdd_unnormed_prob, QMDD, BDDVAR, BDDVAR); - -/** - * Get amplitude of given basis state. - * - * @param qmdd A QMDD encoding some quantum state |psi>. - * @param basis_state A bitstring x of some computational basis state |x>. - * - * @return The amplitude . - */ -complex_t qmdd_get_amplitude(QMDD qmdd, bool *basis_state); - -/** - * Computes the probability from a given edge weight index. - * - * @param a Edge weight index - * @return The probability |value(a)|^2 - */ -double qmdd_amp_to_prob(AMP a); - -/** - * Turns a given probability into a quantum amplitude, and stores it in the - * edge weight table. - * - * @param a Probability - * @return Edge weight table index of sqrt(a) - */ -AMP qmdd_amp_from_prob(double a); - -/*********************************************/ - - - - - -/**************************************************************/ - -// For now we have at most 3 control qubits -#define MAX_CONTROLS 3 - -/* Applies given (single qubit) gate to |q>. (Wrapper function) */ -#define qmdd_gate(qmdd,gate,target) (RUN(qmdd_gate,qmdd,gate,target)) -TASK_DECL_3(QMDD, qmdd_gate, QMDD, gate_id_t, BDDVAR); - -/* Applies given controlled gate to |q>. (Wrapper function) */ -#define qmdd_cgate(qmdd,gate,c,t) (RUN(qmdd_cgate,qmdd,gate,c,t)) -#define qmdd_cgate2(qmdd,gate,c1,c2,t) (RUN(qmdd_cgate2,qmdd,gate,c1,c2,t)) -#define qmdd_cgate3(qmdd,gate,c1,c2,c3,t) (RUN(qmdd_cgate3,qmdd,gate,c1,c2,c3,t)) -TASK_DECL_4(QMDD, qmdd_cgate, QMDD, gate_id_t, BDDVAR, BDDVAR); -TASK_DECL_5(QMDD, qmdd_cgate2, QMDD, gate_id_t, BDDVAR, BDDVAR, BDDVAR); -TASK_DECL_6(QMDD, qmdd_cgate3, QMDD, gate_id_t, BDDVAR, BDDVAR, BDDVAR, BDDVAR); - -/* Applies given controlled gate to |q>. (Wrapper function) */ -#define qmdd_cgate_range(qmdd,gate,c_first,c_last,t) (RUN(qmdd_cgate_range,qmdd,gate,c_first,c_last,t)) -TASK_DECL_5(QMDD, qmdd_cgate_range, QMDD, gate_id_t, BDDVAR, BDDVAR, BDDVAR); - -/** - * Recursive implementation of applying single qubit gates - */ -#define qmdd_gate_rec(q,gate,target) (RUN(qmdd_gate_rec,q,gate,target)) -TASK_DECL_3(QMDD, qmdd_gate_rec, QMDD, gate_id_t, BDDVAR); - -/** - * Recursive implementation of applying controlled gates - */ -#define qmdd_cgate_rec(q,gate,cs,t) (RUN(qmdd_cgate_rec,q,gate,cs,0,t)) -TASK_DECL_5(QMDD, qmdd_cgate_rec, QMDD, gate_id_t, BDDVAR*, uint32_t, BDDVAR); - -/** - * Recursive implementation of applying controlled gates where the controlles - * are defined by a range 'c_first' through 'c_last'. - */ -#define qmdd_cgate_range_rec(q,gate,c_first,c_last,t) (RUN(qmdd_cgate_range_rec,q,gate,c_first,c_last,t,0)) -TASK_DECL_6(QMDD, qmdd_cgate_range_rec, QMDD, gate_id_t, BDDVAR, BDDVAR, BDDVAR, BDDVAR); - -/*************************************************************/ - - - - - -/******************************************/ - -// Circuit IDs -typedef enum circuit_id { - CIRCID_swap, - CIRCID_reverse_range, - CIRCID_QFT, - CIRCID_QFT_inv -} circuit_id_t; - -/** - * Circuit which implements a SWAP gate from single-qubit and controlled gates. - */ -QMDD qmdd_circuit_swap(QMDD qmdd, BDDVAR qubit1, BDDVAR qubit2); - -/** - * Circuit which reverses the order of the qubits in the given range. - */ -QMDD qmdd_circuit_reverse_range(QMDD qmdd, BDDVAR first, BDDVAR last); - -/** - * Executes the QFT circuit on qubits `first` through `last`. - */ -QMDD qmdd_circuit_QFT(QMDD qmdd, BDDVAR first, BDDVAR last); - -/** - * Executes the inverse QFT circuit on qubits `first` through `last`. - */ -QMDD qmdd_circuit_QFT_inv(QMDD qmdd, BDDVAR first, BDDVAR last); - -/** - * Applies the given circuit (parameters can be two targets or a range - * depending on the circuit.) - */ -QMDD qmdd_circuit(QMDD qmdd, circuit_id_t circ_id, BDDVAR t1, BDDVAR t2); - -/** - * Generalized implementation of applying controlled versions of sub-circuit - * functions defined here. - * @param circ_id CIRCID_something - * @param cs BDDVAR[] of control qubits. Needs to be length 3. If using fewer - * controls use e.g. cs = [c1, c2, AADD_INVALID_VAR] - * @param t1 BDDVAR. Parameter 1 for given circuit. - * @param t2 BDDVAR. Parameter 2 for given circuit. - */ -#define qmdd_ccircuit(qmdd, circ_id, cs, t1, t2) (RUN(qmdd_ccircuit,qmdd,circ_id,cs,0,t1,t2)); -TASK_DECL_6(QMDD, qmdd_ccircuit, QMDD, circuit_id_t, BDDVAR*, uint32_t, BDDVAR, BDDVAR); - -/** - * Applies a phase of -1 to a single basis state |x>. - * This is a CZ gate where we control on all qubits and when x_k = 0 we control - * qubit k on 0, and where x_k = 1 we control qubit k on 1. - * - * @param qmdd A QMDD encoding some quantum state |psi>. - * @param n Number of qubits. - * @param x A bitstring x of some computational basis state |x>. - * - * TODO: generalize this to control on some but not all qubits. - */ -QMDD qmdd_all_control_phase(QMDD qmdd, BDDVAR n, bool *x); - -/*****************************************/ - - - - - -/***************************************************************/ - -/** - * Removes any global phase from state vector, assumes is the root edge. - */ -QMDD qmdd_remove_global_phase(QMDD qmdd); - -/** - * Trigger for gc of node table every n gates - */ -void qmdd_set_periodic_gc_nodetable(int every_n_gates); - -/**************************************************************/ - - - - - -/***************************************************************/ - -void qmdd_stats_start(FILE *out); -void qmdd_stats_set_granularity(uint32_t g); // log every 'g' gates (default 1) -void qmdd_stats_log(QMDD qmdd); -uint64_t qmdd_stats_get_nodes_peak(); -double qmdd_stats_get_nodes_avg(); -uint64_t qmdd_stats_get_logcounter(); -void qmdd_stats_finish(); - -/**************************************************************/ - - - - - -/*********************************************************/ - -/** - * Turns on/off (expensive) sanity checks - */ -void qmdd_set_testing_mode(bool on); - -/** - * Brute force sanity check to verify a given QMDD encodes (something close to) - * a unit vector. - */ -bool qmdd_is_close_to_unitvector(QMDD qmdd, BDDVAR n, double tol); -bool qmdd_is_unitvector(QMDD qmdd, BDDVAR n); - -/** - * Get the magnitude of a given qmdd. (Should equal 1 if the qmdd represents a - * (pure) quantum state. -*/ -double qmdd_get_magnitude(QMDD qmdd, BDDVAR n); - -/********************************************************/ - - -#ifdef __cplusplus -} -#endif /* __cplusplus */ - -#endif diff --git a/src/qsylvan_simulator.c b/src/qsylvan_simulator.c index eb3eb5b..53dac5c 100644 --- a/src/qsylvan_simulator.c +++ b/src/qsylvan_simulator.c @@ -1,8 +1,23 @@ -#include +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ +#include #include - static bool testing_mode = 0; // turns on/off (expensive) sanity checks static int granularity = 1; // operation cache access granularity diff --git a/src/qsylvan_simulator.h b/src/qsylvan_simulator.h index 7bf4bbe..cc2c601 100644 --- a/src/qsylvan_simulator.h +++ b/src/qsylvan_simulator.h @@ -1,11 +1,26 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ /** -* -* This code is created by Sebastiaan Brand based on the earlier developed MTBDD method by Tom van Dijk. -* -* It is currently not used in the QSylvan top-layer, and its unit tests. -* -*/ + * This code is created by Sebastiaan Brand based on the earlier developed MTBDD method by Tom van Dijk. + * + * It is currently not used in the QSylvan top-layer, and its unit tests. + * + */ #ifndef QSYLVAN_SIMULATOR_H #define QSYLVAN_SIMULATOR_H diff --git a/src/qsylvan_simulator_mtbdd.c b/src/qsylvan_simulator_mtbdd.c new file mode 100644 index 0000000..6c79108 --- /dev/null +++ b/src/qsylvan_simulator_mtbdd.c @@ -0,0 +1,290 @@ +/* + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +#include + +#include +#include +#include +//#include + +/** + * Convert zero state vector x[] = {0.0, 0.0, ..., 0.0} to MTBDD + * + * This results in an MTBDD with mpc complex zero leafs, deepness n + * + * x0 + * x1 0.0+i0.0 + * ... + * x(n-1) + * + * 0.0+i0.0 0.0+i0.0 + * + */ +MTBDD +mtbdd_create_all_zero_state_double(BDDVAR n) +{ + bool x[n]; + for (BDDVAR k=0; k = |(0,0,0,0,1,0,0,0) = |4+1> = |5> + * + * Equivalent with: + * + * x[n] = {false, false, ..., false} + * + * This results in an MTBDD with mpc complex leafs, deepness n, + * so number of leafs 2 ** n. + * + * The vars are even {0,2,4,...} to align linear algebra operations. + * + * x0 + * + * x2 x2 + * + * .... .... + * + * x(2**(n-1)) x(2**(n-1)) x(2**(n-1)) x(2**(n-1)) + * + * 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 + * + * This will be reduced to + * + * x0 + * + * x2 0.0 + * + * / 0.0 + * + * x(2**(n-1)) + * + * 0.0 0.0 + * + * We will build the reduced mtbdd directly starting at the bottom where we place a 1 leaf. + * + * n highest var + * 0 0 + * 1 n-1 = 0 + * 2 n+0 = 2 + * 3 n+1 = 4 + * 4 n+2 = 6 + * ... + * n n+n-2 = 2n-2, n > 0 + * + */ +MTBDD +mtbdd_create_basis_state_double(BDDVAR n, bool* x) // TODO: convert double to complex_t +{ + if(n==0) + return MTBDD_ZERO; + + uint32_t var = 2*n-2; + + double double_zero = 0.0; + double double_one = 1.0; + + MTBDD zero = mtbdd_double(double_zero); + MTBDD one = mtbdd_double(double_one); + + MTBDD node = one; + + // + // Start with the least significant qubit first, x[n-1] + // + // Build a path from the bottom (place one leaf = 1.0) leaf to the root. + // + // If the qubit is 0 choose the low edge, otherwise the high edge + // + for(int i = (int)n - 1; i>=0; i--) + { + //printf("x[%d] = %d, var = %d\n", i, x[i], var); + + if(x[i] == 0) + node = mtbdd_makenode(var, node, zero); + + if(x[i] == 1) + node = mtbdd_makenode(var, zero, node); + + // var of node always even + var = var - 2; + } + + return node; +} + +MTBDD +mtbdd_create_basis_state_mpc(BDDVAR n, bool* x) +{ + if(n==0) + return MTBDD_ZERO; + + uint32_t var = 2*n-2; + + mpc_t mpc_zero; + mpc_init2(mpc_zero, MPC_PRECISION); + mpc_assign(mpc_zero, 0.0, 0.0); + + mpc_t mpc_one; + mpc_init2(mpc_one, MPC_PRECISION); + mpc_assign(mpc_one, 1.0, 0.0); + + MTBDD zero = mtbdd_makeleaf(MPC_TYPE, (size_t)mpc_zero); + MTBDD one = mtbdd_makeleaf(MPC_TYPE, (size_t)mpc_one); + + MTBDD node = one; + + // + // Start with the least significant qubit first, x[n-1] + // + // Build a path from the bottom (place one leaf = 1.0) leaf to the root. + // + // If the qubit is 0 choose the low edge, otherwise the high edge + // + for(int i = (int)n - 1; i>=0; i--) + { + //printf("x[%d] = %d, var = %d\n", i, x[i], var); + + if(x[i] == 0) + node = mtbdd_makenode(var, node, zero); + + if(x[i] == 1) + node = mtbdd_makenode(var, zero, node); + + // var of node always even + var = var - 2; + } + + mpc_clear(mpc_zero); + mpc_clear(mpc_one); + + return node; +} + +/** + * + * Create a single Qubit gate surrounded by I gates: + * + * q(n-1) ----- I(n-1) ----- + * q(n-2) ----- I(n-2) ----- + * q( . ) ----- I( . ) ----- + * q( t ) ----- G( t ) ----- + * q( . ) ----- I( . ) ----- + * q( 0 ) ----- I( 0 ) ----- + * + * q(i) = {0,1}, a qubit + * + * n is number of qubits + * t is index of qubit connected to the gate G + * + * (I(0) x ... x G(t) x ... x I(n-1) ) | q(n-1)q(n-2)...q(t)...q0> + * + * a x b is the tensor of a and b, with a and b matrices / vectors. + * + */ +MTBDD +mtbdd_create_single_gate_for_qubits_mpc(BDDVAR n, BDDVAR t, MTBDD I_dd, MTBDD G_dd) //gate_id_t gateid) +{ + MTBDD dd = I_dd; + + for(uint32_t k=0; k= 0; k--) { + + if ((unsigned int)k == t) + node = mtbdd_stack_matrix(k, gateid); + else + node = mtbdd_stack_matrix(k, GATEID_I); + } + return node; +} + +MTBDD +mtbdd_stack_matrix(BDDVAR k, gate_id_t gateid) +{ + // This function effectively does a Kronecker product gate \tensor below + BDDVAR s, t; + MTBDD u00, u01, u10, u11, low, high, res; + + // Even + uneven variable are used to encode the 4 values + s = 2 * k; + t = s + 1; + + // Matrix U = [u00 u01 + // u10 u11] encoded in a small tree + u00 = mtbdd_makeleaf(MPC_TYPE, (size_t)gates[gateid][0]); + u10 = mtbdd_makeleaf(MPC_TYPE, (size_t)gates[gateid][2]); + u01 = mtbdd_makeleaf(MPC_TYPE, (size_t)gates[gateid][1]); + u11 = mtbdd_makeleaf(MPC_TYPE, (size_t)gates[gateid][3]); + + low = mtbdd_makenode(t, u00, u10); + high = mtbdd_makenode(t, u01, u11); + res = mtbdd_makenode(s, low, high); + + // Propagate common factor on previous root amp to new root amp + // AMP new_root_amp = wgt_mul(AADD_WEIGHT(below), AADD_WEIGHT(res)); + // res = aadd_bundle(AADD_TARGET(res), new_root_amp); + return res; +} +*/ + diff --git a/src/qsylvan_simulator_mtbdd.h b/src/qsylvan_simulator_mtbdd.h new file mode 100644 index 0000000..08a5fa3 --- /dev/null +++ b/src/qsylvan_simulator_mtbdd.h @@ -0,0 +1,71 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +#ifndef QSYLVAN_SIMULATOR_MTBDD_H +#define QSYLVAN_SIMULATOR_MTBDD_H + +//#include +#include +//#include +//#include + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +//static const uint64_t MTBDD_TERMINAL = 1; +//static const BDDVAR MTBDD_INVALID_VAR = UINT8_MAX; + +/** + * Creates an MTBDD for an n-qubit state |00...0>. + * + * @param n Number of qubits. + * + * @return An MTBDD encoding the n-qubit state |00..0>. + */ +MTBDD mtbdd_create_all_zero_state_double(BDDVAR n); //TODO: extend with complex_t +MTBDD mtbdd_create_all_zero_state_mpc(BDDVAR n); + +/** + * Creates an MTBDD for an n-qubit state |x>. + * + * @param n Number of qubits. + * @param x A qubitstring x \in {0,1}^n. + * + * @return An MTBDD encoding of the n-qubit state |x>. + */ +MTBDD mtbdd_create_basis_state_double(BDDVAR n, bool* x); //TODO: extend with complex_t +MTBDD mtbdd_create_basis_state_mpc(BDDVAR n, bool* x); + +/** + * Creates an MTBDD matrix which applies gate G to qubit t and I to all others. + * + * @param n Total number of qubits. + * @param t Target qubit for given gate (index connected first qubit to gate G). + * @param I_dd mtbdd of predefined single qubit identity gate I. + * @param G_dd mtbdd of predefined single qubit gate G. + * + * @return An MTBDD encoding of I(0) x ... x G(t) x ... x I(n-1) + * + */ +MTBDD mtbdd_create_single_gate_for_qubits_mpc(BDDVAR n, BDDVAR t, MTBDD I_dd, MTBDD G_dd); + +#ifdef __cplusplus +} +#endif /* __cplusplus */ + +#endif diff --git a/src/sylvan_mpc.c b/src/sylvan_mpc.c index 376f6a1..7731e46 100644 --- a/src/sylvan_mpc.c +++ b/src/sylvan_mpc.c @@ -232,13 +232,66 @@ mpc_init() * Assign a complex number based on real and imaginair double */ void -mpc_assign(mpc_ptr complexnumber, double real, double imag) +mpc_assign(mpc_ptr z, double real, double imag) { - mpc_init2(complexnumber, MPC_PRECISION); - mpc_set_d_d(complexnumber, real, imag, MPC_ROUNDING); + mpc_init2(z, MPC_PRECISION); + mpc_set_d_d(z, real, imag, MPC_ROUNDING); return; } +/** + * Assign the constant Pi + */ +void +mpc_assign_const_pi(mpc_ptr mpc_const_pi) +{ + mpfr_t mpfr_pi; + mpfr_init2(mpfr_pi, MPC_PRECISION); + mpfr_const_pi(mpfr_pi, MPC_ROUNDING); + + mpc_init2(mpc_const_pi, MPC_PRECISION); + mpc_set_fr(mpc_const_pi, mpfr_pi, MPC_ROUNDING); + mpfr_clear(mpfr_pi); + + return; + + // + // Typical use: + // + // mpc_ptr Pi; + // mpc_assign_const_pi(Pi); + // + // ... + // + // mpc_clear(mpc_const_pi); + // +} + +/** + * Compute the sqrt of a complex number + */ +void +mpc_sqrt_assign(mpc_ptr mpc_sqrt_z, double real, double imag) +{ + mpc_ptr z = NULL; + mpc_init2(z, MPC_PRECISION); + mpc_set_ld_ld(z, real, imag, MPC_ROUNDING); + mpc_sqrt(mpc_sqrt_z, z, MPC_ROUNDING); + + return; + + // + // Typical use: + // + // mpc_ptr sqrt_2; + // mpc_sqrt_assign(sqrt_2, 2.0, 0.0); + // + // ... + // + // mpc_clear(sqrt_2); + // +} + /** * Add two complex numbers with multiprecision */ @@ -263,6 +316,18 @@ mpc_multiplication(mpc_ptr x, mpc_ptr z1, mpc_ptr z2) return; } +/** + * Divide two complex numbers with multiprecision + */ +void +mpc_divide(mpc_ptr x, mpc_ptr z1, mpc_ptr z2) +{ + mpc_init2(x, MPC_PRECISION); + mpc_div(x, z1, z2, MPC_ROUNDING); + + return; +} + /** * Compare mpc leafs, re1 == re2 and im1 == im2 */ @@ -364,8 +429,6 @@ free_matrix_array_mpc(mpc_ptr **W_arr, int n) int print_vector_array_mpc(mpc_ptr *v_arr, int n) { - return 0; // Switch printing off - if(n<0) return 1; @@ -384,8 +447,6 @@ print_vector_array_mpc(mpc_ptr *v_arr, int n) int print_matrix_array_mpc(mpc_ptr **W_arr, int n) { - return 0; // Switch printing off - if(n<0) return 1; diff --git a/src/sylvan_mpc.h b/src/sylvan_mpc.h index 72edf31..8830269 100644 --- a/src/sylvan_mpc.h +++ b/src/sylvan_mpc.h @@ -59,7 +59,19 @@ mpc_init(void); * Assign a MPC complex number */ void -mpc_assign(mpc_ptr complexnumber, double real, double imag); +mpc_assign(mpc_ptr z, double real, double imag); + +/** + * Compute the sqrt of a complex number + */ +void +mpc_sqrt_assign(mpc_ptr mpc_sqrt_z, double real, double imag); + +/** + * Assign the constant Pi + */ +void +mpc_assign_const_pi(mpc_ptr mpc_const_pi); /** * Add two complex numbers with multiprecision @@ -73,6 +85,12 @@ mpc_addition(mpc_ptr x, mpc_ptr z1, mpc_ptr z2); void mpc_multiplication(mpc_ptr x, mpc_ptr z1, mpc_ptr z2); +/** + * Divide two complex numbers with multiprecision + */ +void +mpc_divide(mpc_ptr x, mpc_ptr z1, mpc_ptr z2); + /** * Compare MPC variables, re1 == re2 and im1 == im2 */ diff --git a/src/sylvan_mtbdd.c b/src/sylvan_mtbdd.c index 00ed737..c7fac76 100644 --- a/src/sylvan_mtbdd.c +++ b/src/sylvan_mtbdd.c @@ -3862,8 +3862,6 @@ free_matrix_array(MatArr_t **W_arr, int n) int print_vector_array(VecArr_t *v_arr, int n) { - return 0; // switch off printing - if(n<0) return 1; @@ -3879,8 +3877,6 @@ print_vector_array(VecArr_t *v_arr, int n) int print_matrix_array(MatArr_t **W_arr, int n) { - return 0; // switch off printing - if(n<0) return 1; @@ -4660,6 +4656,7 @@ MTBDD mtbdd_matmat_mult_alt(MTBDD M1, MTBDD M2, int n) */ MTBDD mtbdd_matvec_mult(MTBDD M, MTBDD v, int nvars, int currentvar) { +/* int maxvar = -1; int minvar = 100; int leafcount = 0; @@ -4669,6 +4666,7 @@ MTBDD mtbdd_matvec_mult(MTBDD M, MTBDD v, int nvars, int currentvar) minvar = 100; leafcount = 0; determine_top_var_and_leafcount(v, &minvar, &maxvar, &leafcount); +*/ MTBDD result = MTBDD_ZERO; @@ -4974,7 +4972,7 @@ MTBDD mtbdd_M(uint32_t row, uint32_t col, uint32_t n, MTBDD a, MTBDD b) * Tensor product (same function for matrix or vector MTBDDs) * TODO: this function should be a Lace TASK so that it can be parallelized. */ -MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_var_a) // leaf_var_a == n +MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_depth_of_a) // leaf_depth_of_a == nr of variables in a { mtbddnode_t na = MTBDD_GETNODE(a); mtbddnode_t nb = MTBDD_GETNODE(b); @@ -4986,7 +4984,7 @@ MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_var_a) // leaf_var_a == n // Check cache MTBDD result; - if (cache_get3(CACHE_MTBDD_TENSOR, a, b, (uint64_t) leaf_var_a, &result)) { + if (cache_get3(CACHE_MTBDD_TENSOR, a, b, (uint64_t) leaf_depth_of_a, &result)) { return result; } @@ -4994,20 +4992,47 @@ MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_var_a) // leaf_var_a == n BDDVAR var; // Recurse over A first (A \tensor B != B \tensor A) if (!mtbddnode_isleaf(na)) { - low = mtbdd_tensor_prod(mtbddnode_getlow(na), b, leaf_var_a); - high = mtbdd_tensor_prod(mtbddnode_gethigh(na), b, leaf_var_a); + low = mtbdd_tensor_prod(mtbddnode_getlow(na), b, leaf_depth_of_a); + high = mtbdd_tensor_prod(mtbddnode_gethigh(na), b, leaf_depth_of_a); var = mtbddnode_getvariable(na); } else { // A is a leaf and B is not - low = mtbdd_tensor_prod(a, mtbddnode_getlow(nb), leaf_var_a); - high = mtbdd_tensor_prod(a, mtbddnode_gethigh(nb), leaf_var_a); - var = mtbddnode_getvariable(nb) + leaf_var_a; // vars in B are offset by the depth of A + low = mtbdd_tensor_prod(a, mtbddnode_getlow(nb), leaf_depth_of_a); + high = mtbdd_tensor_prod(a, mtbddnode_gethigh(nb), leaf_depth_of_a); + var = mtbddnode_getvariable(nb) + leaf_depth_of_a; // vars in B are offset by the depth of A } result = mtbdd_makenode(var, low, high); // Cache put - cache_put3(CACHE_MTBDD_TENSOR, a, b, (uint64_t) leaf_var_a, result); + cache_put3(CACHE_MTBDD_TENSOR, a, b, (uint64_t) leaf_depth_of_a, result); 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; +} + diff --git a/src/sylvan_mtbdd.h b/src/sylvan_mtbdd.h index 05d8b42..6fcf192 100644 --- a/src/sylvan_mtbdd.h +++ b/src/sylvan_mtbdd.h @@ -48,6 +48,8 @@ extern "C" { #endif /* __cplusplus */ +#include + /** * 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. @@ -1336,14 +1338,20 @@ MTBDD mtbdd_mat_tensor_prod(MTBDD M1, MTBDD M2, int n); * * @param a MTBDD encoding of a matrix or vector * @param b MTBDD encoding of a matrix or vector - * @param leaf_var_a The 'level' of the leaves of 'a'. If 'a' encodes a matrix - * over DD variables {0,1,2,3,4}, then 'leaf_var_a' should be 5. If 'a' - * encodes a vector over DD variables {0,2,4}, then 'leaf_var_a' should - * be 6. + * @param leaf_var_a The 'level' of the leaves of 'a'. + * If 'a' encodes a matrix over DD variables {0,1,2,3,4}, then 'leaf_depth_of_a' should be 5 = nr of variables. + * If 'a' encodes a vector over DD variables {0,2,4}, then 'leaf_depth_of_a' should be 6 = 2 * nr of variables. * * @returns MTBDD encoding of 'a' tensor 'b' */ -MTBDD mtbdd_tensor_prod(MTBDD a, MTBDD b, int leaf_var_a); +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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e9407ea..b3d121e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -13,6 +13,12 @@ target_link_libraries(test_mtbdd PRIVATE qsylvan::sylvan) target_compile_features(test_mtbdd PRIVATE c_std_11) target_compile_options(test_mtbdd PRIVATE -Wall -Wextra -Werror -Wno-deprecated) +add_executable(test_simulator_mtbdd) +target_sources(test_simulator_mtbdd PRIVATE test_simulator_mtbdd.c) +target_link_libraries(test_simulator_mtbdd PRIVATE qsylvan::sylvan) +target_compile_features(test_simulator_mtbdd PRIVATE c_std_11) +target_compile_options(test_simulator_mtbdd PRIVATE -Wall -Wextra -Werror -Wno-deprecated) + add_executable(test_cxx) target_sources(test_cxx PRIVATE test_cxx.cpp) target_link_libraries(test_cxx PRIVATE qsylvan::sylvan) @@ -25,9 +31,10 @@ target_link_libraries(test_zdd PRIVATE qsylvan::sylvan) target_compile_features(test_zdd PRIVATE c_std_11) target_compile_options(test_zdd PRIVATE -Wall -Wextra -Werror -Wno-deprecated) -# tests for mtb, bdd and ldd +# tests for mtbdd, bdd and ldd add_test(test_basic test_basic) add_test(test_mtbdd test_mtbdd) +add_test(test_simulator_mtbdd test_simulator_mtbdd) add_test(test_cxx test_cxx) add_test(test_zdd test_zdd) diff --git a/test/test_mtbdd.c b/test/test_mtbdd.c index 5dfd25b..f41fe6f 100644 --- a/test/test_mtbdd.c +++ b/test/test_mtbdd.c @@ -1,5 +1,4 @@ - -/* +/** * Copyright 2023 System Verification Lab, LIACS, Leiden University * * Licensed under the Apache License, Version 2.0 (the "License"); @@ -689,6 +688,11 @@ test_mtbdd_arithmic_plus_sub_times_functions_complex() mpc_clear(sub_10); mpc_clear(sub_11); + mpc_clear(value_00); + mpc_clear(value_01); + mpc_clear(value_10); + mpc_clear(value_11); + return 0; } @@ -860,6 +864,13 @@ test_mtbdd_arithmic_min_max_functions_complex() test_assert( mpc_compare( mtbdd_getvalue(mtbdd_getlow( mtbdd_gethigh( dd_max ))), (uint64_t)max_10)); test_assert( mpc_compare( mtbdd_getvalue(mtbdd_gethigh(mtbdd_gethigh( dd_max ))), (uint64_t)max_11)); + mpc_clear(value_00); + mpc_clear(value_01); + mpc_clear(value_10); + mpc_clear(value_11); + + mpc_clear(value_0); + mpc_clear(min_00); mpc_clear(min_01); mpc_clear(min_10); @@ -1368,7 +1379,7 @@ test_vector_array_to_mtbdd_double() MatArr_t **W_arr = NULL; allocate_matrix_array(&W_arr, n); mtbdd_to_matrix_array(v, n, ALTERNATE_ROW_FIRST_WISE_MODE, W_arr); - print_matrix_array(W_arr, n); + //print_matrix_array(W_arr, n); test_assert(W_arr[0][0] == 1.0); test_assert(W_arr[0][1] == 1.0); @@ -1382,7 +1393,7 @@ test_vector_array_to_mtbdd_double() allocate_matrix_array(&W_arr, n); mtbdd_to_matrix_array(v, n, ALTERNATE_ROW_FIRST_WISE_MODE, W_arr); - print_matrix_array(W_arr, n); + //print_matrix_array(W_arr, n); test_assert(W_arr[0][0] == 1.0); test_assert(W_arr[0][1] == 1.0); @@ -1395,7 +1406,7 @@ test_vector_array_to_mtbdd_double() mtbdd_to_vector_array(v, n, ALTERNATE_ROW_FIRST_WISE_MODE, v_arr); - print_vector_array(v_arr, n); + //print_vector_array(v_arr, n); test_assert(v_arr[0] == 1.0); test_assert(v_arr[1] == 2.0); @@ -1404,7 +1415,7 @@ test_vector_array_to_mtbdd_double() mtbdd_to_vector_array(v, n, ALTERNATE_ROW_FIRST_WISE_MODE, v_arr); - print_vector_array(v_arr, n); + //print_vector_array(v_arr, n); test_assert(v_arr[0] == 1.0); test_assert(v_arr[1] == 2.0); @@ -1441,7 +1452,7 @@ test_matrix_array_to_mtbdd_double() mtbdd_to_matrix_array(K, n, ALTERNATE_ROW_FIRST_WISE_MODE, K_arr); - print_matrix_array(K_arr, n); + //print_matrix_array(K_arr, n); test_assert(K_arr[0][0] == 1.0); test_assert(K_arr[0][1] == 3.0); @@ -1452,7 +1463,7 @@ test_matrix_array_to_mtbdd_double() mtbdd_to_matrix_array(K, n, ALTERNATE_ROW_FIRST_WISE_MODE, K_arr); - print_matrix_array(K_arr, n); + //print_matrix_array(K_arr, n); test_assert(K_arr[0][0] == 1.0); test_assert(K_arr[0][1] == 3.0); @@ -1492,7 +1503,7 @@ test_mtbdd_to_matrix_array_double() mtbdd_to_matrix_array(K, 1, ALTERNATE_ROW_FIRST_WISE_MODE, K_arr); - print_matrix_array(K_arr, 1); + //print_matrix_array(K_arr, 1); test_assert(K_arr[0][0] == 1.0); test_assert(K_arr[0][1] == 3.0); @@ -1511,7 +1522,7 @@ test_mtbdd_to_matrix_array_double() mtbdd_to_matrix_array(M, n, ALTERNATE_ROW_FIRST_WISE_MODE, W_arr); - print_matrix_array(W_arr, n); + //print_matrix_array(W_arr, n); test_assert(W_arr[0][0] == 1.0); // First row test_assert(W_arr[0][1] == 0.5); @@ -1535,7 +1546,7 @@ test_mtbdd_to_matrix_array_double() mtbdd_to_matrix_array(M, n, COLUMN_WISE_MODE, W_arr); - print_matrix_array(W_arr, n); + //print_matrix_array(W_arr, n); test_assert(W_arr[0][0] == 1.0); // Right Down test_assert(W_arr[0][1] == 0.5); @@ -1694,7 +1705,7 @@ test_mtbdd_matrix_kronecker_multiplication_complex() mtbdd_to_matrix_array_mpc(K, 1, ALTERNATE_ROW_FIRST_WISE_MODE, K_arr); - print_matrix_array_mpc(K_arr, 1); + //print_matrix_array_mpc(K_arr, 1); test_assert( mpc_compare( (uint64_t)K_arr[0][0], (uint64_t)K00) ); test_assert( mpc_compare( (uint64_t)K_arr[0][1], (uint64_t)K01) ); @@ -1709,7 +1720,7 @@ test_mtbdd_matrix_kronecker_multiplication_complex() mtbdd_to_matrix_array_mpc(M, n, COLUMN_WISE_MODE, W_arr); - print_matrix_array_mpc(W_arr, n); + //print_matrix_array_mpc(W_arr, n); // Row 0 mpc_t W00, W01, W02, W03; @@ -1758,7 +1769,7 @@ test_mtbdd_matrix_kronecker_multiplication_complex() mtbdd_to_matrix_array_mpc(M, n, ALTERNATE_ROW_FIRST_WISE_MODE, W_arr); - print_matrix_array_mpc(W_arr, n); + //print_matrix_array_mpc(W_arr, n); // Left Up mpc_multiplication(W00, K00, L00); @@ -1805,7 +1816,17 @@ test_mtbdd_matrix_kronecker_multiplication_complex() test_assert( mpc_compare( (uint64_t)W_arr[3][3], (uint64_t)W03 ) ); free_matrix_array_mpc(W_arr, n); - + + mpc_clear(K00); + mpc_clear(K01); + mpc_clear(K10); + mpc_clear(K11); + + mpc_clear(L00); + mpc_clear(L01); + mpc_clear(L10); + mpc_clear(L11); + mpc_clear(W00); mpc_clear(W01); mpc_clear(W02); @@ -2070,6 +2091,14 @@ test_mtbdd_matrix_vector_multiplication_complex() free_matrix_array_mpc(K_arr, n); + mpc_clear(K00); + mpc_clear(K01); + mpc_clear(K10); + mpc_clear(K11); + + mpc_clear(L00); + mpc_clear(L10); + mpc_clear(X00); mpc_clear(X10); @@ -2113,9 +2142,9 @@ test_mtbdd_matrix_matrix_multiplication_1_double() allocate_matrix_array(&W_arr, n); mtbdd_to_matrix_array(product, n, COLUMN_WISE_MODE, W_arr); - print_matrix_array(K_arr, n); - print_matrix_array(M_arr, n); - print_matrix_array(W_arr, n); + //print_matrix_array(K_arr, n); + //print_matrix_array(M_arr, n); + //print_matrix_array(W_arr, n); if(false) { printf("W[0][0]: %lf %lf\n", W_arr[0][0], K_arr[0][0] * M_arr[0][0] + K_arr[0][1] * M_arr[1][0]); @@ -2166,9 +2195,9 @@ test_mtbdd_matrix_matrix_multiplication_2_double() allocate_matrix_array(&W_arr, n); mtbdd_to_matrix_array(product, n, COLUMN_WISE_MODE, W_arr); - print_matrix_array(K_arr, n); - print_matrix_array(M_arr, n); - print_matrix_array(W_arr, n); + //print_matrix_array(K_arr, n); + //print_matrix_array(M_arr, n); + //print_matrix_array(W_arr, n); if(false) { printf("W[0][0]: %lf %lf\n", W_arr[0][0], K_arr[0][0] * M_arr[0][0] + K_arr[0][1] * M_arr[1][0]); @@ -2241,8 +2270,8 @@ test_matrix_matrix_multiplication_4x4_double() allocate_matrix_array(&M_, n); mtbdd_to_matrix_array(M1, n, ALTERNATE_ROW_FIRST_WISE_MODE, M_); - print_matrix_array(W_, n); - print_matrix_array(M_, n); + //print_matrix_array(W_, n); + //print_matrix_array(M_, n); for(int col=0; col<(1< -#include - -//#include "qsylvan.h" - -#include -#include -#include - -#include "test_assert.h" - -uint64_t MTBDD_WEIGHT_ZERO = 0; -uint64_t MTBDD_WEIGHT_ONE = 1; - -/** - * - * Basis test for MTBDD method. - * -*/ -int test_basis_state_creation_mtbdd_1() -{ - MTBDD dd0, dd1; - - bool qubits[] = {0}; - qubits[0] = 0; dd0 = mtbdd_create_basis_state(1, qubits); // |0> function in QSimulator - qubits[0] = 1; dd1 = mtbdd_create_basis_state(1, qubits); // |1> - - test_assert( mtbdd_countnodes(dd0) == 3 ); // Only MTBDD have more terminals, one default, one for 0 and one for 1, so 3 - test_assert( mtbdd_countnodes(dd1) == 3 ); - test_assert( mtbdd_is_unitvector(dd0, 1) ); - test_assert( mtbdd_is_unitvector(dd1, 1) ); - - test_assert( mtbdd_is_ordered(dd0, 1) ); // Not needed? For sanity check! - test_assert( mtbdd_is_ordered(dd1, 1) ); - - //AMP a; // Index of edge weight, integers on the leaves (terminals) - MTBDDMAP index; // ComplexIndex index; in case of classes // ITCN? C_index? Points only to index as integer to fraction - qubits[0] = 0; - index = mtbdd_getvalue(dd0); //, qubits); - test_assert(index == MTBDD_WEIGHT_ONE); - - qubits[0] = 1; - index = mtbdd_getvalue(dd0); //, qubits); - test_assert(index == MTBDD_WEIGHT_ZERO); - - qubits[0] = 0; - index = mtbdd_getvalue(dd1); //, qubits); - test_assert(index == MTBDD_WEIGHT_ZERO); - - qubits[0] = 1; - index = mtbdd_getvalue(dd1); //, qubits); - test_assert(index == MTBDD_WEIGHT_ONE); - - printf("basis state creation mtrdd 1: ok\n"); - return 0; -} - -/** - * - * - * -*/ -int test_basis_state_creation_mtbdd_2() -{ - MTBDD dd0, dd1; // Test on 3 qubits + terminal node = 4 in number - - bool qubits[] = {0, 0, 0}; - dd0 = qmdd_create_basis_state(3, qubits); // |000> - - qubits[2] = 1; - qubits[0] = 1; - dd1 = qmdd_create_basis_state(3, qubits); // |101> - - test_assert( mtbdd_countnodes(dd0) == 4 ); - test_assert( mtbdd_countnodes(dd1) == 4 ); - test_assert( mtbdd_is_unitvector(dd0, 3) ); - test_assert( mtbdd_is_unitvector(dd1, 3) ); - - test_assert( mtbdd_is_ordered(dd0, 3) ); - test_assert( mtbdd_is_ordered(dd1, 3) ); - -/* - for(qubits, n, dirac_str) - { - for(k=2; k!=0; k--) - { - if(dirac_str[k] == 1) - qubits[k] = 1 - } - } -*/ - - x3[2] = 0; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ONE); - x3[2] = 0; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 0; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 0; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); - - x3[2] = 0; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 0; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 0; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 0; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ONE); // |101> -> 2^2*1 + 2^1*0 + 2^0*1 = index 5, starts with 0 - x3[2] = 1; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - x3[2] = 1; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); - - printf("basis state creation mtrdd 2: ok\n"); - return 0; -} diff --git a/test/test_simulator_mtbdd.c b/test/test_simulator_mtbdd.c new file mode 100644 index 0000000..a5d5448 --- /dev/null +++ b/test/test_simulator_mtbdd.c @@ -0,0 +1,295 @@ +/** + * Copyright 2024 System Verification Lab, LIACS, Leiden University + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include "test_assert.h" +#include +#include + +#include "qsylvan_simulator_mtbdd.h" + +int +test_create_all_zero_state_double() // TODO: change into complex_t +{ + BDDVAR n = 2; + MTBDD a = mtbdd_create_all_zero_state_double(n); + + VecArr_t v_arr[(1 << n)]; + + mtbdd_to_vector_array(a, n, COLUMN_WISE_MODE, v_arr); + + print_vector_array(v_arr, n); + + test_assert(v_arr[0]==(double)1.0); + + return 0; +} + +int +test_create_all_zero_state_complex() +{ + BDDVAR n = 2; + MTBDD a = mtbdd_create_all_zero_state_mpc(n); + + mpc_ptr v_arr[(1 << n)]; + + mtbdd_to_vector_array_mpc(a, n, COLUMN_WISE_MODE, v_arr); + + print_vector_array_mpc(v_arr, n); + + mpc_t mpc_one; + mpc_init2(mpc_one, MPC_PRECISION); + mpc_assign(mpc_one, 1.0, 0.0); + + test_assert(mpc_compare( (uint64_t)v_arr[0], (uint64_t)mpc_one)); + + mpc_clear(mpc_one); + + return 0; +} + +int +test_create_basis_state_double() // TODO: change into complex_t +{ + BDDVAR n = 4; + + bool x[n]; + for(int i=0; i<(int)n; i++) x[i] = 0; + x[0] = 1; + + MTBDD a = mtbdd_create_basis_state_double(n, x); + + VecArr_t v_arr[1 << n]; + + mtbdd_to_vector_array(a, n, COLUMN_WISE_MODE, v_arr); + + print_vector_array(v_arr, n); + + test_assert(v_arr[(1 << 3)] == (double)1.0); + + return 0; +} + +int +test_create_basis_state_complex() +{ + BDDVAR n = 4; + + bool x[n]; + for(int i=0; i<(int)n; i++) x[i] = 0; + x[0] = 1; + + MTBDD a = mtbdd_create_basis_state_mpc(n, x); + + mpc_ptr v_arr[1 << n]; + + mtbdd_to_vector_array_mpc(a, n, COLUMN_WISE_MODE, v_arr); + + print_vector_array_mpc(v_arr, n); + + mpc_t mpc_one; + mpc_init2(mpc_one, MPC_PRECISION); + mpc_assign(mpc_one, 1.0, 0.0); + + test_assert(mpc_compare( (uint64_t)v_arr[(1 << 3)], (uint64_t)mpc_one)); + + mpc_clear(mpc_one); + + return 0; +} + +TASK_0(int, runtests) +{ + // We are not testing garbage collection + sylvan_gc_disable(); + + // Test 1 + printf("\nTesting create all zero state double.\n"); + if (test_create_all_zero_state_double()) return 1; + + printf("\nTesting create all zero state complex.\n"); + if (test_create_all_zero_state_complex()) return 1; + + // Test 2 + printf("\nTesting create basis state double.\n"); + if (test_create_basis_state_double()) return 1; + + printf("\nTesting create basis state complex.\n"); + if (test_create_basis_state_complex()) return 1; + + return 0; +} + +int main() +{ + // Standard Lace initialization with 1 worker + lace_start(1, 0); + + // Simple Sylvan initialization, also initialize BDD, MTBDD and LDD support + sylvan_set_sizes(1LL<<20, 1LL<<20, 1LL<<16, 1LL<<16); + sylvan_init_package(); + sylvan_init_mtbdd(); + + printf("Mtbdd initialization complete.\n\n"); + + uint32_t mpc_type = mpc_init(); + printf("Mtbdd mpc type initialization complete, mpc_type = %d.\n\n", mpc_type); + + test_assert(mpc_type == MPC_TYPE); + + int result = RUN(runtests); + + sylvan_quit(); + lace_stop(); + + return result; +} + + + + + + + + + +/* +//#include +//#include + +//#include "qsylvan.h" + +//#include +//#include +//#include + +//#include "test_assert.h" + +uint64_t MTBDD_WEIGHT_ZERO = 0; +uint64_t MTBDD_WEIGHT_ONE = 1; + +** + * + * Basis test for MTBDD method. + * +/ +int test_basis_state_creation_mtbdd_1() +{ + MTBDD dd0, dd1; + + bool qubits[] = {0}; + qubits[0] = 0; dd0 = mtbdd_create_basis_state(1, qubits); // |0> function in QSimulator + qubits[0] = 1; dd1 = mtbdd_create_basis_state(1, qubits); // |1> + + test_assert( mtbdd_countnodes(dd0) == 3 ); // Only MTBDD have more terminals, one default, one for 0 and one for 1, so 3 + test_assert( mtbdd_countnodes(dd1) == 3 ); + test_assert( mtbdd_is_unitvector(dd0, 1) ); + test_assert( mtbdd_is_unitvector(dd1, 1) ); + + test_assert( mtbdd_is_ordered(dd0, 1) ); // Not needed? For sanity check! + test_assert( mtbdd_is_ordered(dd1, 1) ); + + //AMP a; // Index of edge weight, integers on the leaves (terminals) + MTBDDMAP index; // ComplexIndex index; in case of classes // ITCN? C_index? Points only to index as integer to fraction + qubits[0] = 0; + index = mtbdd_getvalue(dd0); //, qubits); + test_assert(index == MTBDD_WEIGHT_ONE); + + qubits[0] = 1; + index = mtbdd_getvalue(dd0); //, qubits); + test_assert(index == MTBDD_WEIGHT_ZERO); + + qubits[0] = 0; + index = mtbdd_getvalue(dd1); //, qubits); + test_assert(index == MTBDD_WEIGHT_ZERO); + + qubits[0] = 1; + index = mtbdd_getvalue(dd1); //, qubits); + test_assert(index == MTBDD_WEIGHT_ONE); + + printf("basis state creation mtrdd 1: ok\n"); + return 0; +} + +** + * + * + * +/ +int test_basis_state_creation_mtbdd_2() +{ + MTBDD dd0, dd1; // Test on 3 qubits + terminal node = 4 in number + + bool qubits[] = {0, 0, 0}; + dd0 = qmdd_create_basis_state(3, qubits); // |000> + + qubits[2] = 1; + qubits[0] = 1; + dd1 = qmdd_create_basis_state(3, qubits); // |101> + + test_assert( mtbdd_countnodes(dd0) == 4 ); + test_assert( mtbdd_countnodes(dd1) == 4 ); + test_assert( mtbdd_is_unitvector(dd0, 3) ); + test_assert( mtbdd_is_unitvector(dd1, 3) ); + + test_assert( mtbdd_is_ordered(dd0, 3) ); + test_assert( mtbdd_is_ordered(dd1, 3) ); + + + // for(qubits, n, dirac_str) + // { + // for(k=2; k!=0; k--) + // { + // if(dirac_str[k] == 1) + // qubits[k] = 1 + // } + // } + + + x3[2] = 0; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ONE); + x3[2] = 0; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 0; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 0; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q2, x3); test_assert(a == AADD_ZERO); + + x3[2] = 0; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 0; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 0; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 0; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 0; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 0; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ONE); // |101> -> 2^2*1 + 2^1*0 + 2^0*1 = index 5, starts with 0 + x3[2] = 1; x3[1] = 1; x3[0] = 0; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + x3[2] = 1; x3[1] = 1; x3[0] = 1; a = aadd_getvalue(q3, x3); test_assert(a == AADD_ZERO); + + printf("basis state creation mtrdd 2: ok\n"); + return 0; +} +*/ \ No newline at end of file