From 56f58739355cfa688c66f9c32b984358f09e4ca1 Mon Sep 17 00:00:00 2001 From: Maurice Jamieson Date: Thu, 6 Nov 2025 10:03:33 +0000 Subject: [PATCH 1/2] Migration to QuEST v4 API --- CMakeLists.txt | 5 +- pyquest/core.pxd | 4 +- pyquest/core.pyx | 199 +++++++++------- pyquest/decoherence.pxd | 2 - pyquest/decoherence.pyx | 5 +- pyquest/gates.pyx | 6 +- pyquest/operators.pxd | 20 +- pyquest/operators.pyx | 443 ++++++++++++++++++++++++++---------- pyquest/quest_interface.pxd | 414 +++++++++++++++++++-------------- pyquest/unitaries.pxd | 11 +- pyquest/unitaries.pyx | 162 +++++++------ 11 files changed, 783 insertions(+), 488 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 20171f4..2c9fc88 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,10 @@ set(BUILD_SHARED_LIBS TRUE) add_subdirectory(QuEST) # To get proper exception/error handling, we must override `invalidQuESTInputError()` # in QuEST, which we can do by adding `pyquest/quest_exception.c` to its sources. -target_include_directories(QuEST PRIVATE ${CMAKE_SOURCE_DIR}/pyquest) +target_include_directories(QuEST PRIVATE + $ + $ +) target_sources(QuEST PRIVATE ${CMAKE_SOURCE_DIR}/pyquest/quest_exception.cpp) # For Win32 `dll`s every exported symbol must usually be marked explicitly # as such. Fortunately, CMake supports automatic export of all symbols (yay). diff --git a/pyquest/core.pxd b/pyquest/core.pxd index 17732ef..31e4826 100644 --- a/pyquest/core.pxd +++ b/pyquest/core.pxd @@ -5,7 +5,7 @@ from libc.stdint cimport uintptr_t from cpython.ref cimport PyObject, Py_XDECREF from cpython.pycapsule cimport PyCapsule_New, PyCapsule_GetPointer cimport pyquest.quest_interface as quest -from pyquest.quest_interface cimport qreal, qcomp, Complex +from pyquest.quest_interface cimport qreal, qcomp from pyquest.quest_interface cimport OP_TYPES, Qureg, QuESTEnv from pyquest.operators cimport BaseOperator, GlobalOperator from pyquest.gates cimport M @@ -30,7 +30,7 @@ cdef class Register: cdef Qureg c_register cdef object _borrowed_from cdef object _borrowers - cdef Complex _scaling_factor + cdef qcomp _scaling_factor cpdef init_blank_state(self) cpdef apply_circuit(self, Circuit circ) cpdef apply_operator(self, BaseOperator op) diff --git a/pyquest/core.pyx b/pyquest/core.pyx index 816d5fa..50ed6df 100644 --- a/pyquest/core.pyx +++ b/pyquest/core.pyx @@ -59,12 +59,18 @@ cdef class QuESTEnvironment: def __cinit__(self): """Create internals and extract environment properties.""" - self.c_env = quest.createQuESTEnv() - logger.info("Created QuEST Environment at " + hex(&self.c_env)) + try: + quest.initQuESTEnv() # v4 API: no return value, manages environment internally + except QuESTError: + # Environment already initialized - in v4 there is + # only one global environment per process) + pass + self.c_env = quest.getQuESTEnv() # v4 API: get local copy + logger.info("Initialized QuEST Environment") self._env_capsule = PyCapsule_New(&self.c_env, NULL, NULL) self._logged_registers = WeakSet() cdef char[200] env_str - quest.getEnvironmentString(self.c_env, env_str) + quest.getEnvironmentString(env_str) py_string = env_str.decode('UTF-8') prop_dict = dict([prop.split("=") for prop in py_string.split()]) self._cuda = prop_dict['CUDA'] != '0' @@ -143,33 +149,26 @@ cdef class QuESTEnvironment: processes. This is useful for printing messages and accessing the file system. """ - return self.c_env.rank + cdef int r = self.c_env.rank + return r @property def seeds(self): - seed_list = [] - cdef int k - for k in range(self.c_env.numSeeds): - seed_list.append(self.c_env.seeds[k]) - return seed_list + cdef int num_seeds = quest.getNumSeeds() + cdef unsigned[:] seed_array = np.ndarray(num_seeds, dtype=np.uint32) + quest.getSeeds(&seed_array[0]) + return seed_array.base.tolist() @seeds.setter - def seeds(self, val): - cdef int num_seeds - try: - num_seeds = len(val) - except TypeError: - val = [val] - num_seeds = 1 - cdef unsigned long *seeds = malloc(sizeof(seeds[0]) * num_seeds) - cdef int k - for k in range(num_seeds): - seeds[k] = val[k] - if seeds[k] != val[k]: - free(seeds) - raise ValueError("Seeds for QuEST may only be positive integers.") - quest.seedQuEST(&(self.c_env), seeds, num_seeds) - free(seeds) + def seeds(self, value): + cdef unsigned[:] seed_array = np.asarray(value, dtype=np.uint32) + if seed_array.size != quest.getNumSeeds(): + raise ValueError( + f"Expected {quest.getNumSeeds()} seeds, got {seed_array.size}") + quest.setSeeds(&seed_array[0], seed_array.size) + + def reset_seeds_to_default(self): + quest.setSeedsToDefault() def close_env(self): """Close the QuEST environment. @@ -188,8 +187,8 @@ cdef class QuESTEnvironment: cdef Register reg for reg in self._logged_registers: reg._destroy() - logger.info("Closing QuEST Environment at " + hex(&self.c_env)) - quest.destroyQuESTEnv(self.c_env) + logger.info("Finalizing QuEST Environment") + quest.finalizeQuESTEnv() # v4 API: no parameter needed cdef class Register: @@ -224,19 +223,16 @@ cdef class Register: self._borrowed_from = None self._borrowers = WeakSet() if num_qubits == 0: - self.c_register.numAmpsTotal = 0 + self.c_register.numAmps = 0 return if copy_reg is None: if density_matrix: - self.c_register = quest.createDensityQureg( - num_qubits, (pyquest.env).c_env) + self.c_register = quest.createDensityQureg(num_qubits) else: - self.c_register = quest.createQureg( - num_qubits, (pyquest.env).c_env) + self.c_register = quest.createQureg(num_qubits) else: copy_reg._apply_delayed_operations() - self.c_register = quest.createCloneQureg( - copy_reg.c_register, (pyquest.env).c_env) + self.c_register = quest.createCloneQureg(copy_reg.c_register) logger.info("Created quantum register at " + hex(id(self))) (pyquest.env).log_register(self) @@ -338,29 +334,48 @@ cdef class Register: def __add__(left, right): if not (isinstance(left, Register) and isinstance(right, Register)): return NotImplemented - cdef Register res_reg = Register.zero_like(left) - cdef Complex zero - zero.real = 0 - zero.imag = 0 - quest.setWeightedQureg( - (left)._scaling_factor, (left).c_register, - (right)._scaling_factor, (right).c_register, - zero, res_reg.c_register) + cdef Register left_reg = left + cdef Register right_reg = right + cdef Register res_reg = Register.zero_like(left_reg) + cdef qcomp[2] coeffs + cdef Qureg[2] quregs + + # Create result = 1.0 * left_scaling * left + 1.0 * right_scaling * right + coeffs[0] = left_reg._scaling_factor + coeffs[1] = right_reg._scaling_factor + quregs[0] = left_reg.c_register + quregs[1] = right_reg.c_register + + quest.setQuregToWeightedSum(res_reg.c_register, &(coeffs[0]), &(quregs[0]), 2) + + # Result has scaling factor 1 since it's the sum + res_reg._scaling_factor.real = 1.0 + res_reg._scaling_factor.imag = 0.0 + return res_reg def __sub__(left, right): if not (isinstance(left, Register) and isinstance(right, Register)): return NotImplemented - cdef Register res_reg = Register.zero_like(left) - cdef Complex zero, right_fac - zero.real = 0 - zero.imag = 0 - right_fac.real = -(right)._scaling_factor.real - right_fac.imag = -(right)._scaling_factor.imag - quest.setWeightedQureg( - (left)._scaling_factor, (left).c_register, - right_fac, (right).c_register, - zero, res_reg.c_register) + cdef Register left_reg = left + cdef Register right_reg = right + cdef Register res_reg = Register.zero_like(left_reg) + cdef qcomp[2] coeffs + cdef Qureg[2] quregs + + # Create result = 1.0 * left_scaling * left - 1.0 * right_scaling * right + coeffs[0] = left_reg._scaling_factor + coeffs[1].real = -right_reg._scaling_factor.real + coeffs[1].imag = -right_reg._scaling_factor.imag + quregs[0] = left_reg.c_register + quregs[1] = right_reg.c_register + + quest.setQuregToWeightedSum(res_reg.c_register, &(coeffs[0]), &(quregs[0]), 2) + + # Result has scaling factor 1 since it's the difference + res_reg._scaling_factor.real = 1.0 + res_reg._scaling_factor.imag = 0.0 + return res_reg def __getitem__(self, index): @@ -457,6 +472,7 @@ cdef class Register: cdef qreal val_imag, val_real cdef bool_t from_scalar cdef const qcomp[:] value_arr + cdef qcomp amp try: value[0] from_scalar = False @@ -468,7 +484,7 @@ cdef class Register: "Manually setting elements is only supported " "for state vectors.") if isinstance(index, slice): - start, stop, step = index.indices(self.c_register.numAmpsTotal) + start, stop, step = index.indices(self.c_register.numAmps) num_index = len(range(start, stop, step)) k = start if not from_scalar: @@ -479,7 +495,8 @@ cdef class Register: for m in range(num_index): val_real = value_arr[m].real val_imag = value_arr[m].imag - quest.setAmps(self.c_register, k, &val_real, &val_imag, 1) + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, k, &, 1) k += step except (TypeError, ValueError): for m in range(num_index): @@ -488,14 +505,16 @@ cdef class Register: # Because QuEST needs separate arrays for real # and imaginary parts, we cannot hand off a # pointer to a single numpy memoryview. Thus we - # call setAmps for each element separately. - quest.setAmps(self.c_register, k, &val_real, &val_imag, 1) + # call setQuregAmps for each element separately. + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, k, &, 1) k += step else: val_real = value.real val_imag = value.imag for m in range(num_index): - quest.setAmps(self.c_register, k, &val_real, &val_imag, 1) + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, k, &, 1) k += step else: try: @@ -504,26 +523,29 @@ cdef class Register: for m in range(num_index): val_real = value[m].real val_imag = value[m].imag - quest.setAmps(self.c_register, index[m], &val_real, &val_imag, 1) + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, index[m], &, 1) else: val_real = value.real val_imag = value.imag for m in range(num_index): - quest.setAmps(self.c_register, index[m], &val_real, &val_imag, 1) + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, index[m], &, 1) except TypeError: # Last guess is we got a scalar index. val_real = value.real val_imag = value.imag - quest.setAmps(self.c_register, index, &val_real, &val_imag, 1) + amp = val_real + val_imag * 1j + quest.setQuregAmps(self.c_register, index, &, 1) @property def is_alive(self): """Return whether the underlying QuEST structure is valid.""" - return self.c_register.numAmpsTotal > 0 + return self.c_register.numAmps > 0 # v4 API: direct field access @property def num_qubits(self): """Return the number of qubits in the register.""" - return quest.getNumQubits(self.c_register) + return self.c_register.numQubits # v4 API: direct field access @property def num_amps(self): @@ -532,7 +554,7 @@ cdef class Register: This is given by 2 ** num_qubits for pure states, and 2 ** (2 * num_qubits) for density matrices. """ - return quest.getNumAmps(self.c_register) + return self.c_register.numAmps # v4 API: direct field access @cython.boundscheck(False) @cython.wraparound(False) @@ -583,7 +605,7 @@ cdef class Register: cdef int num_qubits = arr_qubits.size cdef qreal[:] outcome_probs = np.ndarray(1 << num_qubits, dtype=np_qreal) - quest.calcProbOfAllOutcomes(&outcome_probs[0], self.c_register, + quest.calcProbsOfAllMultiQubitOutcomes(&outcome_probs[0], self.c_register, &arr_qubits[0], num_qubits) return outcome_probs.base @@ -622,7 +644,7 @@ cdef class Register: of the same type (state vector or density matrix). """ self._apply_delayed_operations() - quest.cloneQureg(other.c_register, self.c_register) + quest.setQuregToClone(other.c_register, self.c_register) cpdef void copy_from(self, Register other): """Copy the state from another ``Register``. @@ -631,7 +653,7 @@ cdef class Register: of the same type (state vector or density matrix). """ self._apply_delayed_operations() - quest.cloneQureg(self.c_register, other.c_register) + quest.setQuregToClone(self.c_register, other.c_register) def destroy_reg(self): """Free the resources of the underlying data structure. @@ -656,9 +678,9 @@ cdef class Register: """ self._apply_delayed_operations() other._apply_delayed_operations() - cdef Complex prod + cdef qcomp prod if self.c_register.isDensityMatrix: - return quest.calcDensityInnerProduct( + return quest.calcInnerProduct( self.c_register, other.c_register) else: prod = quest.calcInnerProduct(self.c_register, other.c_register) @@ -716,7 +738,7 @@ cdef class Register: created register are identical to ``other``, but the returned register is initialised to the zero product state. """ - cdef Register new_reg = Register(other.c_register.numQubitsRepresented, + cdef Register new_reg = Register(other.c_register.numQubits, other.c_register.isDensityMatrix) return new_reg @@ -730,7 +752,7 @@ cdef class Register: logger.info("Destroying quantum register at " + hex(id(self))) # Only call destroyQureg if this is a valid Qureg; # otherwise destroyQureg will segfault. - if self.c_register.numAmpsTotal == 0: + if self.c_register.numAmps == 0: logger.debug("Underlying Qureg already destroyed") return # A borrowed c_register must not be destroyed. @@ -747,9 +769,8 @@ cdef class Register: logger.debug("Transferred ownership of Qureg to Register " "at " + hex(id(new_owner))) return - quest.destroyQureg(self.c_register, - (pyquest.env).c_env) - self.c_register.numAmpsTotal = 0 + quest.destroyQureg(self.c_register) + self.c_register.numAmps = 0 @staticmethod cdef Register _create_with_borrowed_reference(Register original_reg): @@ -778,8 +799,7 @@ cdef class Register: if borrowee is None: if self._borrowed_from is not None: self.c_register = quest.createCloneQureg( - (self._borrowed_from()).c_register, - (pyquest.env).c_env) + (self._borrowed_from()).c_register) (self._borrowed_from())._unregister_borrower(self) self._borrowed_from = None return @@ -813,14 +833,21 @@ cdef class Register: self._apply_scaling() cdef void _apply_scaling(self): - cdef Complex zero + cdef qcomp zero + cdef qcomp[3] coeffs + cdef Qureg[3] quregs if self._scaling_factor.real != 1 or self._scaling_factor.imag != 0: self._ensure_no_borrow() zero.real = 0. zero.imag = 0. - quest.setWeightedQureg( - zero, self.c_register, zero, self.c_register, - self._scaling_factor, self.c_register) + # Use setQuregToWeightedSum instead of setWeightedQureg + coeffs[0] = self._scaling_factor + coeffs[1] = zero + coeffs[2] = zero + quregs[0] = self.c_register + quregs[1] = self.c_register + quregs[2] = self.c_register + quest.setQuregToWeightedSum(self.c_register, &(coeffs[0]), &(quregs[0]), 3) self._scaling_factor.real = 1 self._scaling_factor.imag = 0 @@ -833,19 +860,19 @@ cdef class Register: at position ``row`` is returned. Density matrices use both indices as the index of the amplitude to fetch. """ - cdef Complex amp + cdef qcomp amp if self.c_register.isDensityMatrix: - amp = quest.getDensityAmp(self.c_register, row, col) + amp = quest.getDensityQuregAmp(self.c_register, row, col) else: - amp = quest.getAmp(self.c_register, row) - return amp.real + 1j * amp.imag + amp = quest.getQuregAmp(self.c_register, row) # v4 API: renamed function + return amp # These specialised functions speed up the state retrieval by having # C-loops over the sliced dimensions. @cython.boundscheck(False) @cython.wraparound(False) cdef qcomp[:, :] _get_state_from_slices(self, slice row_slice, slice col_slice): - cdef size_t mat_dim = 1LL << self.c_register.numQubitsRepresented + cdef size_t mat_dim = 1LL << self.c_register.numQubits cdef size_t c_start, c_stop, num_cols, cur_col, r_start, r_stop, num_rows, cur_row, k, m cdef int c_step, r_step c_start, c_stop, c_step = col_slice.indices(mat_dim) @@ -865,7 +892,7 @@ cdef class Register: @cython.boundscheck(False) @cython.wraparound(False) cdef qcomp[:, :] _get_state_from_col_slice(self, row_index, slice col_slice): - cdef size_t mat_dim = 1LL << self.c_register.numQubitsRepresented + cdef size_t mat_dim = 1LL << self.c_register.numQubits cdef size_t start, stop, num_cols, cur_row, cur_col, k, m cdef int step start, stop, step = col_slice.indices(mat_dim) @@ -882,7 +909,7 @@ cdef class Register: @cython.boundscheck(False) @cython.wraparound(False) cdef qcomp[:, :] _get_state_from_row_slice(self, slice row_slice, col_index): - cdef size_t mat_dim = 1LL << self.c_register.numQubitsRepresented + cdef size_t mat_dim = 1LL << self.c_register.numQubits cdef size_t start, stop, num_rows, cur_row, cur_col, k, m cdef int step start, stop, step = row_slice.indices(mat_dim) diff --git a/pyquest/decoherence.pxd b/pyquest/decoherence.pxd index 98d8bd6..24a2ee6 100644 --- a/pyquest/decoherence.pxd +++ b/pyquest/decoherence.pxd @@ -1,8 +1,6 @@ from libc.stdlib cimport malloc, free cimport pyquest.quest_interface as quest from pyquest.quest_interface cimport qreal, Qureg -from pyquest.quest_interface cimport ComplexMatrix2, ComplexMatrix4, ComplexMatrixN -from pyquest.quest_interface cimport createComplexMatrixN, destroyComplexMatrixN from pyquest.operators cimport SingleQubitOperator, MultiQubitOperator, GlobalOperator from pyquest.core cimport OP_TYPES, Register diff --git a/pyquest/decoherence.pyx b/pyquest/decoherence.pyx index 8ffc669..cf42fb9 100644 --- a/pyquest/decoherence.pyx +++ b/pyquest/decoherence.pyx @@ -59,7 +59,7 @@ cdef class PauliNoise(SingleQubitOperator): self._prob_z = probs[2] cdef int apply_to(self, Qureg c_register) except -1: - quest.mixPauli(c_register, self._target, + quest.mixPaulis(c_register, self._target, self._prob_x, self._prob_y, self._prob_z) @@ -76,5 +76,4 @@ cdef class MixDensityMatrix(GlobalOperator): self._other_register = density_matrix cdef int apply_to(self, Qureg c_register) except -1: - quest.mixDensityMatrix(c_register, self._prob, - self._other_register.c_register) + quest.mixQureg(c_register, self._other_register.c_register, self._prob) diff --git a/pyquest/gates.pyx b/pyquest/gates.pyx index 571e495..e93e27e 100644 --- a/pyquest/gates.pyx +++ b/pyquest/gates.pyx @@ -97,12 +97,12 @@ cdef class M(MultiQubitOperator): self._probabilities[k] = -1 for k in range(self._num_targets): if self._force == NULL or self._force[k] == -1: - self._results[k] = quest.measureWithStats( + self._results[k] = quest.applyQubitMeasurementAndGetProb( c_register, self._targets[k], &(self._probabilities[k])) else: - self._probabilities[k] = quest.collapseToOutcome( + self._probabilities[k] = 1.0 + self._results[k] = quest.applyForcedQubitMeasurement( c_register, self._targets[k], self._force[k]) - self._results[k] = self._force[k] @property def results(self): diff --git a/pyquest/operators.pxd b/pyquest/operators.pxd index cb6172e..8336dc0 100644 --- a/pyquest/operators.pxd +++ b/pyquest/operators.pxd @@ -1,12 +1,12 @@ from libcpp cimport bool from libc.stdlib cimport malloc, calloc, free +from libc.math cimport sqrt, sin, cos from cpython.pycapsule cimport PyCapsule_GetPointer cimport pyquest.quest_interface as quest -from pyquest.quest_interface cimport qreal, qcomp, OP_TYPES, Qureg, pauliOpType -from pyquest.quest_interface cimport phaseFunc, bitEncoding -from pyquest.quest_interface cimport QuESTEnv, Complex -from pyquest.quest_interface cimport ComplexMatrix2, ComplexMatrix4, ComplexMatrixN -from pyquest.quest_interface cimport createComplexMatrixN, destroyComplexMatrixN +from pyquest.quest_interface cimport qreal, qcomp, qindex, OP_TYPES, Qureg, pauliOpType +from pyquest.quest_interface cimport QuESTEnv +from pyquest.quest_interface cimport CompMatr1, CompMatr2, CompMatr +from pyquest.quest_interface cimport DiagMatr1, DiagMatr2, DiagMatr, FullStateDiagMatr cimport numpy as np @@ -47,10 +47,12 @@ cdef class MatrixOperator(MultiQubitOperator): cdef _copy_csingle_array(self, float complex[:, :] arr) cdef _copy_cdouble_array(self, double complex[:, :] arr) cdef _copy_clongdouble_array(self, long double complex[:, :] arr) + cdef qcomp _get_matrix_element(self, size_t k, size_t n) + cdef void _set_matrix_element(self, size_t k, size_t n, qcomp val) cdef class DiagonalOperator(GlobalOperator): - cdef quest.DiagonalOp _diag_op + cdef FullStateDiagMatr _diag_op cdef class PauliProduct(GlobalOperator): @@ -74,9 +76,9 @@ cdef class TrotterCircuit(GlobalOperator): cdef class PhaseFunc(GlobalOperator): - cdef phaseFunc _phase_func_type + cdef int _phase_func_type cdef bool _is_poly - cdef bitEncoding _bit_encoding + cdef int _bit_encoding cdef int _num_overrides cdef long long int *_override_inds cdef qreal *_override_phases @@ -88,6 +90,8 @@ cdef class PhaseFunc(GlobalOperator): cdef int *_num_terms_per_reg cdef qreal *_coeffs cdef qreal *_exponents + cdef qcomp _compute_for_indices_c(self, qindex* indices) + cdef int apply_to(self, Qureg c_register) except -1 cdef class QFT(MultiQubitOperator): diff --git a/pyquest/operators.pyx b/pyquest/operators.pyx index 4afa53d..04e78b4 100644 --- a/pyquest/operators.pyx +++ b/pyquest/operators.pyx @@ -77,23 +77,19 @@ cdef class BaseOperator: except NotImplementedError: raise ValueError( "Number of qubits must be given for this operator.") - cdef QuESTEnv *env_ptr = PyCapsule_GetPointer( - pyquest.env.env_capsule, NULL) - cdef Qureg tmp_reg = quest.createQureg( - num_qubits, env_ptr[0]) + cdef Qureg tmp_reg = quest.createQureg(num_qubits) cdef long long k, m cdef long long mat_dim = 1LL << num_qubits - cdef Complex amp + cdef qcomp amp cdef qcomp[:, :] res_mat = np.ndarray( (mat_dim, mat_dim), dtype=pyquest.core.np_qcomp) for k in range(mat_dim): quest.initClassicalState(tmp_reg, k) self.apply_to(tmp_reg) for m in range(mat_dim): - amp = quest.getAmp(tmp_reg, m) - res_mat[m, k].real = amp.real - res_mat[m, k].imag = amp.imag - quest.destroyQureg(tmp_reg, env_ptr[0]) + amp = quest.getQuregAmp(tmp_reg, m) + res_mat[m, k] = amp + quest.destroyQureg(tmp_reg) return res_mat.base cdef int apply_to(self, Qureg c_register) except -1: @@ -291,13 +287,7 @@ cdef class MatrixOperator(MultiQubitOperator): # _matrix attribute, so they are freed togehter with _matrix. if self._num_targets > 2 or self._num_controls > 0: if self._matrix != NULL: - destroyComplexMatrixN((self._matrix)[0]) - else: - # The arrays of the row-pointers for 1 and 2 qubits - # are allocated manually, because the arrays themselves - # are in the _matrix variable. - free(self._real) - free(self._imag) + quest.destroyCompMatr((self._matrix)[0]) free(self._matrix) def __repr__(self): @@ -306,11 +296,15 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t matrix_dim = 1 # Assigning a 1 first prevents integer overflows of the bit shift matrix_dim = matrix_dim << self._num_targets cdef size_t i, j + cdef qcomp element res += "\n array(\n [" for i in range(matrix_dim): res += "[" for j in range(matrix_dim): - res += f'{self._real[i][j]:.15f}{self._imag[i][j]:+.15f}j, ' + element = self._get_matrix_element(i, j) + # Convert qcomp to Python complex for formatting + py_complex = element + res += f'{py_complex.real:.15f}{py_complex.imag:+.15f}j, ' res = res[:-2] + "],\n " res = res[:-11] + "])" if self._num_controls > 0: @@ -322,58 +316,64 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t mat_dim = 1 mat_dim = mat_dim << self._num_targets cdef size_t k, n + cdef qcomp element cdef qcomp[:, :] np_mat = np.ndarray( (mat_dim, mat_dim), dtype=pyquest.core.np_qcomp) for k in range(mat_dim): for n in range(mat_dim): - np_mat[k, n] = self._real[k][n] + 1j * self._imag[k][n] + element = self._get_matrix_element(k, n) + np_mat[k, n] = element return np_mat.base cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: if self._num_targets == 1: - quest.applyMatrix2( + quest.leftapplyCompMatr1( c_register, self._targets[0], - (self._matrix)[0]) + (self._matrix)[0]) elif self._num_targets == 2: - quest.applyMatrix4( + quest.leftapplyCompMatr2( c_register, self._targets[0], self._targets[1], - (self._matrix)[0]) + (self._matrix)[0]) else: - quest.applyMatrixN( + quest.leftapplyCompMatr( c_register, self._targets, self._num_targets, - (self._matrix)[0]) + (self._matrix)[0]) else: - quest.applyMultiControlledMatrixN( + quest.applyMultiControlledCompMatr( c_register, self._controls, self._num_controls, self._targets, - self._num_targets, (self._matrix)[0]) + self._num_targets, (self._matrix)[0]) cdef _create_array_property(self): cdef size_t matrix_dim = 1 # Assigning a 1 first prevents integer overflows of the bit shift matrix_dim = matrix_dim << self._num_targets - cdef size_t k - # The only cases where core QuEST supports ComplexMatrix2 - # or ComplexMatrix4 for generic matrices are non-controlled - # cases. + # The only cases where core QuEST supports CompMatr1/CompMatr2 + # for generic matrices are non-controlled cases. if self._num_targets == 1 and self._num_controls == 0: - self._matrix = malloc(sizeof(ComplexMatrix2)) - self._real = malloc(matrix_dim * sizeof(self._real[0])) - self._imag = malloc(matrix_dim * sizeof(self._imag[0])) - for k in range(matrix_dim): - self._real[k] = (self._matrix).real[k] - self._imag[k] = (self._matrix).imag[k] + self._matrix = malloc(sizeof(quest.CompMatr1)) elif self._num_targets == 2 and self._num_controls == 0: - self._matrix = malloc(sizeof(ComplexMatrix4)) - self._real = malloc(matrix_dim * sizeof(self._real[0])) - self._imag = malloc(matrix_dim * sizeof(self._imag[0])) - for k in range(matrix_dim): - self._real[k] = (self._matrix).real[k] - self._imag[k] = (self._matrix).imag[k] + self._matrix = malloc(sizeof(quest.CompMatr2)) + else: + self._matrix = malloc(sizeof(quest.CompMatr)) + (self._matrix)[0] = quest.createCompMatr(self._num_targets) + + cdef qcomp _get_matrix_element(self, size_t k, size_t n): + # Get matrix element accounting for different matrix types. + if self._num_targets == 1: + return (self._matrix).elems[k][n] + elif self._num_targets == 2: + return (self._matrix).elems[k][n] else: - self._matrix = malloc(sizeof(ComplexMatrixN)) - (self._matrix)[0] = createComplexMatrixN(self._num_targets) - self._real = (self._matrix).real - self._imag = (self._matrix).imag + return (self._matrix).cpuElems[k][n] + + cdef void _set_matrix_element(self, size_t k, size_t n, qcomp val): + # Set matrix element accounting for different matrix types. + if self._num_targets == 1: + (self._matrix).elems[k][n] = val + elif self._num_targets == 2: + (self._matrix).elems[k][n] = val + else: + (self._matrix).cpuElems[k][n] = val cdef _numpy_array_to_matrix_attribute(self, np.ndarray arr): # For typed memoryviews we need to call different methods @@ -401,8 +401,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m].real - self._imag[k][m] = arr[k, m].imag + self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) @cython.boundscheck(False) @cython.wraparound(False) @@ -410,8 +409,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m] - self._imag[k][m] = 0. + self._set_matrix_element(k, m, (arr[k, m] + 0j)) @cython.boundscheck(False) @cython.wraparound(False) @@ -419,8 +417,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m] - self._imag[k][m] = 0. + self._set_matrix_element(k, m, (arr[k, m] + 0j)) @cython.boundscheck(False) @cython.wraparound(False) @@ -428,8 +425,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m] - self._imag[k][m] = 0. + self._set_matrix_element(k, m, (arr[k, m] + 0j)) @cython.boundscheck(False) @cython.wraparound(False) @@ -437,8 +433,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m].real - self._imag[k][m] = arr[k, m].imag + self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) @cython.boundscheck(False) @cython.wraparound(False) @@ -446,8 +441,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m].real - self._imag[k][m] = arr[k, m].imag + self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) @cython.boundscheck(False) @cython.wraparound(False) @@ -455,8 +449,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._real[k][m] = arr[k, m].real - self._imag[k][m] = arr[k, m].imag + self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) cdef class PauliSum(GlobalOperator): @@ -475,11 +468,10 @@ cdef class TrotterCircuit(GlobalOperator): cdef class DiagonalOperator(GlobalOperator): def __cinit__(self, num_qubits, diag_elements=None): - cdef QuESTEnv *env_ptr = PyCapsule_GetPointer( - pyquest.env.env_capsule, NULL) - cdef qreal[:] real_el, imag_el + cdef qcomp[:] qcomp_el + cdef qindex num_elems self.TYPE = OP_TYPES.OP_DIAGONAL - self._diag_op = quest.createDiagonalOp(num_qubits, env_ptr[0]) + self._diag_op = quest.createFullStateDiagMatr(num_qubits) if diag_elements is not None: # We'll try to handle anything that can be an ndarray by # just giving it to the np.array constructor. @@ -492,22 +484,13 @@ cdef class DiagonalOperator(GlobalOperator): if diag_elements.size != 1 << self._diag_op.numQubits: raise ValueError("'diag_elements' must have length of " "2**num_qubits") - # FIXME This possibly creates copies of real and imaginary - # parts (QuEST requires contiguous arrays for real and - # imag separately), so doing it in chunks or write - # straight to memory would be better. - real_el = np.require( - diag_elements.real, dtype=pyquest.core.np_qreal, - requirements='CW') - imag_el = np.require( - diag_elements.imag, dtype=pyquest.core.np_qreal, - requirements='CW') - quest.initDiagonalOp(self._diag_op, &real_el[0], &imag_el[0]) + # Convert Python complex array to qcomp array + num_elems = diag_elements.size + qcomp_el = np.ascontiguousarray(diag_elements, dtype=np.complex128) + quest.setFullStateDiagMatr(self._diag_op, 0, &qcomp_el[0], num_elems) def __dealloc__(self): - cdef QuESTEnv *env_ptr = PyCapsule_GetPointer( - pyquest.env.env_capsule, NULL) - quest.destroyDiagonalOp(self._diag_op, env_ptr[0]) + quest.destroyFullStateDiagMatr(self._diag_op) def __repr__(self): return ("<" + type(self).__name__ + " on " @@ -515,34 +498,34 @@ cdef class DiagonalOperator(GlobalOperator): + hex(id(self)) + ">") cdef int apply_to(self, Qureg c_register) except -1: - quest.applyDiagonalOp(c_register, self._diag_op) + quest.leftapplyFullStateDiagMatr(c_register, self._diag_op) class _BitEncoding(enum.IntEnum): - UNSIGNED = quest.bitEncoding.UNSIGNED - TWOS_COMPLEMENT = quest.bitEncoding.TWOS_COMPLEMENT + UNSIGNED = 0 # Standard unsigned binary interpretation + TWOS_COMPLEMENT = 1 # Two's complement signed interpretation class _PhaseFuncType(enum.IntEnum): - # Careful when adding new functions to this enum, the PhaseFunc - # constructor, __repr__, and inverse rely on the "SCALED" and - # "INVERSE" naming conventions to determine the correct structure - # of parameters to pass to QuEST. - NORM = quest.phaseFunc.NORM - SCALED_NORM = quest.phaseFunc.SCALED_NORM - INVERSE_NORM = quest.phaseFunc.INVERSE_NORM - SCALED_INVERSE_NORM = quest.phaseFunc.SCALED_INVERSE_NORM - SCALED_INVERSE_SHIFTED_NORM = quest.phaseFunc.SCALED_INVERSE_SHIFTED_NORM - PRODUCT = quest.phaseFunc.PRODUCT - SCALED_PRODUCT = quest.phaseFunc.SCALED_PRODUCT - INVERSE_PRODUCT = quest.phaseFunc.INVERSE_PRODUCT - SCALED_INVERSE_PRODUCT = quest.phaseFunc.SCALED_INVERSE_PRODUCT - DISTANCE = quest.phaseFunc.DISTANCE - SCALED_DISTANCE = quest.phaseFunc.SCALED_DISTANCE - INVERSE_DISTANCE = quest.phaseFunc.INVERSE_DISTANCE - SCALED_INVERSE_DISTANCE = quest.phaseFunc.SCALED_INVERSE_DISTANCE - SCALED_INVERSE_SHIFTED_DISTANCE = quest.phaseFunc.SCALED_INVERSE_SHIFTED_DISTANCE - EXPONENTIAL_POLYNOMIAL = enum.auto() + # Phase function types available in QuEST v4. + # These are mapped to v4's function-based phase API via callbacks. + # The naming conventions for "SCALED", "INVERSE", and "SHIFTED" are used + # by PhaseFunc to determine parameter structures for QuEST calls. + NORM = 0 + SCALED_NORM = 1 + INVERSE_NORM = 2 + SCALED_INVERSE_NORM = 3 + SCALED_INVERSE_SHIFTED_NORM = 4 + PRODUCT = 5 + SCALED_PRODUCT = 6 + INVERSE_PRODUCT = 7 + SCALED_INVERSE_PRODUCT = 8 + DISTANCE = 9 + SCALED_DISTANCE = 10 + INVERSE_DISTANCE = 11 + SCALED_INVERSE_DISTANCE = 12 + SCALED_INVERSE_SHIFTED_DISTANCE = 13 + EXPONENTIAL_POLYNOMIAL = 14 # Extension types do not allow nested classes, so they are defined as @@ -554,6 +537,26 @@ _BitEncoding.__name__ = "BitEncoding" _PhaseFuncType.__name__ = "FuncType" +# Global reference to current PhaseFunc instance for callback routing +# This is needed because QuEST v4 uses bare C function pointers. +cdef object _current_phase_func = None + + +# Helper function to convert a real phase value to qcomp (exp(1.0j * phase)) +cdef qcomp _phase_to_qcomp(qreal phase_value): + return cos(phase_value) + 1.0j * sin(phase_value) + + +# Module-level callback function invoked by QuEST for each basis state +cdef qcomp _compute_phase_callback(qindex* indices): + # Callback function passed to QuEST v4's setFullStateDiagMatrFromMultiVarFunc. + global _current_phase_func + if _current_phase_func is None: + raise RuntimeError("PhaseFunc callback invoked with no active instance") + # Call the C-level compute method directly + return (_current_phase_func)._compute_for_indices_c(indices) + + cdef class PhaseFunc(GlobalOperator): BitEncoding = _BitEncoding @@ -941,28 +944,222 @@ cdef class PhaseFunc(GlobalOperator): constructor_args['shifts'] = self.shifts return PhaseFunc(**constructor_args) - cdef int apply_to(self, Qureg c_register) except -1: - pass - if self._is_poly: - if self._num_regs == 0: - quest.applyPhaseFuncOverrides( - c_register, self._qubits_in_regs, self._num_qubits_per_reg[0], - self._bit_encoding, self._coeffs, self._exponents, - self._num_terms, self._override_inds, - self._override_phases, self._num_overrides) + cdef qcomp _compute_for_indices_c(self, qindex* indices): + # Compute the phase value for a given set of basis state indices. + # This is called by QuEST v4 for each basis state when initializing + # the diagonal matrix. Routes to the appropriate phase function + # implementation based on self._phase_func_type. + cdef qreal phase_value = 0.0 + cdef qreal norm_val, product_val, distance_val, shift_val + cdef int i, j, k + cdef qreal scaled_factor = 1.0 + cdef qreal divergence_val = 0.0 + + # Compute base phase value depending on function type + func_type = _PhaseFuncType(self._phase_func_type) + func_name = func_type.name + + if func_type == _PhaseFuncType.NORM: + # sqrt(sum of squared indices) + norm_val = 0.0 + for i in range(self._num_regs): + norm_val += indices[i] * indices[i] + phase_value = sqrt(norm_val) if norm_val > 0 else 0.0 + elif func_type == _PhaseFuncType.SCALED_NORM: + # Scaling * sqrt(sum of squared indices) + norm_val = 0.0 + for i in range(self._num_regs): + norm_val += indices[i] * indices[i] + phase_value = sqrt(norm_val) if norm_val > 0 else 0.0 + elif func_type == _PhaseFuncType.INVERSE_NORM: + # 1 / sqrt(sum of squared indices) + norm_val = 0.0 + for i in range(self._num_regs): + norm_val += indices[i] * indices[i] + if norm_val > 0: + phase_value = 1.0 / sqrt(norm_val) else: - quest.applyMultiVarPhaseFuncOverrides( - c_register, self._qubits_in_regs, self._num_qubits_per_reg, - self._num_regs, self._bit_encoding, self._coeffs, - self._exponents, self._num_terms_per_reg, - self._override_inds, self._override_phases, - self._num_overrides) - else: - quest.applyParamNamedPhaseFuncOverrides( - c_register, self._qubits_in_regs, self._num_qubits_per_reg, - self._num_regs, self._bit_encoding, self._phase_func_type, - self._parameters, self._num_parameters, self._override_inds, - self._override_phases, self._num_overrides) + phase_value = divergence_val + elif func_type == _PhaseFuncType.SCALED_INVERSE_NORM: + # Scaling / sqrt(sum of squared indices) + norm_val = 0.0 + for i in range(self._num_regs): + norm_val += indices[i] * indices[i] + if norm_val > 0: + phase_value = 1.0 / sqrt(norm_val) + else: + phase_value = divergence_val + elif func_type == _PhaseFuncType.SCALED_INVERSE_SHIFTED_NORM: + # Scaling / sqrt(sum of (indices[i] - shift[i])^2) + norm_val = 0.0 + for i in range(self._num_regs): + shift_val = self._parameters[2 + i] if i + 2 < self._num_parameters else 0.0 + norm_val += (indices[i] - shift_val) ** 2 + if norm_val > 0: + phase_value = 1.0 / sqrt(norm_val) + else: + phase_value = divergence_val + elif func_type == _PhaseFuncType.PRODUCT: + # Product of all indices + product_val = 1.0 + for i in range(self._num_regs): + product_val *= indices[i] + phase_value = product_val + elif func_type == _PhaseFuncType.SCALED_PRODUCT: + # Scaling * product of indices + product_val = 1.0 + for i in range(self._num_regs): + product_val *= indices[i] + phase_value = product_val + elif func_type == _PhaseFuncType.INVERSE_PRODUCT: + # 1 / product of indices + product_val = 1.0 + for i in range(self._num_regs): + product_val *= indices[i] + if product_val != 0: + phase_value = 1.0 / product_val + else: + phase_value = divergence_val + elif func_type == _PhaseFuncType.SCALED_INVERSE_PRODUCT: + # Scaling / product of indices + product_val = 1.0 + for i in range(self._num_regs): + product_val *= indices[i] + if product_val != 0: + phase_value = 1.0 / product_val + else: + phase_value = divergence_val + + elif func_type == _PhaseFuncType.DISTANCE: + # Distance between pairs of indices + if self._num_regs >= 2: + distance_val = 0.0 + for i in range(self._num_regs - 1): + distance_val += (indices[i] - indices[i + 1]) ** 2 + phase_value = sqrt(distance_val) if distance_val > 0 else 0.0 + else: + phase_value = 0.0 + + elif func_type == _PhaseFuncType.SCALED_DISTANCE: + # Scaling * distance between indices + if self._num_regs >= 2: + distance_val = 0.0 + for i in range(self._num_regs - 1): + distance_val += (indices[i] - indices[i + 1]) ** 2 + phase_value = sqrt(distance_val) if distance_val > 0 else 0.0 + else: + phase_value = 0.0 + elif func_type == _PhaseFuncType.INVERSE_DISTANCE: + # 1 / distance between indices + if self._num_regs >= 2: + distance_val = 0.0 + for i in range(self._num_regs - 1): + distance_val += (indices[i] - indices[i + 1]) ** 2 + if distance_val > 0: + phase_value = 1.0 / sqrt(distance_val) + else: + phase_value = divergence_val + else: + phase_value = divergence_val + elif func_type == _PhaseFuncType.SCALED_INVERSE_DISTANCE: + # SCALED_INVERSE_DISTANCE: scaling / distance between indices + if self._num_regs >= 2: + distance_val = 0.0 + for i in range(self._num_regs - 1): + distance_val += (indices[i] - indices[i + 1]) ** 2 + if distance_val > 0: + phase_value = 1.0 / sqrt(distance_val) + else: + phase_value = divergence_val + else: + phase_value = divergence_val + + elif func_type == _PhaseFuncType.SCALED_INVERSE_SHIFTED_DISTANCE: + # SCALED_INVERSE_SHIFTED_DISTANCE: scaling / sqrt(sum of (indices[i] - shift[i])^2) + if self._num_regs >= 2: + distance_val = 0.0 + for i in range(self._num_regs - 1): + shift_val = self._parameters[2 + i] if i + 2 < self._num_parameters else 0.0 + distance_val += (indices[i] - indices[i + 1] - shift_val) ** 2 + if distance_val > 0: + phase_value = 1.0 / sqrt(distance_val) + else: + phase_value = divergence_val + else: + phase_value = divergence_val + elif func_type == _PhaseFuncType.EXPONENTIAL_POLYNOMIAL: + # Sum of c_k * (indices[r_k])^{e_k} + phase_value = 0.0 + term_idx = 0 + for i in range(self._num_regs): + num_terms_in_reg = self._num_terms_per_reg[i] + for j in range(num_terms_in_reg): + coeff = self._coeffs[term_idx] + exponent = self._exponents[term_idx] + if indices[i] == 0 and exponent < 0: + # Handle negative exponents of zero + phase_value += 0.0 # or raise error? + else: + phase_value += coeff * (indices[i] ** exponent) + term_idx += 1 + + if "SCALED" in func_name: + # First parameter is always the scaling factor + phase_value *= self._parameters[0] + + # Apply inverse if "INVERSE" in function name + if "INVERSE" in func_name: + # Divergence override is always after scaling (if present) + div_idx = 1 if "SCALED" in func_name else 0 + divergence_val = self._parameters[div_idx] if div_idx < self._num_parameters else 0.0 + if phase_value != 0: + phase_value = 1.0 / phase_value + else: + phase_value = divergence_val + # Check for manual override + if self._num_overrides: + for override_idx in range(self._num_overrides): + # Check if current indices match this override + match = True + for i in range(self._num_regs): + override_ind = self._override_inds[override_idx * self._num_regs + i] + if indices[i] != override_ind: + match = False + break + if match: + # Use override phase instead + phase_value = self._override_phases[override_idx] + break + # Return as exp(1.0j * phase_value) + return _phase_to_qcomp(phase_value) + + cdef int apply_to(self, Qureg c_register) except -1: + global _current_phase_func + # Set global reference so callback can reach this instance + _current_phase_func = self + + cdef FullStateDiagMatr diag_matr + try: + # Create diagonal matrix + diag_matr = quest.createFullStateDiagMatr(c_register.numQubits) + + # Populate matrix using v4's function-based API + quest.setFullStateDiagMatrFromMultiVarFunc( + diag_matr, + _compute_phase_callback, # Function pointer to callback + self._num_qubits_per_reg, + self._num_regs, + self._bit_encoding + ) + + # Apply diagonal matrix to qubit register + quest.leftapplyFullStateDiagMatr(c_register, diag_matr) + quest.destroyFullStateDiagMatr(diag_matr) + + finally: + _current_phase_func = None + + return 0 cdef class QFT(MultiQubitOperator): @@ -971,7 +1168,7 @@ cdef class QFT(MultiQubitOperator): self.TYPE = OP_TYPES.OP_QFT cdef int apply_to(self, Qureg c_register) except -1: - quest.applyQFT(c_register, self._targets, self._num_targets) + quest.applyQuantumFourierTransform(c_register, self._targets, self._num_targets) cdef class FullQFT(GlobalOperator): @@ -980,4 +1177,4 @@ cdef class FullQFT(GlobalOperator): self.TYPE = OP_TYPES.OP_FULL_QFT cdef int apply_to(self, Qureg c_register) except -1: - quest.applyFullQFT(c_register) + quest.applyFullQuantumFourierTransform(c_register) diff --git a/pyquest/quest_interface.pxd b/pyquest/quest_interface.pxd index 2784972..1d71207 100644 --- a/pyquest/quest_interface.pxd +++ b/pyquest/quest_interface.pxd @@ -1,6 +1,6 @@ -"""C-level interface to QuEST functionality. +"""C-level interface to QuEST v4 functionality. -This pxd-only module provides the low-level interface to the QuEST API. +This pxd-only module provides the low-level interface to the QuEST v4 API. It also contains the type definitions for the qreal and qcomp data types and structs used by QuEST so they can be passed as arguments to the API functions. @@ -10,11 +10,10 @@ The QuEST native error handling is extended to throw Python exceptions Attributes: qreal: Floating point data type used by QuEST internals. - qcomp: Complex floating point data type equivalent to the QuEST - data type of the same name. Rarely used, the QuEST API uses - the ``Complex`` struct for complex numbers. + qcomp: Complex floating point data type used by QuEST. """ +# Precision types IF QuEST_PREC == 1: ctypedef float qreal ctypedef float complex qcomp @@ -25,83 +24,139 @@ ELIF QuEST_PREC == 4: ctypedef long double qreal ctypedef long double complex qcomp +ctypedef long long int qindex cdef extern from "quest_error.h": - # Place ``#include "error_handler.h"`` in the .cxx files. pass +# Vector structure (needed for rotation functions) +# Defined as pure Cython struct, not from quest.h (doesn't exist in v4) +ctypedef struct Vector: + qreal x + qreal y + qreal z -cdef extern from "QuEST.h": - # Data structures +cdef extern from "quest.h": + # Opaque environment structure (v4 API) ctypedef struct QuESTEnv: + int isMultithreaded + int isGpuAccelerated + int isDistributed + int isCuQuantumEnabled + int isGpuSharingEnabled int rank - int numRanks - int numSeeds - unsigned long int *seeds - ctypedef struct Complex: - qreal real - qreal imag - ctypedef struct ComplexArray: - qreal *real - qreal *imag - ctypedef struct ComplexMatrix2: - qreal real[2][2] - qreal imag[2][2] - ctypedef struct ComplexMatrix4: - qreal real[4][4] - qreal imag[4][4] - ctypedef struct ComplexMatrixN: - int numQubits - qreal **real - qreal **imag + int numNodes + + # Quantum register structure (v4 API) ctypedef struct Qureg: - int chunkId - ComplexArray stateVec - long long int numAmpsTotal - int numQubitsRepresented + int isMultithreaded + int isGpuAccelerated + int isDistributed + int rank + int numNodes + int logNumNodes int isDensityMatrix - ctypedef struct DiagonalOp: int numQubits - ctypedef struct PauliHamil: - pass - ctypedef struct Vector: - qreal x, y, z + qindex numAmps + qindex logNumAmps + qindex numAmpsPerNode + qindex logNumAmpsPerNode + +cdef extern from "quest.h": + # Pauli operator enumeration ctypedef enum pauliOpType: - pass - ctypedef enum bitEncoding: - UNSIGNED - TWOS_COMPLEMENT - ctypedef enum phaseFunc: - NORM - SCALED_NORM - INVERSE_NORM - SCALED_INVERSE_NORM - SCALED_INVERSE_SHIFTED_NORM - PRODUCT - SCALED_PRODUCT - INVERSE_PRODUCT - SCALED_INVERSE_PRODUCT - DISTANCE - SCALED_DISTANCE - INVERSE_DISTANCE - SCALED_INVERSE_DISTANCE - SCALED_INVERSE_SHIFTED_DISTANCE + PAULI_I=0 + PAULI_X=1 + PAULI_Y=2 + PAULI_Z=3 + +cdef extern from "quest.h": + # Fixed-size 1-qubit matrix (2x2) + ctypedef struct CompMatr1: + int numQubits + qindex numRows + qcomp elems[2][2] + + # Fixed-size 2-qubit matrix (4x4) + ctypedef struct CompMatr2: + int numQubits + qindex numRows + qcomp elems[4][4] + + # Variable-size N-qubit matrix + ctypedef struct CompMatr: + int numQubits + qindex numRows + int* isApproxUnitary + int* isApproxHermitian + int* wasGpuSynced + qcomp** cpuElems + qcomp* cpuElemsFlat + qcomp* gpuElemsFlat + + # Fixed-size 1-qubit diagonal matrix + ctypedef struct DiagMatr1: + int numQubits + qindex numElems + qcomp elems[2] + + # Fixed-size 2-qubit diagonal matrix + ctypedef struct DiagMatr2: + int numQubits + qindex numElems + qcomp elems[4] + + # Variable-size diagonal matrix + ctypedef struct DiagMatr: + int numQubits + qindex numElems + int* isApproxUnitary + int* isApproxHermitian + int* isApproxNonZero + int* isStrictlyNonNegative + int* wasGpuSynced + qcomp* cpuElems + qcomp* gpuElems + + # Full-state diagonal matrix + ctypedef struct FullStateDiagMatr: + int numQubits + qindex numElems + int isGpuAccelerated + int isMultithreaded + int isDistributed + qindex numElemsPerNode + int* isApproxUnitary + int* isApproxHermitian + int* isApproxNonZero + int* isStrictlyNonNegative + int* wasGpuSynced + qcomp* cpuElems + qcomp* gpuElems - # QuESTEnv methods - QuESTEnv createQuESTEnv() except + - void destroyQuESTEnv(QuESTEnv env) except + - void getEnvironmentString(QuESTEnv env, char str[200]) - void getQuESTSeeds(QuESTEnv env, unsigned long int **seeds, int *numSeeds) - void seedQuEST(QuESTEnv *env, unsigned long int *seedArray, int numSeeds) + # Environment functions (v4 API) + void initQuESTEnv() except + + QuESTEnv getQuESTEnv() except + + void finalizeQuESTEnv() except + + + # Environment queries (v4 API) + void getEnvironmentString(char str[200]) except + + + # Seeding functions (v4 API) + void setSeeds(unsigned* seeds, int numSeeds) except + + void setSeedsToDefault() except + + void getSeeds(unsigned* seeds) except + + int getNumSeeds() except + - # Quantum register methods - Qureg createQureg(int numQubits, QuESTEnv env) except + - Qureg createDensityQureg(int numQubits, QuESTEnv env) except + - void destroyQureg(Qureg qureg, QuESTEnv env) except + + # Qureg functions + Qureg createQureg(int numQubits) except + + Qureg createDensityQureg(int numQubits) except + + Qureg createCloneQureg(Qureg qureg) except + + void destroyQureg(Qureg qureg) except + + void cloneQureg(Qureg targetQureg, Qureg copyQureg) except + + void setQuregToClone(Qureg targetQureg, Qureg copyQureg) except + # State initializations - Qureg createCloneQureg(Qureg qureg, QuESTEnv env) except + - void cloneQureg(Qureg targetQureg, Qureg copyQureg) except + void initBlankState(Qureg qureg) except + void initClassicalState(Qureg qureg, long long int stateInd) except + void initPlusState(Qureg qureg) except + @@ -110,64 +165,56 @@ cdef extern from "QuEST.h": void initZeroState(Qureg qureg) except + void setAmps(Qureg qureg, long long int startInd, qreal* reals, qreal* imags, long long int numAmps) except + - void setWeightedQureg(Complex fac1, Qureg qureg1, Complex fac2, - Qureg qureg2, Complex facOut, Qureg out) except + + void setQuregAmps(Qureg qureg, long long int startIdx, qcomp* amps, long long int numAmps) except + + void setWeightedQureg(qcomp fac1, Qureg qureg1, qcomp fac2, + Qureg qureg2, qcomp facOut, Qureg out) except + + void setQuregToWeightedSum(Qureg out, qcomp* coeffs, Qureg* quregs, int numQuregs) except + - # Generic operators - ComplexMatrixN createComplexMatrixN(int numQubits) except + - void destroyComplexMatrixN(ComplexMatrixN m) except + - DiagonalOp createDiagonalOp(int numQubits, QuESTEnv env) except + - void initDiagonalOp(DiagonalOp op, qreal* real, qreal* imag) except + - void destroyDiagonalOp(DiagonalOp op, QuESTEnv env) except + - void applyDiagonalOp(Qureg qureg, DiagonalOp op) except + - void applyMatrix2(Qureg qureg, int targetQubit, ComplexMatrix2 u) except + - void applyMatrix4(Qureg qureg, int targetQubit1, int targetQubit2, - ComplexMatrix4 u) except + - void applyMatrixN(Qureg qureg, int* targs, int numTargs, - ComplexMatrixN u) except + + # Generic operators (v4 API) + CompMatr1 getCompMatr1(qcomp** in_) except + + CompMatr2 getCompMatr2(qcomp* in_) except + + CompMatr createCompMatr(int numQubits) except + + void destroyCompMatr(CompMatr m) except + + void setCompMatr(CompMatr out, qcomp** in_) except + + + # Diagonal operators (v4 API) + FullStateDiagMatr createFullStateDiagMatr(int numQubits) except + + void destroyFullStateDiagMatr(FullStateDiagMatr matr) except + + void setFullStateDiagMatr(FullStateDiagMatr out, qindex startIdx, qcomp* in_, qindex numElems) except + + void syncFullStateDiagMatr(FullStateDiagMatr matr) except + + void leftapplyFullStateDiagMatr(Qureg qureg, FullStateDiagMatr matrix) except + + + # Diagonal operators from functions (v4 API - new for PhaseFunc refactoring) + void setFullStateDiagMatrFromMultiVarFunc( + FullStateDiagMatr out, + qcomp (*func)(qindex*), + int* numQubitsPerVar, + int numVars, + int areSigned + ) except + + + void setDiagMatrFromMultiVarFunc( + DiagMatr out, + qcomp (*func)(qindex*), + int* numQubitsPerVar, + int numVars, + int areSigned + ) except + + # Matrix application (v4 API) + void leftapplyCompMatr1(Qureg qureg, int target, CompMatr1 matrix) except + + void leftapplyCompMatr2(Qureg qureg, int target1, int target2, CompMatr2 matrix) except + + void leftapplyCompMatr(Qureg qureg, int* targets, int numTargets, + CompMatr matrix) except + void applyMultiControlledMatrixN(Qureg qureg, int* ctrls, int numCtrls, - int* targs, int numTargs, ComplexMatrixN u) except + + int* targs, int numTargs, CompMatr u) except + qreal calcExpecPauliSum( Qureg qureg, pauliOpType* allPauliCodes, qreal* termCoeffs, int numSumTerms, Qureg workspace) except + - void applyPhaseFunc( - Qureg qureg, int* qubits, int numQubits, bitEncoding encoding, - qreal* coeffs, qreal* exponents, int numTerms) except + - void applyPhaseFuncOverrides( - Qureg qureg, int* qubits, int numQubits, bitEncoding encoding, - qreal* coeffs, qreal* exponents, int numTerms, long long int* overrideInds, - qreal* overridePhases, int numOverrides) except + - void applyMultiVarPhaseFunc( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, qreal* coeffs, qreal* exponents, - int* numTermsPerReg) except + - void applyMultiVarPhaseFuncOverrides( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, qreal* coeffs, qreal* exponents, - int* numTermsPerReg, long long int* overrideInds, - qreal* overridePhases, int numOverrides) except + - void applyNamedPhaseFunc( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, phaseFunc functionNameCode) except + - void applyNamedPhaseFuncOverrides( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, phaseFunc functionNameCode, - long long int* overrideInds, qreal* overridePhases, int numOverrides) except + - void applyParamNamedPhaseFunc( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, phaseFunc functionNameCode, - qreal* params, int numParams) except + - void applyParamNamedPhaseFuncOverrides( - Qureg qureg, int* qubits, int* numQubitsPerReg, int numRegs, - bitEncoding encoding, phaseFunc functionNameCode, - qreal* params, int numParams, long long int* overrideInds, - qreal* overridePhases, int numOverrides) except + void applyPauliSum(Qureg inQureg, pauliOpType* allPauliCodes, qreal* termCoeffs, int numSumTerms, Qureg outQureg) except + - void applyTrotterCircuit(Qureg qureg, PauliHamil hamil, qreal time, - int order, int reps) except + - void applyQFT(Qureg qureg, int* qubits, int numQubits) except + - void applyFullQFT(Qureg qureg) except + + void applyTrotterCircuit(Qureg qureg, int time, qreal period) except + + void applyQuantumFourierTransform(Qureg qureg, int* qubits, int numQubits) except + + void applyFullQuantumFourierTransform(Qureg qureg) except + # Gates (measurements) qreal collapseToOutcome(Qureg qureg, int measureQubit, @@ -185,96 +232,119 @@ cdef extern from "QuEST.h": void controlledPauliY( Qureg qureg, int controlQubit, int targetQubit) except + void pauliZ(Qureg qureg, int targetQubit) except + - void controlledPhaseFlip(Qureg qureg, int idQubit1, int idQubit2) except + - void multiControlledPhaseFlip( + void applyTwoQubitPhaseFlip(Qureg qureg, int idQubit1, int idQubit2) except + + void applyMultiQubitPhaseFlip( Qureg qureg, int *controlQubits, int numControlQubits) except + - void phaseShift(Qureg qureg, int targetQubit, qreal angle) except + - void controlledPhaseShift( + void applyPhaseShift(Qureg qureg, int targetQubit, qreal angle) except + + void applyTwoQubitPhaseShift( Qureg qureg, int idQubit1, int idQubit2, qreal angle) except + - void multiControlledPhaseShift( + void applyMultiQubitPhaseShift( Qureg qureg, int *controlQubits, int numControlQubits, qreal angle) except + - void rotateX(Qureg qureg, int rotQubit, qreal angle) except + - void controlledRotateX( + void applyRotateX(Qureg qureg, int rotQubit, qreal angle) except + + void applyControlledRotateX( Qureg qureg, int controlQubit, int targetQubit, qreal angle) except + - void rotateY(Qureg qureg, int rotQubit, qreal angle) except + - void controlledRotateY( + void applyRotateY(Qureg qureg, int rotQubit, qreal angle) except + + void applyControlledRotateY( Qureg qureg, int controlQubit, int targetQubit, qreal angle) except + - void rotateZ(Qureg qureg, int rotQubit, qreal angle) except + - void controlledRotateZ( + void applyRotateZ(Qureg qureg, int rotQubit, qreal angle) except + + void applyControlledRotateZ( Qureg qureg, int controlQubit, int targetQubit, qreal angle) except + - void rotateAroundAxis( - Qureg qureg, int rotQubit, qreal angle, Vector axis) except + - void controlledRotateAroundAxis( + void applyRotateAroundAxis( + Qureg qureg, int rotQubit, qreal angle, qreal vx, qreal vy, qreal vz) except + + void applyControlledRotateAroundAxis( Qureg qureg, int controlQubit, int targetQubit, qreal angle, - Vector axis) except + + qreal vx, qreal vy, qreal vz) except + void multiRotateZ( Qureg qureg, int* qubits, int numQubits, qreal angle) except + void multiRotatePauli( Qureg qureg, int* targetQubits, pauliOpType* targetPaulis, int numTargets, qreal angle) except + - void hadamard(Qureg qureg, int targetQubit) except + - void sGate(Qureg qureg, int targetQubit) except + - void tGate(Qureg qureg, int targetQubit) except + - void swapGate(Qureg qureg, int qubit1, int qubit2) except + - void sqrtSwapGate(Qureg qureg, int qb1, int qb2) except + + +cdef extern from "quest.h": + void applyHadamard(Qureg qureg, int targetQubit) except + + void applyS(Qureg qureg, int targetQubit) except + + void applyT(Qureg qureg, int targetQubit) except + + void applySwap(Qureg qureg, int qubit1, int qubit2) except + + void applySqrtSwap(Qureg qureg, int qb1, int qb2) except + void compactUnitary( - Qureg qureg, int targetQubit, Complex alpha, Complex beta) except + + Qureg qureg, int targetQubit, qcomp alpha, qcomp beta) except + void controlledCompactUnitary( Qureg qureg, int controlQubit, int targetQubit, - Complex alpha, Complex beta) except + - void unitary(Qureg qureg, int targetQubit, ComplexMatrix2 u) except + - void controlledUnitary( - Qureg qureg, int controlQubit, int targetQubit, ComplexMatrix2 u) except + - void multiControlledUnitary( - Qureg qureg, int* controlQubits, int numControlQubits, - int targetQubit, ComplexMatrix2 u) except + - void multiStateControlledUnitary( - Qureg qureg, int* controlQubits, int* controlState, - int numControlQubits, int targetQubit, ComplexMatrix2 u) except + - void twoQubitUnitary( - Qureg qureg, int targetQubit1, int targetQubit2, ComplexMatrix4 u) except + - void controlledTwoQubitUnitary( - Qureg qureg, int controlQubit, int targetQubit1, int targetQubit2, - ComplexMatrix4 u) except + - void multiControlledTwoQubitUnitary( - Qureg qureg, int* controlQubits, int numControlQubits, - int targetQubit1, int targetQubit2, ComplexMatrix4 u) except + - void multiQubitUnitary( - Qureg qureg, int* targs, int numTargs, ComplexMatrixN u) except + - void controlledMultiQubitUnitary( - Qureg qureg, int ctrl, int* targs, int numTargs, ComplexMatrixN u) except + - void multiControlledMultiQubitUnitary( - Qureg qureg, int* ctrls, int numCtrls, int* targs, - int numTargs, ComplexMatrixN u) except + + qcomp alpha, qcomp beta) except + + # Unitaries (v4 API) + void applyCompMatr1(Qureg qureg, int targetQubit, CompMatr1 u) except + + void applyControlledCompMatr1(Qureg qureg, int control, int target, CompMatr1 matrix) except + + void applyMultiControlledCompMatr1(Qureg qureg, int* controls, int numControls, int target, CompMatr1 matrix) except + + void applyMultiStateControlledCompMatr1(Qureg qureg, int* controls, int* states, int numControls, int target, CompMatr1 matrix) except + + void applyCompMatr2(Qureg qureg, int targetQubit1, int targetQubit2, + CompMatr2 u) except + + void applyControlledCompMatr2(Qureg qureg, int control, int target1, int target2, CompMatr2 matrix) except + + void applyMultiControlledCompMatr2(Qureg qureg, int* controls, int numControls, int target1, int target2, CompMatr2 matrix) except + + void applyMultiStateControlledCompMatr2(Qureg qureg, int* controls, int* states, int numControls, int target1, int target2, CompMatr2 matrix) except + + void applyCompMatr(Qureg qureg, int* targets, int numTargets, + CompMatr u) except + + void applyControlledCompMatr(Qureg qureg, int control, int* targets, int numTargets, CompMatr matrix) except + + void applyMultiControlledCompMatr(Qureg qureg, int* controls, int numControls, int* targets, int numTargets, CompMatr matrix) except + + void applyMultiStateControlledCompMatr(Qureg qureg, int* controls, int* states, int numControls, int* targets, int numTargets, CompMatr matrix) except + + + # Pauli Gates (v4 API) + void applyPauliX(Qureg qureg, int target) except + + void applyPauliY(Qureg qureg, int target) except + + void applyPauliZ(Qureg qureg, int target) except + + void applyControlledPauliX(Qureg qureg, int control, int target) except + + void applyControlledPauliY(Qureg qureg, int control, int target) except + + void applyControlledPauliZ(Qureg qureg, int control, int target) except + + void applyMultiControlledPauliX(Qureg qureg, int* controls, int numControls, int target) except + + void applyMultiControlledPauliY(Qureg qureg, int* controls, int numControls, int target) except + + void applyMultiControlledPauliZ(Qureg qureg, int* controls, int numControls, int target) except + + + # NOT Gates (v4 API) + void applyMultiQubitNot(Qureg qureg, int* targets, int numTargets) except + + void applyControlledMultiQubitNot(Qureg qureg, int control, int* targets, int numTargets) except + + void applyMultiControlledMultiQubitNot(Qureg qureg, int* controls, int numControls, int* targets, int numTargets) except + + + # Phase Flip Gates (v4 API) + void applyPhaseFlip(Qureg qureg, int target) except + + void applyTwoQubitPhaseFlip(Qureg qureg, int control, int target) except + + void applyMultiQubitPhaseFlip(Qureg qureg, int* qubits, int numQubits) except + + + # Measurement Functions (v4 API) + int applyQubitMeasurementAndGetProb(Qureg qureg, int qubit, qreal *outcomeProb) except + + int applyForcedQubitMeasurement(Qureg qureg, int qubit, int outcome) except + + void mixDamping(Qureg qureg, int targetQubit, qreal prob) except + void mixDensityMatrix(Qureg combineQureg, qreal prob, Qureg otherQureg) except + + void mixQureg(Qureg qureg, Qureg otherQureg, qreal prob) except + void mixDephasing(Qureg qureg, int targetQubit, qreal prob) except + void mixDepolarising(Qureg qureg, int targetQubit, qreal prob) except + - void mixKrausMap(Qureg qureg, int target, ComplexMatrix2 *ops, int numOps) except + + void mixKrausMap(Qureg qureg, int target, CompMatr1 *ops, int numOps) except + void mixMultiQubitKrausMap( Qureg qureg, int* targets, int numTargets, - ComplexMatrixN* ops, int numOps) except + + CompMatr* ops, int numOps) except + void mixPauli(Qureg qureg, int targetQubit, qreal probX, qreal probY, qreal probZ) except + + void mixPaulis(Qureg qureg, int targetQubit, qreal probX, qreal probY, qreal probZ) except + void mixTwoQubitDephasing(Qureg qureg, int qubit1, int qubit2, qreal prob) except + void mixTwoQubitDepolarising(Qureg qureg, int qubit1, int qubit2, qreal prob) except + void mixTwoQubitKrausMap( - Qureg qureg, int target1, int target2, ComplexMatrix4 *ops, int numOps) except + + Qureg qureg, int target1, int target2, CompMatr2 *ops, int numOps) except + # Calculations int getNumQubits(Qureg qureg) except + long long int getNumAmps(Qureg qureg) except + qreal calcPurity(Qureg qureg) except + qreal calcTotalProb(Qureg qureg) except + - Complex calcInnerProduct(Qureg bra, Qureg ket) except + + qcomp calcInnerProduct(Qureg bra, Qureg ket) except + qreal calcDensityInnerProduct(Qureg rho1, Qureg rho2) except + qreal calcFidelity(Qureg qureg, Qureg pureState) except + - Complex getAmp(Qureg qureg, long long int index) except + - Complex getDensityAmp(Qureg qureg, long long int row, long long int col) except + - void setAmps(Qureg qureg, long long int startInd, - qreal* reals, qreal* imags, long long int numAmps) except + + qcomp getQuregAmp(Qureg qureg, long long int index) except + + qcomp getDensityQuregAmp(Qureg qureg, long long int row, long long int col) except + + void setQuregAmps(Qureg qureg, long long int startIdx, qcomp* amps, long long int numAmps) except + + void setQuregToClone(Qureg targetQureg, Qureg copyQureg) except + qreal calcProbOfOutcome(Qureg qureg, int measureQubit, int outcome) except + void calcProbOfAllOutcomes(qreal* outcomeProbs, Qureg qureg, int* qubits, int numQubits) except + + void calcProbsOfAllMultiQubitOutcomes(qreal* outcomeProbs, Qureg qureg, + int* qubits, int numQubits) except + cdef enum OP_TYPES: diff --git a/pyquest/unitaries.pxd b/pyquest/unitaries.pxd index 6ee2113..4f434dd 100644 --- a/pyquest/unitaries.pxd +++ b/pyquest/unitaries.pxd @@ -1,8 +1,7 @@ from libc.stdlib cimport malloc, free cimport pyquest.quest_interface as quest -from pyquest.quest_interface cimport Complex, Vector, qreal, OP_TYPES, Qureg -from pyquest.quest_interface cimport ComplexMatrix2, ComplexMatrix4, ComplexMatrixN -from pyquest.quest_interface cimport createComplexMatrixN, destroyComplexMatrixN +from pyquest.quest_interface cimport qcomp, Vector, qreal, OP_TYPES, Qureg +from pyquest.quest_interface cimport CompMatr1, CompMatr2, CompMatr from pyquest.operators cimport SingleQubitOperator, MultiQubitOperator, MatrixOperator @@ -11,8 +10,8 @@ cdef class U(MatrixOperator): cdef class CompactU(SingleQubitOperator): - cdef Complex _alpha - cdef Complex _beta + cdef qcomp _alpha + cdef qcomp _beta cdef class X(MultiQubitOperator): @@ -75,4 +74,4 @@ cdef class MultiRotatePauli(BaseRotate): pass cdef class ParamSwap(BaseRotate): - cdef ComplexMatrix4 _u + cdef CompMatr2 _u diff --git a/pyquest/unitaries.pyx b/pyquest/unitaries.pyx index 432208d..32f3e35 100644 --- a/pyquest/unitaries.pyx +++ b/pyquest/unitaries.pyx @@ -84,15 +84,7 @@ cdef class U(MatrixOperator): # deallocation must also be performed separately. if self._num_targets > 2: if self._matrix != NULL: - destroyComplexMatrixN((self._matrix)[0]) - else: - free(self._real) - free(self._imag) - # The parent destructor might also call free on - # self._matrix, self._real, and self._imag. Setting them to - # NULL avoids freeing the same pointer twice. - self._real = NULL - self._imag = NULL + quest.destroyCompMatr((self._matrix)[0]) free(self._matrix) self._matrix = NULL free(self._control_pattern) @@ -117,84 +109,71 @@ cdef class U(MatrixOperator): cdef _create_array_property(self): # We need to override this method, because core QuEST - # supports ComplexMatrix* with controls for unitaries, but + # supports CompMatr1/CompMatr2 with controls for unitaries, but # not for generic matrices. cdef size_t matrix_dim = 1 # Assigning a 1 first prevents integer overflows of the bit shift matrix_dim = matrix_dim << self._num_targets - cdef size_t k if self._num_targets == 1: - self._matrix = malloc(sizeof(ComplexMatrix2)) - self._real = malloc(matrix_dim * sizeof(self._real[0])) - self._imag = malloc(matrix_dim * sizeof(self._imag[0])) - for k in range(matrix_dim): - self._real[k] = &(self._matrix).real[k][0] - self._imag[k] = &(self._matrix).imag[k][0] + self._matrix = malloc(sizeof(quest.CompMatr1)) elif self._num_targets == 2: - self._matrix = malloc(sizeof(ComplexMatrix4)) - self._real = malloc(matrix_dim * sizeof(self._real[0])) - self._imag = malloc(matrix_dim * sizeof(self._imag[0])) - for k in range(matrix_dim): - self._real[k] = &(self._matrix).real[k][0] - self._imag[k] = &(self._matrix).imag[k][0] + self._matrix = malloc(sizeof(quest.CompMatr2)) else: - self._matrix = malloc(sizeof(ComplexMatrixN)) - (self._matrix)[0] = createComplexMatrixN(self._num_targets) - self._real = (self._matrix).real - self._imag = (self._matrix).imag + self._matrix = malloc(sizeof(quest.CompMatr)) + (self._matrix)[0] = quest.createCompMatr(self._num_targets) cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: if self._num_targets == 1: - quest.unitary( + quest.applyCompMatr1( c_register, self._targets[0], - (self._matrix)[0]) + (self._matrix)[0]) elif self._num_targets == 2: - quest.twoQubitUnitary( + quest.applyCompMatr2( c_register, self._targets[0], self._targets[1], - (self._matrix)[0]) + (self._matrix)[0]) else: - quest.multiQubitUnitary( + quest.applyCompMatr( c_register, self._targets, self._num_targets, - (self._matrix)[0]) + (self._matrix)[0]) elif self._num_controls == 1: if self._num_targets == 1: if self._control_pattern == NULL: - quest.controlledUnitary( + quest.applyControlledCompMatr1( c_register, self._controls[0], self._targets[0], - (self._matrix)[0]) + (self._matrix)[0]) else: - quest.multiStateControlledUnitary( + quest.applyMultiStateControlledCompMatr1( c_register, self._controls, self._control_pattern, self._num_controls, self._targets[0], - (self._matrix)[0]) + (self._matrix)[0]) elif self._num_targets == 2: - quest.controlledTwoQubitUnitary( + quest.applyControlledCompMatr2( c_register, self._controls[0], self._targets[0], - self._targets[1], (self._matrix)[0]) + self._targets[1], (self._matrix)[0]) else: - quest.controlledMultiQubitUnitary( + quest.applyControlledCompMatr( c_register, self._controls[0], self._targets, - self._num_targets, (self._matrix)[0]) + self._num_targets, (self._matrix)[0]) else: if self._num_targets == 1: if self._control_pattern == NULL: - quest.multiControlledUnitary( + quest.applyMultiControlledCompMatr1( c_register, self._controls, self._num_controls, - self._targets[0], (self._matrix)[0]) + self._targets[0], (self._matrix)[0]) else: - quest.multiStateControlledUnitary( + quest.applyMultiStateControlledCompMatr1( c_register, self._controls, self._control_pattern, self._num_controls, self._targets[0], - (self._matrix)[0]) + (self._matrix)[0]) elif self._num_targets == 2: - quest.multiControlledTwoQubitUnitary( + quest.applyMultiControlledCompMatr2( c_register, self._controls, self._num_controls, self._targets[0], self._targets[1], - (self._matrix)[0]) + (self._matrix)[0]) else: - quest.multiControlledMultiQubitUnitary( + quest.applyMultiControlledCompMatr( c_register, self._controls, self._num_controls, self._targets, - self._num_targets, (self._matrix)[0]) + self._num_targets, (self._matrix)[0]) cdef class CompactU(SingleQubitOperator): @@ -210,13 +189,31 @@ cdef class CompactU(SingleQubitOperator): "Multi-controlled compact unitary not yet supported.") cdef int apply_to(self, Qureg c_register) except -1: + cdef CompMatr1 matrix + cdef qcomp[2][2] elems + cdef qcomp beta_conj + cdef qcomp alpha_conj + cdef qcomp *elem_ptrs[2] + # Compute conjugates manually + beta_conj.real = self._beta.real + beta_conj.imag = -self._beta.imag + alpha_conj.real = self._alpha.real + alpha_conj.imag = -self._alpha.imag + # Build matrix + elems[0][0] = self._alpha + elems[0][1] = -beta_conj + elems[1][0] = self._beta + elems[1][1] = alpha_conj + # Set up row pointers for getCompMatr1 + elem_ptrs[0] = &(elems[0][0]) + elem_ptrs[1] = &(elems[1][0]) + matrix = quest.getCompMatr1(elem_ptrs) + if self._num_controls == 0: - quest.compactUnitary( - c_register, self._target, self._alpha, self._beta) + quest.applyCompMatr1(c_register, self._target, matrix) if self._num_controls == 1: - quest.controlledCompactUnitary( - c_register, self._controls[0], self._target, - self._alpha, self._beta) + quest.applyControlledCompMatr1( + c_register, self._controls[0], self._target, matrix) cdef class X(MultiQubitOperator): @@ -227,15 +224,15 @@ cdef class X(MultiQubitOperator): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: if self._num_targets == 1: - quest.pauliX(c_register, self._targets[0]) + quest.applyPauliX(c_register, self._targets[0]) else: - quest.multiQubitNot(c_register, self._targets, + quest.applyMultiQubitNot(c_register, self._targets, self._num_targets) elif self._num_controls == 1 and self._num_targets == 1: - quest.controlledNot( + quest.applyControlledPauliX( c_register, self._controls[0], self._targets[0]) else: - quest.multiControlledMultiQubitNot( + quest.applyMultiControlledMultiQubitNot( c_register, self._controls, self._num_controls, self._targets, self._num_targets) @@ -254,9 +251,9 @@ cdef class Y(SingleQubitOperator): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: - quest.pauliY(c_register, self._target) + quest.applyPauliY(c_register, self._target) else: - quest.controlledPauliY( + quest.applyControlledPauliY( c_register, self._controls[0], self._target) @property @@ -275,20 +272,20 @@ cdef class Z(SingleQubitOperator): cdef int *controls cdef int m if self._num_controls == 0: - quest.pauliZ(c_register, self._target) + quest.applyPauliZ(c_register, self._target) elif self._num_controls == 1: - quest.controlledPhaseFlip( + quest.applyTwoQubitPhaseFlip( c_register, self._controls[0], self._target) else: # We need a single array containing all qubits the # multi-controlled Pauli-Z is acting on, because of - # the signature of multiControlledPhaseFlip(...). + # the signature of applyMultiQubitPhaseFlip(...). controls = malloc((self._num_controls + 1) * sizeof(controls[0])) controls[0] = self._target for m in range(self._num_controls): controls[m + 1] = self._controls[m] - quest.multiControlledPhaseFlip( + quest.applyMultiQubitPhaseFlip( c_register, controls, self._num_controls + 1) free(controls) @@ -305,7 +302,7 @@ cdef class Swap(MultiQubitOperator): raise ValueError("Swap gate must act on exactly two qubits") cdef int apply_to(self, Qureg c_register) except -1: - quest.swapGate(c_register, self._targets[0], self._targets[1]) + quest.applySwap(c_register, self._targets[0], self._targets[1]) @property def inverse(self): @@ -320,7 +317,7 @@ cdef class SqrtSwap(MultiQubitOperator): raise ValueError("Sqrt-swap gate must act on exactly two qubits.") cdef int apply_to(self, Qureg c_register) except -1: - quest.sqrtSwapGate(c_register, self._targets[0], self._targets[1]) + quest.applySqrtSwap(c_register, self._targets[0], self._targets[1]) cdef class H(SingleQubitOperator): @@ -329,7 +326,7 @@ cdef class H(SingleQubitOperator): self.TYPE = OP_TYPES.OP_HADAMARD cdef int apply_to(self, Qureg c_register) except -1: - quest.hadamard(c_register, self._target) + quest.applyHadamard(c_register, self._target) @property def inverse(self): @@ -342,7 +339,7 @@ cdef class S(SingleQubitOperator): self.TYPE = OP_TYPES.OP_S cdef int apply_to(self, Qureg c_register) except -1: - quest.sGate(c_register, self._target) + quest.applyS(c_register, self._target) cdef class T(SingleQubitOperator): @@ -351,7 +348,7 @@ cdef class T(SingleQubitOperator): self.TYPE = OP_TYPES.OP_T cdef int apply_to(self, Qureg c_register) except -1: - quest.tGate(c_register, self._target) + quest.applyT(c_register, self._target) cdef class BaseRotate(SingleQubitOperator): @@ -404,9 +401,9 @@ cdef class Rx(BaseRotate): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: - quest.rotateX(c_register, self._target, self._angle) + quest.applyRotateX(c_register, self._target, self._angle) else: - quest.controlledRotateX( + quest.applyControlledRotateX( c_register, self._controls[0], self._target, self._angle) @@ -420,9 +417,9 @@ cdef class Ry(BaseRotate): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: - quest.rotateY(c_register, self._target, self._angle) + quest.applyRotateY(c_register, self._target, self._angle) else: - quest.controlledRotateY( + quest.applyControlledRotateY( c_register, self._controls[0], self._target, self._angle) @@ -436,9 +433,9 @@ cdef class Rz(BaseRotate): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: - quest.rotateZ(c_register, self._target, self._angle) + quest.applyRotateZ(c_register, self._target, self._angle) else: - quest.controlledRotateZ( + quest.applyControlledRotateZ( c_register, self._controls[0], self._target, self._angle) @@ -451,9 +448,9 @@ cdef class Phase(BaseRotate): cdef size_t m cdef int *controls if self._num_controls == 0: - quest.phaseShift(c_register, self._target, self._angle) + quest.applyPhaseShift(c_register, self._target, self._angle) elif self._num_controls == 1: - quest.controlledPhaseShift( + quest.applyTwoQubitPhaseShift( c_register, self._controls[0], self._target, self._angle) else: # The same as for multi-controlled PauliZ applies here. @@ -462,7 +459,7 @@ cdef class Phase(BaseRotate): controls[0] = self._target for m in range(self._num_controls): controls[m + 1] = self._controls[m] - quest.multiControlledPhaseShift( + quest.applyMultiQubitPhaseShift( c_register, controls, self._num_controls + 1, self._angle) free(controls) @@ -480,12 +477,13 @@ cdef class RotateAroundAxis(BaseRotate): cdef int apply_to(self, Qureg c_register) except -1: if self._num_controls == 0: - quest.rotateAroundAxis( - c_register, self._target, self._angle, self._axis) + quest.applyRotateAroundAxis( + c_register, self._target, self._angle, + self._axis.x, self._axis.y, self._axis.z) else: - quest.controlledRotateAroundAxis( + quest.applyControlledRotateAroundAxis( c_register, self._controls[0], self._target, - self._angle, self._axis) + self._angle, self._axis.x, self._axis.y, self._axis.z) cdef class MultiRotatePauli(BaseRotate): From 0cc5c17f5530eea4aff96ec0ac9e853307ae5a0f Mon Sep 17 00:00:00 2001 From: Maurice Jamieson Date: Mon, 8 Dec 2025 17:31:02 +0000 Subject: [PATCH 2/2] Updates for PR 17 comments --- pyquest/core.pyx | 90 +++++++++++++------------------------------ pyquest/gates.pyx | 4 +- pyquest/operators.pyx | 43 +++++++++++++++------ pyquest/unitaries.pyx | 22 +++++++++++ 4 files changed, 83 insertions(+), 76 deletions(-) diff --git a/pyquest/core.pyx b/pyquest/core.pyx index 50ed6df..d883dc2 100644 --- a/pyquest/core.pyx +++ b/pyquest/core.pyx @@ -161,10 +161,14 @@ cdef class QuESTEnvironment: @seeds.setter def seeds(self, value): - cdef unsigned[:] seed_array = np.asarray(value, dtype=np.uint32) - if seed_array.size != quest.getNumSeeds(): + seed_values = np.asarray(value) + + # Validate values are non-negative + if np.any(seed_values < 0): raise ValueError( - f"Expected {quest.getNumSeeds()} seeds, got {seed_array.size}") + f"Seed values must be non-negative {value}.") + + cdef unsigned[:] seed_array = np.ascontiguousarray(seed_values, dtype=np.uint32) quest.setSeeds(&seed_array[0], seed_array.size) def reset_seeds_to_default(self): @@ -348,10 +352,6 @@ cdef class Register: quest.setQuregToWeightedSum(res_reg.c_register, &(coeffs[0]), &(quregs[0]), 2) - # Result has scaling factor 1 since it's the sum - res_reg._scaling_factor.real = 1.0 - res_reg._scaling_factor.imag = 0.0 - return res_reg def __sub__(left, right): @@ -362,20 +362,15 @@ cdef class Register: cdef Register res_reg = Register.zero_like(left_reg) cdef qcomp[2] coeffs cdef Qureg[2] quregs - + # Create result = 1.0 * left_scaling * left - 1.0 * right_scaling * right coeffs[0] = left_reg._scaling_factor - coeffs[1].real = -right_reg._scaling_factor.real - coeffs[1].imag = -right_reg._scaling_factor.imag + coeffs[1] = -right_reg._scaling_factor quregs[0] = left_reg.c_register quregs[1] = right_reg.c_register quest.setQuregToWeightedSum(res_reg.c_register, &(coeffs[0]), &(quregs[0]), 2) - # Result has scaling factor 1 since it's the difference - res_reg._scaling_factor.real = 1.0 - res_reg._scaling_factor.imag = 0.0 - return res_reg def __getitem__(self, index): @@ -471,8 +466,9 @@ cdef class Register: cdef int step cdef qreal val_imag, val_real cdef bool_t from_scalar - cdef const qcomp[:] value_arr + cdef qcomp[:] value_arr cdef qcomp amp + cdef qcomp[1] amp_arr try: value[0] from_scalar = False @@ -493,49 +489,34 @@ cdef class Register: try: value_arr = value for m in range(num_index): - val_real = value_arr[m].real - val_imag = value_arr[m].imag - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, k, &, 1) + quest.setQuregAmps(self.c_register, k, &value_arr[m], 1) k += step except (TypeError, ValueError): for m in range(num_index): - val_real = value[m].real - val_imag = value[m].imag - # Because QuEST needs separate arrays for real - # and imaginary parts, we cannot hand off a - # pointer to a single numpy memoryview. Thus we - # call setQuregAmps for each element separately. - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, k, &, 1) + quest.setQuregAmps(self.c_register, k, &value_arr[m], 1) k += step else: - val_real = value.real - val_imag = value.imag + amp = value + amp_arr[0] = amp for m in range(num_index): - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, k, &, 1) + quest.setQuregAmps(self.c_register, k, &_arr[0], 1) k += step else: try: num_index = len(index) if not from_scalar: + value_arr = value for m in range(num_index): - val_real = value[m].real - val_imag = value[m].imag - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, index[m], &, 1) + quest.setQuregAmps(self.c_register, index[m], &value_arr[m], 1) else: - val_real = value.real - val_imag = value.imag + amp = value + amp_arr[0] = amp for m in range(num_index): - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, index[m], &, 1) + quest.setQuregAmps(self.c_register, index[m], &_arr[0], 1) except TypeError: # Last guess is we got a scalar index. - val_real = value.real - val_imag = value.imag - amp = val_real + val_imag * 1j - quest.setQuregAmps(self.c_register, index, &, 1) + amp = value + amp_arr[0] = amp + quest.setQuregAmps(self.c_register, index, &_arr[0], 1) @property def is_alive(self): @@ -678,13 +659,8 @@ cdef class Register: """ self._apply_delayed_operations() other._apply_delayed_operations() - cdef qcomp prod - if self.c_register.isDensityMatrix: - return quest.calcInnerProduct( - self.c_register, other.c_register) - else: - prod = quest.calcInnerProduct(self.c_register, other.c_register) - return prod.real + 1j * prod.imag + return quest.calcInnerProduct(self.c_register, other.c_register) + cpdef qreal fidelity(self, Register other): """Calculate fidelity with the pure state in another register. @@ -833,21 +809,9 @@ cdef class Register: self._apply_scaling() cdef void _apply_scaling(self): - cdef qcomp zero - cdef qcomp[3] coeffs - cdef Qureg[3] quregs if self._scaling_factor.real != 1 or self._scaling_factor.imag != 0: self._ensure_no_borrow() - zero.real = 0. - zero.imag = 0. - # Use setQuregToWeightedSum instead of setWeightedQureg - coeffs[0] = self._scaling_factor - coeffs[1] = zero - coeffs[2] = zero - quregs[0] = self.c_register - quregs[1] = self.c_register - quregs[2] = self.c_register - quest.setQuregToWeightedSum(self.c_register, &(coeffs[0]), &(quregs[0]), 3) + quest.setQuregToWeightedSum(self.c_register, &self._scaling_factor, &self.c_register, 1) self._scaling_factor.real = 1 self._scaling_factor.imag = 0 diff --git a/pyquest/gates.pyx b/pyquest/gates.pyx index e93e27e..a8d0fbe 100644 --- a/pyquest/gates.pyx +++ b/pyquest/gates.pyx @@ -100,9 +100,9 @@ cdef class M(MultiQubitOperator): self._results[k] = quest.applyQubitMeasurementAndGetProb( c_register, self._targets[k], &(self._probabilities[k])) else: - self._probabilities[k] = 1.0 - self._results[k] = quest.applyForcedQubitMeasurement( + self._probabilities[k] = quest.applyForcedQubitMeasurement( c_register, self._targets[k], self._force[k]) + self._results[k] = self._force[k] @property def results(self): diff --git a/pyquest/operators.pyx b/pyquest/operators.pyx index 04e78b4..0ce8892 100644 --- a/pyquest/operators.pyx +++ b/pyquest/operators.pyx @@ -80,15 +80,13 @@ cdef class BaseOperator: cdef Qureg tmp_reg = quest.createQureg(num_qubits) cdef long long k, m cdef long long mat_dim = 1LL << num_qubits - cdef qcomp amp cdef qcomp[:, :] res_mat = np.ndarray( (mat_dim, mat_dim), dtype=pyquest.core.np_qcomp) for k in range(mat_dim): quest.initClassicalState(tmp_reg, k) self.apply_to(tmp_reg) for m in range(mat_dim): - amp = quest.getQuregAmp(tmp_reg, m) - res_mat[m, k] = amp + res_mat[m, k] = quest.getQuregAmp(tmp_reg, m) quest.destroyQureg(tmp_reg) return res_mat.base @@ -347,12 +345,34 @@ cdef class MatrixOperator(MultiQubitOperator): cdef _create_array_property(self): cdef size_t matrix_dim = 1 # Assigning a 1 first prevents integer overflows of the bit shift matrix_dim = matrix_dim << self._num_targets + cdef size_t k + cdef quest.CompMatr1* matr1 + cdef quest.CompMatr2* matr2 + cdef qcomp* row_ptr # The only cases where core QuEST supports CompMatr1/CompMatr2 # for generic matrices are non-controlled cases. if self._num_targets == 1 and self._num_controls == 0: self._matrix = malloc(sizeof(quest.CompMatr1)) + matr1 = self._matrix + matr1.numQubits = self._num_targets + matr1.numRows = 1 << self._num_targets + self._real = malloc(matrix_dim * sizeof(self._real[0])) + self._imag = malloc(matrix_dim * sizeof(self._imag[0])) + for k in range(matrix_dim): + row_ptr = &(matr1.elems[k][0]) + self._real[k] = row_ptr + self._imag[k] = row_ptr + 1 elif self._num_targets == 2 and self._num_controls == 0: self._matrix = malloc(sizeof(quest.CompMatr2)) + matr2 = self._matrix + matr2.numQubits = self._num_targets + matr2.numRows = 1 << self._num_targets + self._real = malloc(matrix_dim * sizeof(self._real[0])) + self._imag = malloc(matrix_dim * sizeof(self._imag[0])) + for k in range(matrix_dim): + row_ptr = &(matr2.elems[k][0]) + self._real[k] = row_ptr + self._imag[k] = row_ptr + 1 else: self._matrix = malloc(sizeof(quest.CompMatr)) (self._matrix)[0] = quest.createCompMatr(self._num_targets) @@ -401,7 +421,8 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) + # Direct cast to qcomp handles complex conversion + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -409,7 +430,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m] + 0j)) + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -417,7 +438,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m] + 0j)) + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -425,7 +446,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m] + 0j)) + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -433,7 +454,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -441,7 +462,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) + self._set_matrix_element(k, m, arr[k, m]) @cython.boundscheck(False) @cython.wraparound(False) @@ -449,7 +470,7 @@ cdef class MatrixOperator(MultiQubitOperator): cdef size_t k, m for k in range(arr.shape[0]): for m in range(arr.shape[1]): - self._set_matrix_element(k, m, (arr[k, m].real + 1j * arr[k, m].imag)) + self._set_matrix_element(k, m, arr[k, m]) cdef class PauliSum(GlobalOperator): @@ -486,7 +507,7 @@ cdef class DiagonalOperator(GlobalOperator): "2**num_qubits") # Convert Python complex array to qcomp array num_elems = diag_elements.size - qcomp_el = np.ascontiguousarray(diag_elements, dtype=np.complex128) + qcomp_el = np.ascontiguousarray(diag_elements, dtype=pyquest.core.np_qcomp) quest.setFullStateDiagMatr(self._diag_op, 0, &qcomp_el[0], num_elems) def __dealloc__(self): diff --git a/pyquest/unitaries.pyx b/pyquest/unitaries.pyx index 32f3e35..de0d202 100644 --- a/pyquest/unitaries.pyx +++ b/pyquest/unitaries.pyx @@ -113,10 +113,32 @@ cdef class U(MatrixOperator): # not for generic matrices. cdef size_t matrix_dim = 1 # Assigning a 1 first prevents integer overflows of the bit shift matrix_dim = matrix_dim << self._num_targets + cdef size_t k + cdef quest.CompMatr1* matr1 + cdef quest.CompMatr2* matr2 + cdef qcomp* row_ptr if self._num_targets == 1: self._matrix = malloc(sizeof(quest.CompMatr1)) + matr1 = self._matrix + matr1.numQubits = self._num_targets + matr1.numRows = matrix_dim + self._real = malloc(matrix_dim * sizeof(self._real[0])) + self._imag = malloc(matrix_dim * sizeof(self._imag[0])) + for k in range(matrix_dim): + row_ptr = &(matr1.elems[k][0]) + self._real[k] = row_ptr + self._imag[k] = row_ptr + 1 elif self._num_targets == 2: self._matrix = malloc(sizeof(quest.CompMatr2)) + matr2 = self._matrix + matr2.numQubits = self._num_targets + matr2.numRows = matrix_dim + self._real = malloc(matrix_dim * sizeof(self._real[0])) + self._imag = malloc(matrix_dim * sizeof(self._imag[0])) + for k in range(matrix_dim): + row_ptr = &(matr2.elems[k][0]) + self._real[k] = row_ptr + self._imag[k] = row_ptr + 1 else: self._matrix = malloc(sizeof(quest.CompMatr)) (self._matrix)[0] = quest.createCompMatr(self._num_targets)