diff --git a/hafnian/__init__.py b/hafnian/__init__.py index 1510b969c..c1ffa4b68 100644 --- a/hafnian/__init__.py +++ b/hafnian/__init__.py @@ -54,11 +54,14 @@ .. autosummary:: hafnian hafnian_repeated + hafnian_batched + hermite_multidimensional tor perm permanent_repeated reduction version + """ # pylint: disable=wrong-import-position import os @@ -66,13 +69,6 @@ import numpy as np - -if platform.system() == "Windows": # pragma: no cover - extra_dll_dir = os.path.join(os.path.dirname(__file__), ".libs") - if os.path.isdir(extra_dll_dir): - os.environ["PATH"] += os.pathsep + extra_dll_dir - - from ._hafnian import ( haf_complex, haf_int, @@ -83,10 +79,16 @@ hafnian_repeated, reduction, ) +from ._hermite_multidimensional import hafnian_batched, hermite_multidimensional from ._permanent import perm, perm_complex, perm_real, permanent_repeated from ._torontonian import tor from ._version import __version__ +if platform.system() == "Windows": # pragma: no cover + extra_dll_dir = os.path.join(os.path.dirname(__file__), ".libs") + if os.path.isdir(extra_dll_dir): + os.environ["PATH"] += os.pathsep + extra_dll_dir + __all__ = [ "hafnian", @@ -95,6 +97,8 @@ "perm", "permanent_repeated", "reduction", + "hafnian_batched", + "hermite_multidimensional", "version", ] diff --git a/hafnian/_hafnian.py b/hafnian/_hafnian.py index 435ced030..f3428a4d5 100644 --- a/hafnian/_hafnian.py +++ b/hafnian/_hafnian.py @@ -174,7 +174,10 @@ def hafnian( def hafnian_repeated(A, rpt, mu=None, loop=False, tol=1e-12): - r"""Returns the hafnian of matrix with repeated rows/columns. + r""" + THIS NEEDS TO BE WRITTEN + + Returns the hafnian of matrix with repeated rows/columns. The :func:`reduction` function may be used to show the resulting matrix with repeated rows and columns as per ``rpt``. diff --git a/hafnian/_hermite_multidimensional.py b/hafnian/_hermite_multidimensional.py new file mode 100644 index 000000000..425ae9f8b --- /dev/null +++ b/hafnian/_hermite_multidimensional.py @@ -0,0 +1,145 @@ +# Copyright 2019 Xanadu Quantum Technologies Inc. + +# 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. +""" +Hermite Multidimensional Python interface +""" +from itertools import product +import numpy as np + +from .lib.libhaf import hermite_multidimensional as hm +from ._hafnian import input_validation + + +def return_prod(C, index): + r"""Given an array :math:`C_{i,j}` and an array or list of indices :math:`index = [i_1,i_2,i_3,\dots,i_n] `, returns :math:`prod_{k=1}^n C_{1,i_1}`. + + Args: + C (array): An array + index (array): A set of indices + + Returns: + complex: The product of the array elements determines by index + """ + return np.prod([C[mode, val] for mode, val in enumerate(index)]) + + +def expansion_coeff(alpha, resolution, renorm=True): + r"""Returns the (quasi) geometric series as a vector with components :math:`alpha^i/sqrt(i!)` for :math:`0 \leq i < resolution`. + + If renorm is false it omits the division by factorials + + Args: + alpha (complex): Ratio of the geometric series + resoluton (int): Cutoff of the geometric series + renorm (bool): Decides whether to normalize by the factorials + + Returns: + array: The power series + """ + vals = np.empty([resolution], dtype=type(alpha)) + vals[0] = 1.0 + if renorm: + for i in range(1, resolution): + vals[i] = vals[i - 1] * alpha / np.sqrt(i) + else: + for i in range(1, resolution): + vals[i] = vals[i - 1] * alpha + return vals + + +def hermite_multidimensional(R, resolution, y=None, renorm=False, make_tensor=True): + r"""Returns the multidimensional Hermite polynomials :math:`H_k^{(R)}(y)`. + + Here :math:`R` is an :math:n \times n: square matrix, + :math:`y` is an :math:`n` dimensional vector. The polynomials are parametrized by the multi-index :math:`k=(k_0,k_1,\ldots,k_{n-1}) + and are calculated for all values :math:`0 \leq k_j < \text{resolution}`, thus a tensor of dimensions :math:`\text{resolution}^n` is returned. + This tensor can either be flattened into a vector or returned as an actual tensor with :math:`n` indices. + Note that is R = np.array([[1.0+0.0j]]) then :math:`H_k^{(R)}(y)` are precisely the well known **probabilists' Hermite polynomials** :math:`He_k(y)`:, + and if R = 2*np.array([[1.0+0.0j]]) then :math:`H_k^{(R)}(y)` are precisely the well known **physicists' Hermite polynomials** :math:`H_k(y)`:. + + Args: + R (array): Square matrix parametrizing the Hermite polynomial family + resolution (int): Maximum size of the subindices in the Hermite polynomial + y (array): Vector for the argument of the Hermite polynomial + renorm (bool): If True returns :math:`H_k^{(R)}(y)/\prod(\prod_i k_i!)` + make_tensor: If False returns a flattened one dimensional array instead of a tensor with the values of the polynomial + + Returns: + (array): The multidimensional Hermite polynomials. + """ + input_validation(R) + n, _ = R.shape + if y is None: + y = np.zeros([n], dtype=complex) + + m = y.shape[0] + if m != n: + raise ValueError("The matrix R and vector y have incompatible dimensions") + + values = np.array(hm(R, y, resolution, ren=renorm)) + + if make_tensor: + shape = resolution * np.ones([n], dtype=int) + values = np.reshape(values, shape) + + return values + + +def hafnian_batched(A, resolution, mu=None, tol=1e-12, renorm=False, make_tensor=True): + r"""Returns the haf(reduction(A, k)) where k is a vector of (non-negative) integers with the same dimensions as the square matrix A + :math:`k = (k_0,k_1,\ldots,k_{n-1})` and where :math:`0 \leq k_j < \text{resolution}`. + + If mu is not None the it instead calculates :math:`lhaf(fill_diagonal(reduction(A, k),reduction(mu, k)))`, this calculation can only be performed if + the matrix A has an inverse. + + Args: + R (array): Square matrix parametrizing + resolution (int): Maximum size of the subindices in the Hermite polynomial + y (array): Vector for the argument of the Hermite polynomial + renorm (bool): If True returns :math:`haf(reduction(A, k))/\prod(\prod_i k_i!)` or :math:`lhaf(fill_diagonal(reduction(A, k),reduction(mu, k)))` is mu is not None + make_tensor: If False returns a flattened one dimensional array instead of a tensor with the values of the polynomial + + Returns: + (array): The values of the hafnians. + """ + # pylint: disable=too-many-return-statements,too-many-branches,too-many-arguments + input_validation(A, tol=tol) + n, _ = A.shape + if not np.allclose(A, np.zeros([n, n])): + if mu is not None: + try: + yi = np.linalg.solve(A, mu) + except np.linalg.LinAlgError: + raise ValueError("The matrix does not have an inverse") + return hermite_multidimensional( + -A, resolution, y=-yi, renorm=renorm, make_tensor=make_tensor + ) + yi = np.zeros([n], dtype=complex) + return hermite_multidimensional( + -A, resolution, y=-yi, renorm=renorm, make_tensor=make_tensor + ) + # Note the minus signs in the arguments. Those are intentional + + if mu is None: + tensor = np.zeros([resolution ** n], dtype=complex) + tensor[0] = 1.0 + else: + index = resolution * np.ones([n], dtype=int) + tensor = np.empty(index, dtype=complex) + prim = np.array([expansion_coeff(alpha, resolution, renorm=renorm) for alpha in mu]) + for i in product(range(resolution), repeat=n): + tensor[i] = return_prod(prim, i) + if make_tensor: + return tensor + return tensor.flatten() diff --git a/hafnian/hafnian.cpp b/hafnian/hafnian.cpp index d0e36fd9d..1ffe1c696 100644 --- a/hafnian/hafnian.cpp +++ b/hafnian/hafnian.cpp @@ -1,13 +1,15 @@ -/* Generated by Cython 0.29.8 */ +/* Generated by Cython 0.29.4 */ /* BEGIN: Cython Metadata { "distutils": { "depends": [ + "src/batchhafnian.hpp", "src/eigenvalue_hafnian.hpp", "src/fsum.hpp", "src/hafnian.hpp", "src/hafnian_approx.hpp", + "src/permanent.hpp", "src/recursive_hafnian.hpp", "src/repeated_hafnian.hpp", "src/stdafx.h", @@ -19,33 +21,24 @@ "-Wall", "-fPIC", "-shared", - "-Xpreprocessor", "-fopenmp", - "-lomp", - "-mmacosx-version-min=10.9", - "-I/Users/bgupt/anaconda/lib/python3.6/site-packages/numpy/core/include" + "-I/home/nicolas/.local/lib/python3.6/site-packages/numpy/core/include" ], "extra_link_args": [ - "-Xpreprocessor -fopenmp -lomp" + "-fopenmp" ], "include_dirs": [ "hafnian", - "/opt/local/include/libomp/", - "/Users/bgupt/anaconda/lib/python3.6/site-packages/numpy/core/include", - "/Users/bgupt/eigen-eigen-323c052e1731", + "", + "/home/nicolas/.local/lib/python3.6/site-packages/numpy/core/include", "/usr/local/include/eigen3", "/usr/include/eigen3", - "src", - "/usr/local/opt/libomp/include" + "src" ], "language": "c++", - "libraries": [ - "omp" - ], "library_dirs": [ "/usr/lib", - "/usr/local/lib", - "/usr/local/opt/libomp/lib" + "/usr/local/lib" ], "name": "libhaf", "sources": [ @@ -63,8 +56,8 @@ END: Cython Metadata */ #elif PY_VERSION_HEX < 0x02060000 || (0x03000000 <= PY_VERSION_HEX && PY_VERSION_HEX < 0x03030000) #error Cython requires Python 2.6+ or Python 3.3+. #else -#define CYTHON_ABI "0_29_8" -#define CYTHON_HEX_VERSION 0x001D08F0 +#define CYTHON_ABI "0_29_4" +#define CYTHON_HEX_VERSION 0x001D04F0 #define CYTHON_FUTURE_DIVISION 0 #include #ifndef offsetof @@ -380,13 +373,8 @@ class __Pyx_FakeReference { #define __Pyx_DefaultClassType PyClass_Type #else #define __Pyx_BUILTIN_MODULE_NAME "builtins" -#if PY_VERSION_HEX < 0x030800A4 #define __Pyx_PyCode_New(a, k, l, s, f, code, c, n, v, fv, cell, fn, name, fline, lnos)\ PyCode_New(a, k, l, s, f, code, c, n, v, fv, cell, fn, name, fline, lnos) -#else - #define __Pyx_PyCode_New(a, k, l, s, f, code, c, n, v, fv, cell, fn, name, fline, lnos)\ - PyCode_New(a, 0, k, l, s, f, code, c, n, v, fv, cell, fn, name, fline, lnos) -#endif #define __Pyx_DefaultClassType PyType_Type #endif #ifndef Py_TPFLAGS_CHECKTYPES @@ -421,6 +409,26 @@ class __Pyx_FakeReference { #else #define __Pyx_PyFastCFunction_Check(func) 0 #endif +#if CYTHON_USE_DICT_VERSIONS +#define __PYX_GET_DICT_VERSION(dict) (((PyDictObject*)(dict))->ma_version_tag) +#define __PYX_UPDATE_DICT_CACHE(dict, value, cache_var, version_var)\ + (version_var) = __PYX_GET_DICT_VERSION(dict);\ + (cache_var) = (value); +#define __PYX_PY_DICT_LOOKUP_IF_MODIFIED(VAR, DICT, LOOKUP) {\ + static PY_UINT64_T __pyx_dict_version = 0;\ + static PyObject *__pyx_dict_cached_value = NULL;\ + if (likely(__PYX_GET_DICT_VERSION(DICT) == __pyx_dict_version)) {\ + (VAR) = __pyx_dict_cached_value;\ + } else {\ + (VAR) = __pyx_dict_cached_value = (LOOKUP);\ + __pyx_dict_version = __PYX_GET_DICT_VERSION(DICT);\ + }\ + } +#else +#define __PYX_GET_DICT_VERSION(dict) (0) +#define __PYX_UPDATE_DICT_CACHE(dict, value, cache_var, version_var) +#define __PYX_PY_DICT_LOOKUP_IF_MODIFIED(VAR, DICT, LOOKUP) (VAR) = (LOOKUP); +#endif #if CYTHON_COMPILING_IN_PYPY && !defined(PyObject_Malloc) #define PyObject_Malloc(s) PyMem_Malloc(s) #define PyObject_Free(p) PyMem_Free(p) @@ -686,8 +694,7 @@ typedef struct {PyObject **p; const char *s; const Py_ssize_t n; const char* enc const char is_unicode; const char is_str; const char intern; } __Pyx_StringTabEntry; #define __PYX_DEFAULT_STRING_ENCODING_IS_ASCII 0 -#define __PYX_DEFAULT_STRING_ENCODING_IS_UTF8 0 -#define __PYX_DEFAULT_STRING_ENCODING_IS_DEFAULT (PY_MAJOR_VERSION >= 3 && __PYX_DEFAULT_STRING_ENCODING_IS_UTF8) +#define __PYX_DEFAULT_STRING_ENCODING_IS_DEFAULT 0 #define __PYX_DEFAULT_STRING_ENCODING "" #define __Pyx_PyObject_FromString __Pyx_PyBytes_FromString #define __Pyx_PyObject_FromStringAndSize __Pyx_PyBytes_FromStringAndSize @@ -1274,6 +1281,23 @@ static CYTHON_INLINE int __pyx_sub_acquisition_count_locked( static CYTHON_INLINE void __Pyx_INC_MEMVIEW(__Pyx_memviewslice *, int, int); static CYTHON_INLINE void __Pyx_XDEC_MEMVIEW(__Pyx_memviewslice *, int, int); +/* ListCompAppend.proto */ +#if CYTHON_USE_PYLIST_INTERNALS && CYTHON_ASSUME_SAFE_MACROS +static CYTHON_INLINE int __Pyx_ListComp_Append(PyObject* list, PyObject* x) { + PyListObject* L = (PyListObject*) list; + Py_ssize_t len = Py_SIZE(list); + if (likely(L->allocated > len)) { + Py_INCREF(x); + PyList_SET_ITEM(list, len, x); + Py_SIZE(list) = len+1; + return 0; + } + return PyList_Append(list, x); +} +#else +#define __Pyx_ListComp_Append(L,x) PyList_Append(L,x) +#endif + /* ArgTypeTest.proto */ #define __Pyx_ArgTypeTest(obj, type, none_allowed, name, exact)\ ((likely((Py_TYPE(obj) == type) | (none_allowed && (obj == Py_None)))) ? 1 :\ @@ -1455,32 +1479,6 @@ static CYTHON_INLINE int __Pyx_PyErr_ExceptionMatchesInState(PyThreadState* tsta /* GetAttr3.proto */ static CYTHON_INLINE PyObject *__Pyx_GetAttr3(PyObject *, PyObject *, PyObject *); -/* PyDictVersioning.proto */ -#if CYTHON_USE_DICT_VERSIONS && CYTHON_USE_TYPE_SLOTS -#define __PYX_DICT_VERSION_INIT ((PY_UINT64_T) -1) -#define __PYX_GET_DICT_VERSION(dict) (((PyDictObject*)(dict))->ma_version_tag) -#define __PYX_UPDATE_DICT_CACHE(dict, value, cache_var, version_var)\ - (version_var) = __PYX_GET_DICT_VERSION(dict);\ - (cache_var) = (value); -#define __PYX_PY_DICT_LOOKUP_IF_MODIFIED(VAR, DICT, LOOKUP) {\ - static PY_UINT64_T __pyx_dict_version = 0;\ - static PyObject *__pyx_dict_cached_value = NULL;\ - if (likely(__PYX_GET_DICT_VERSION(DICT) == __pyx_dict_version)) {\ - (VAR) = __pyx_dict_cached_value;\ - } else {\ - (VAR) = __pyx_dict_cached_value = (LOOKUP);\ - __pyx_dict_version = __PYX_GET_DICT_VERSION(DICT);\ - }\ -} -static CYTHON_INLINE PY_UINT64_T __Pyx_get_tp_dict_version(PyObject *obj); -static CYTHON_INLINE PY_UINT64_T __Pyx_get_object_dict_version(PyObject *obj); -static CYTHON_INLINE int __Pyx_object_dict_version_matches(PyObject* obj, PY_UINT64_T tp_dict_version, PY_UINT64_T obj_dict_version); -#else -#define __PYX_GET_DICT_VERSION(dict) (0) -#define __PYX_UPDATE_DICT_CACHE(dict, value, cache_var, version_var) -#define __PYX_PY_DICT_LOOKUP_IF_MODIFIED(VAR, DICT, LOOKUP) (VAR) = (LOOKUP); -#endif - /* GetModuleGlobalName.proto */ #if CYTHON_USE_DICT_VERSIONS #define __Pyx_GetModuleGlobalName(var, name) {\ @@ -1563,23 +1561,6 @@ static CYTHON_INLINE int __Pyx_PyErr_GivenExceptionMatches2(PyObject *err, PyObj #define __Pyx_PyException_Check(obj) __Pyx_TypeCheck(obj, PyExc_Exception) static CYTHON_UNUSED int __pyx_memoryview_getbuffer(PyObject *__pyx_v_self, Py_buffer *__pyx_v_info, int __pyx_v_flags); /*proto*/ -/* ListCompAppend.proto */ -#if CYTHON_USE_PYLIST_INTERNALS && CYTHON_ASSUME_SAFE_MACROS -static CYTHON_INLINE int __Pyx_ListComp_Append(PyObject* list, PyObject* x) { - PyListObject* L = (PyListObject*) list; - Py_ssize_t len = Py_SIZE(list); - if (likely(L->allocated > len)) { - Py_INCREF(x); - PyList_SET_ITEM(list, len, x); - Py_SIZE(list) = len+1; - return 0; - } - return PyList_Append(list, x); -} -#else -#define __Pyx_ListComp_Append(L,x) PyList_Append(L,x) -#endif - /* PyIntBinop.proto */ #if !CYTHON_COMPILING_IN_PYPY static PyObject* __Pyx_PyInt_AddObjC(PyObject *op1, PyObject *op2, long intval, int inplace, int zerodivision_check); @@ -1876,6 +1857,9 @@ __pyx_memoryview_copy_new_contig(const __Pyx_memviewslice *from_mvs, /* CIntFromPy.proto */ static CYTHON_INLINE int __Pyx_PyInt_As_int(PyObject *); +/* CIntFromPy.proto */ +static CYTHON_INLINE size_t __Pyx_PyInt_As_size_t(PyObject *); + /* CIntFromPy.proto */ static CYTHON_INLINE long __Pyx_PyInt_As_long(PyObject *); @@ -1920,6 +1904,7 @@ static PyObject *contiguous = 0; static PyObject *indirect_contiguous = 0; static int __pyx_memoryview_thread_locks_used; static PyThread_type_lock __pyx_memoryview_thread_locks[8]; +static PyObject *__pyx_convert_vector_to_py___pyx_t_double_complex(const std::vector<__pyx_t_double_complex> &); /*proto*/ static struct __pyx_array_obj *__pyx_array_new(PyObject *, Py_ssize_t, char *, char *, char *); /*proto*/ static void *__pyx_align_pointer(void *, size_t); /*proto*/ static PyObject *__pyx_memoryview_new(PyObject *, int, int, __Pyx_TypeInfo *); /*proto*/ @@ -1984,6 +1969,7 @@ static const char __pyx_k_mat[] = "mat"; static const char __pyx_k_new[] = "__new__"; static const char __pyx_k_nud[] = "nud"; static const char __pyx_k_obj[] = "obj"; +static const char __pyx_k_ren[] = "ren"; static const char __pyx_k_rpt[] = "rpt"; static const char __pyx_k_base[] = "base"; static const char __pyx_k_dict[] = "__dict__"; @@ -2000,12 +1986,14 @@ static const char __pyx_k_step[] = "step"; static const char __pyx_k_stop[] = "stop"; static const char __pyx_k_test[] = "__test__"; static const char __pyx_k_ASCII[] = "ASCII"; +static const char __pyx_k_R_mat[] = "R_mat"; static const char __pyx_k_class[] = "__class__"; static const char __pyx_k_error[] = "error"; static const char __pyx_k_flags[] = "flags"; static const char __pyx_k_range[] = "range"; static const char __pyx_k_shape[] = "shape"; static const char __pyx_k_start[] = "start"; +static const char __pyx_k_y_mat[] = "y_mat"; static const char __pyx_k_approx[] = "approx"; static const char __pyx_k_encode[] = "encode"; static const char __pyx_k_format[] = "format"; @@ -2014,6 +2002,7 @@ static const char __pyx_k_libhaf[] = "libhaf"; static const char __pyx_k_name_2[] = "__name__"; static const char __pyx_k_pickle[] = "pickle"; static const char __pyx_k_reduce[] = "__reduce__"; +static const char __pyx_k_renorm[] = "renorm"; static const char __pyx_k_struct[] = "struct"; static const char __pyx_k_unpack[] = "unpack"; static const char __pyx_k_update[] = "update"; @@ -2037,6 +2026,7 @@ static const char __pyx_k_IndexError[] = "IndexError"; static const char __pyx_k_ValueError[] = "ValueError"; static const char __pyx_k_pyx_result[] = "__pyx_result"; static const char __pyx_k_pyx_vtable[] = "__pyx_vtable__"; +static const char __pyx_k_resolution[] = "resolution"; static const char __pyx_k_MemoryError[] = "MemoryError"; static const char __pyx_k_PickleError[] = "PickleError"; static const char __pyx_k_haf_complex[] = "haf_complex"; @@ -2064,6 +2054,7 @@ static const char __pyx_k_MemoryView_of_r_object[] = "" static const char __pyx_k_MemoryView_of_r_at_0x_x[] = ""; static const char __pyx_k_contiguous_and_indirect[] = ""; static const char __pyx_k_Cannot_index_with_type_s[] = "Cannot index with type '%s'"; +static const char __pyx_k_hermite_multidimensional[] = "hermite_multidimensional"; static const char __pyx_k_Invalid_shape_in_axis_d_d[] = "Invalid shape in axis %d: %d."; static const char __pyx_k_itemsize_0_for_cython_array[] = "itemsize <= 0 for cython.array"; static const char __pyx_k_unable_to_allocate_array_data[] = "unable to allocate array data."; @@ -2101,6 +2092,7 @@ static PyObject *__pyx_kp_s_MemoryView_of_r_object; static PyObject *__pyx_n_b_O; static PyObject *__pyx_kp_s_Out_of_bounds_on_buffer_access_a; static PyObject *__pyx_n_s_PickleError; +static PyObject *__pyx_n_s_R_mat; static PyObject *__pyx_n_s_TypeError; static PyObject *__pyx_kp_s_Unable_to_convert_item_to_object; static PyObject *__pyx_n_s_ValueError; @@ -2133,6 +2125,7 @@ static PyObject *__pyx_n_s_haf_real; static PyObject *__pyx_n_s_haf_rpt_complex; static PyObject *__pyx_n_s_haf_rpt_real; static PyObject *__pyx_kp_s_hafnian_hafnian_pyx; +static PyObject *__pyx_n_s_hermite_multidimensional; static PyObject *__pyx_n_s_i; static PyObject *__pyx_n_s_id; static PyObject *__pyx_n_s_import; @@ -2174,6 +2167,9 @@ static PyObject *__pyx_n_s_recursive; static PyObject *__pyx_n_s_reduce; static PyObject *__pyx_n_s_reduce_cython; static PyObject *__pyx_n_s_reduce_ex; +static PyObject *__pyx_n_s_ren; +static PyObject *__pyx_n_s_renorm; +static PyObject *__pyx_n_s_resolution; static PyObject *__pyx_n_s_rpt; static PyObject *__pyx_n_s_setstate; static PyObject *__pyx_n_s_setstate_cython; @@ -2194,6 +2190,7 @@ static PyObject *__pyx_kp_s_unable_to_allocate_array_data; static PyObject *__pyx_kp_s_unable_to_allocate_shape_and_str; static PyObject *__pyx_n_s_unpack; static PyObject *__pyx_n_s_update; +static PyObject *__pyx_n_s_y_mat; static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, PyObject *__pyx_v_fsum); /* proto */ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, PyObject *__pyx_v_fsum); /* proto */ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, __Pyx_memviewslice __pyx_v_rpt, __Pyx_memviewslice __pyx_v_mu, int __pyx_v_loop); /* proto */ @@ -2203,6 +2200,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, int __pyx_v_loop, int __pyx_v_recursive, PyObject *__pyx_v_quad, int __pyx_v_approx, PyObject *__pyx_v_nsamples); /* proto */ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, PyObject *__pyx_v_quad); /* proto */ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, PyObject *__pyx_v_quad, PyObject *__pyx_v_fsum); /* proto */ +static PyObject *__pyx_pf_6libhaf_18hermite_multidimensional(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, __Pyx_memviewslice __pyx_v_d, int __pyx_v_resolution, PyObject *__pyx_v_ren); /* proto */ static int __pyx_array___pyx_pf_15View_dot_MemoryView_5array___cinit__(struct __pyx_array_obj *__pyx_v_self, PyObject *__pyx_v_shape, Py_ssize_t __pyx_v_itemsize, PyObject *__pyx_v_format, PyObject *__pyx_v_mode, int __pyx_v_allocate_buffer); /* proto */ static int __pyx_array___pyx_pf_15View_dot_MemoryView_5array_2__getbuffer__(struct __pyx_array_obj *__pyx_v_self, Py_buffer *__pyx_v_info, int __pyx_v_flags); /* proto */ static void __pyx_array___pyx_pf_15View_dot_MemoryView_5array_4__dealloc__(struct __pyx_array_obj *__pyx_v_self); /* proto */ @@ -2284,11 +2282,12 @@ static PyObject *__pyx_tuple__33; static PyObject *__pyx_tuple__35; static PyObject *__pyx_tuple__37; static PyObject *__pyx_tuple__39; -static PyObject *__pyx_tuple__40; static PyObject *__pyx_tuple__41; static PyObject *__pyx_tuple__42; static PyObject *__pyx_tuple__43; static PyObject *__pyx_tuple__44; +static PyObject *__pyx_tuple__45; +static PyObject *__pyx_tuple__46; static PyObject *__pyx_codeobj__22; static PyObject *__pyx_codeobj__24; static PyObject *__pyx_codeobj__26; @@ -2298,10 +2297,11 @@ static PyObject *__pyx_codeobj__32; static PyObject *__pyx_codeobj__34; static PyObject *__pyx_codeobj__36; static PyObject *__pyx_codeobj__38; -static PyObject *__pyx_codeobj__45; +static PyObject *__pyx_codeobj__40; +static PyObject *__pyx_codeobj__47; /* Late includes */ -/* "hafnian/hafnian.pyx":54 +/* "hafnian/hafnian.pyx":56 * * * def torontonian_complex(double complex[:, :] A, fsum=False): # <<<<<<<<<<<<<< @@ -2347,7 +2347,7 @@ static PyObject *__pyx_pw_6libhaf_1torontonian_complex(PyObject *__pyx_self, PyO } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "torontonian_complex") < 0)) __PYX_ERR(0, 54, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "torontonian_complex") < 0)) __PYX_ERR(0, 56, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -2358,12 +2358,12 @@ static PyObject *__pyx_pw_6libhaf_1torontonian_complex(PyObject *__pyx_self, PyO default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 54, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 56, __pyx_L3_error) __pyx_v_fsum = values[1]; } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("torontonian_complex", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 54, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("torontonian_complex", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 56, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.torontonian_complex", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -2397,7 +2397,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ __pyx_t_double_complex __pyx_t_11; __Pyx_RefNannySetupContext("torontonian_complex", 0); - /* "hafnian/hafnian.pyx":72 + /* "hafnian/hafnian.pyx":74 * np.complex128: the torontonian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -2406,7 +2406,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":74 + /* "hafnian/hafnian.pyx":76 * cdef int i, j, n = A.shape[0] * cdef vector[double complex] mat * cdef int m = n/2 # <<<<<<<<<<<<<< @@ -2415,7 +2415,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ */ __pyx_v_m = __Pyx_div_long(__pyx_v_n, 2); - /* "hafnian/hafnian.pyx":76 + /* "hafnian/hafnian.pyx":78 * cdef int m = n/2 * * for i in range(n): # <<<<<<<<<<<<<< @@ -2427,7 +2427,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":77 + /* "hafnian/hafnian.pyx":79 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -2439,7 +2439,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":78 + /* "hafnian/hafnian.pyx":80 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -2452,22 +2452,22 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ __pyx_v_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 78, __pyx_L1_error) + __PYX_ERR(0, 80, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":80 + /* "hafnian/hafnian.pyx":82 * mat.push_back(A[i, j]) * * if fsum: # <<<<<<<<<<<<<< * return torontonian_fsum(mat) * */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 80, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 82, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":81 + /* "hafnian/hafnian.pyx":83 * * if fsum: * return torontonian_fsum(mat) # <<<<<<<<<<<<<< @@ -2475,13 +2475,13 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ * return torontonian_quad(mat) */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_fsum<__pyx_t_double_complex>(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 81, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_fsum<__pyx_t_double_complex>(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 83, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":80 + /* "hafnian/hafnian.pyx":82 * mat.push_back(A[i, j]) * * if fsum: # <<<<<<<<<<<<<< @@ -2490,7 +2490,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ */ } - /* "hafnian/hafnian.pyx":83 + /* "hafnian/hafnian.pyx":85 * return torontonian_fsum(mat) * * return torontonian_quad(mat) # <<<<<<<<<<<<<< @@ -2499,13 +2499,13 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ */ __Pyx_XDECREF(__pyx_r); __pyx_t_11 = hafnian::torontonian_quad(__pyx_v_mat); - __pyx_t_10 = __pyx_PyComplex_FromComplex(__pyx_t_11); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 83, __pyx_L1_error) + __pyx_t_10 = __pyx_PyComplex_FromComplex(__pyx_t_11); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 85, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":54 + /* "hafnian/hafnian.pyx":56 * * * def torontonian_complex(double complex[:, :] A, fsum=False): # <<<<<<<<<<<<<< @@ -2525,7 +2525,7 @@ static PyObject *__pyx_pf_6libhaf_torontonian_complex(CYTHON_UNUSED PyObject *__ return __pyx_r; } -/* "hafnian/hafnian.pyx":86 +/* "hafnian/hafnian.pyx":88 * * * def torontonian_real(double[:, :] A, fsum=False): # <<<<<<<<<<<<<< @@ -2571,7 +2571,7 @@ static PyObject *__pyx_pw_6libhaf_3torontonian_real(PyObject *__pyx_self, PyObje } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "torontonian_real") < 0)) __PYX_ERR(0, 86, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "torontonian_real") < 0)) __PYX_ERR(0, 88, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -2582,12 +2582,12 @@ static PyObject *__pyx_pw_6libhaf_3torontonian_real(PyObject *__pyx_self, PyObje default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 86, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 88, __pyx_L3_error) __pyx_v_fsum = values[1]; } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("torontonian_real", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 86, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("torontonian_real", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 88, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.torontonian_real", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -2620,7 +2620,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py PyObject *__pyx_t_10 = NULL; __Pyx_RefNannySetupContext("torontonian_real", 0); - /* "hafnian/hafnian.pyx":104 + /* "hafnian/hafnian.pyx":106 * np.float64: the torontonian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -2629,7 +2629,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":106 + /* "hafnian/hafnian.pyx":108 * cdef int i, j, n = A.shape[0] * cdef vector[double] mat * cdef int m = n/2 # <<<<<<<<<<<<<< @@ -2638,7 +2638,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py */ __pyx_v_m = __Pyx_div_long(__pyx_v_n, 2); - /* "hafnian/hafnian.pyx":108 + /* "hafnian/hafnian.pyx":110 * cdef int m = n/2 * * for i in range(n): # <<<<<<<<<<<<<< @@ -2650,7 +2650,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":109 + /* "hafnian/hafnian.pyx":111 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -2662,7 +2662,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":110 + /* "hafnian/hafnian.pyx":112 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -2675,22 +2675,22 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py __pyx_v_mat.push_back((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 110, __pyx_L1_error) + __PYX_ERR(0, 112, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":113 + /* "hafnian/hafnian.pyx":115 * * * if fsum: # <<<<<<<<<<<<<< * return torontonian_fsum(mat) * */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 113, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 115, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":114 + /* "hafnian/hafnian.pyx":116 * * if fsum: * return torontonian_fsum(mat) # <<<<<<<<<<<<<< @@ -2698,13 +2698,13 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py * return torontonian_quad(mat) */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_fsum(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 114, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_fsum(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 116, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":113 + /* "hafnian/hafnian.pyx":115 * * * if fsum: # <<<<<<<<<<<<<< @@ -2713,7 +2713,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py */ } - /* "hafnian/hafnian.pyx":116 + /* "hafnian/hafnian.pyx":118 * return torontonian_fsum(mat) * * return torontonian_quad(mat) # <<<<<<<<<<<<<< @@ -2721,13 +2721,13 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 116, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::torontonian_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 118, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":86 + /* "hafnian/hafnian.pyx":88 * * * def torontonian_real(double[:, :] A, fsum=False): # <<<<<<<<<<<<<< @@ -2747,7 +2747,7 @@ static PyObject *__pyx_pf_6libhaf_2torontonian_real(CYTHON_UNUSED PyObject *__py return __pyx_r; } -/* "hafnian/hafnian.pyx":123 +/* "hafnian/hafnian.pyx":125 * * * def haf_rpt_real(double[:, :] A, int[:] rpt, double[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< @@ -2794,7 +2794,7 @@ static PyObject *__pyx_pw_6libhaf_5haf_rpt_real(PyObject *__pyx_self, PyObject * case 1: if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_rpt)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("haf_rpt_real", 0, 2, 4, 1); __PYX_ERR(0, 123, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_rpt_real", 0, 2, 4, 1); __PYX_ERR(0, 125, __pyx_L3_error) } CYTHON_FALLTHROUGH; case 2: @@ -2810,7 +2810,7 @@ static PyObject *__pyx_pw_6libhaf_5haf_rpt_real(PyObject *__pyx_self, PyObject * } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_rpt_real") < 0)) __PYX_ERR(0, 123, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_rpt_real") < 0)) __PYX_ERR(0, 125, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -2824,23 +2824,23 @@ static PyObject *__pyx_pw_6libhaf_5haf_rpt_real(PyObject *__pyx_self, PyObject * default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 123, __pyx_L3_error) - __pyx_v_rpt = __Pyx_PyObject_to_MemoryviewSlice_ds_int(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_rpt.memview)) __PYX_ERR(0, 123, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 125, __pyx_L3_error) + __pyx_v_rpt = __Pyx_PyObject_to_MemoryviewSlice_ds_int(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_rpt.memview)) __PYX_ERR(0, 125, __pyx_L3_error) if (values[2]) { - __pyx_v_mu = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[2], PyBUF_WRITABLE); if (unlikely(!__pyx_v_mu.memview)) __PYX_ERR(0, 123, __pyx_L3_error) + __pyx_v_mu = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[2], PyBUF_WRITABLE); if (unlikely(!__pyx_v_mu.memview)) __PYX_ERR(0, 125, __pyx_L3_error) } else { __pyx_v_mu = __pyx_k_; __PYX_INC_MEMVIEW(&__pyx_v_mu, 1); } if (values[3]) { - __pyx_v_loop = __Pyx_PyObject_IsTrue(values[3]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 123, __pyx_L3_error) + __pyx_v_loop = __Pyx_PyObject_IsTrue(values[3]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 125, __pyx_L3_error) } else { __pyx_v_loop = ((int)0); } } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("haf_rpt_real", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 123, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_rpt_real", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 125, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.haf_rpt_real", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -2878,7 +2878,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se PyObject *__pyx_t_14 = NULL; __Pyx_RefNannySetupContext("haf_rpt_real", 0); - /* "hafnian/hafnian.pyx":140 + /* "hafnian/hafnian.pyx":142 * np.float64: the hafnian * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -2887,7 +2887,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":144 + /* "hafnian/hafnian.pyx":146 * cdef vector[double] mat, d * * for i in range(n): # <<<<<<<<<<<<<< @@ -2899,7 +2899,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":145 + /* "hafnian/hafnian.pyx":147 * * for i in range(n): * nud.push_back(rpt[i]) # <<<<<<<<<<<<<< @@ -2911,10 +2911,10 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_v_nud.push_back((*((int *) ( /* dim=0 */ (__pyx_v_rpt.data + __pyx_t_4 * __pyx_v_rpt.strides[0]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 145, __pyx_L1_error) + __PYX_ERR(0, 147, __pyx_L1_error) } - /* "hafnian/hafnian.pyx":147 + /* "hafnian/hafnian.pyx":149 * nud.push_back(rpt[i]) * * if mu is None: # <<<<<<<<<<<<<< @@ -2924,7 +2924,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_t_5 = ((((PyObject *) __pyx_v_mu.memview) == Py_None) != 0); if (__pyx_t_5) { - /* "hafnian/hafnian.pyx":148 + /* "hafnian/hafnian.pyx":150 * * if mu is None: * d.push_back(A[i, i]) # <<<<<<<<<<<<<< @@ -2937,10 +2937,10 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_v_d.push_back((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_6 * __pyx_v_A.strides[0]) ) + __pyx_t_7 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 148, __pyx_L1_error) + __PYX_ERR(0, 150, __pyx_L1_error) } - /* "hafnian/hafnian.pyx":147 + /* "hafnian/hafnian.pyx":149 * nud.push_back(rpt[i]) * * if mu is None: # <<<<<<<<<<<<<< @@ -2950,7 +2950,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se goto __pyx_L5; } - /* "hafnian/hafnian.pyx":150 + /* "hafnian/hafnian.pyx":152 * d.push_back(A[i, i]) * else: * d.push_back(mu[i]) # <<<<<<<<<<<<<< @@ -2963,12 +2963,12 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_v_d.push_back((*((double *) ( /* dim=0 */ (__pyx_v_mu.data + __pyx_t_8 * __pyx_v_mu.strides[0]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 150, __pyx_L1_error) + __PYX_ERR(0, 152, __pyx_L1_error) } } __pyx_L5:; - /* "hafnian/hafnian.pyx":152 + /* "hafnian/hafnian.pyx":154 * d.push_back(mu[i]) * * for j in range(n): # <<<<<<<<<<<<<< @@ -2980,7 +2980,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se for (__pyx_t_11 = 0; __pyx_t_11 < __pyx_t_10; __pyx_t_11+=1) { __pyx_v_j = __pyx_t_11; - /* "hafnian/hafnian.pyx":153 + /* "hafnian/hafnian.pyx":155 * * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -2993,12 +2993,12 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_v_mat.push_back((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_12 * __pyx_v_A.strides[0]) ) + __pyx_t_13 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 153, __pyx_L1_error) + __PYX_ERR(0, 155, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":156 + /* "hafnian/hafnian.pyx":158 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3008,7 +3008,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se __pyx_t_5 = (__pyx_v_loop != 0); if (__pyx_t_5) { - /* "hafnian/hafnian.pyx":157 + /* "hafnian/hafnian.pyx":159 * # Exposes a c function to python * if loop: * return loop_hafnian_rpt_quad(mat, d, nud) # <<<<<<<<<<<<<< @@ -3016,13 +3016,13 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se * return hafnian_rpt_quad(mat, nud) */ __Pyx_XDECREF(__pyx_r); - __pyx_t_14 = PyFloat_FromDouble(hafnian::loop_hafnian_rpt_quad(__pyx_v_mat, __pyx_v_d, __pyx_v_nud)); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 157, __pyx_L1_error) + __pyx_t_14 = PyFloat_FromDouble(hafnian::loop_hafnian_rpt_quad(__pyx_v_mat, __pyx_v_d, __pyx_v_nud)); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 159, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_14); __pyx_r = __pyx_t_14; __pyx_t_14 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":156 + /* "hafnian/hafnian.pyx":158 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3031,7 +3031,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se */ } - /* "hafnian/hafnian.pyx":159 + /* "hafnian/hafnian.pyx":161 * return loop_hafnian_rpt_quad(mat, d, nud) * * return hafnian_rpt_quad(mat, nud) # <<<<<<<<<<<<<< @@ -3039,13 +3039,13 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_14 = PyFloat_FromDouble(hafnian::hafnian_rpt_quad(__pyx_v_mat, __pyx_v_nud)); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 159, __pyx_L1_error) + __pyx_t_14 = PyFloat_FromDouble(hafnian::hafnian_rpt_quad(__pyx_v_mat, __pyx_v_nud)); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 161, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_14); __pyx_r = __pyx_t_14; __pyx_t_14 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":123 + /* "hafnian/hafnian.pyx":125 * * * def haf_rpt_real(double[:, :] A, int[:] rpt, double[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< @@ -3067,7 +3067,7 @@ static PyObject *__pyx_pf_6libhaf_4haf_rpt_real(CYTHON_UNUSED PyObject *__pyx_se return __pyx_r; } -/* "hafnian/hafnian.pyx":162 +/* "hafnian/hafnian.pyx":164 * * * def haf_rpt_complex(double complex[:, :] A, int[:] rpt, double complex[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< @@ -3114,7 +3114,7 @@ static PyObject *__pyx_pw_6libhaf_7haf_rpt_complex(PyObject *__pyx_self, PyObjec case 1: if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_rpt)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("haf_rpt_complex", 0, 2, 4, 1); __PYX_ERR(0, 162, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_rpt_complex", 0, 2, 4, 1); __PYX_ERR(0, 164, __pyx_L3_error) } CYTHON_FALLTHROUGH; case 2: @@ -3130,7 +3130,7 @@ static PyObject *__pyx_pw_6libhaf_7haf_rpt_complex(PyObject *__pyx_self, PyObjec } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_rpt_complex") < 0)) __PYX_ERR(0, 162, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_rpt_complex") < 0)) __PYX_ERR(0, 164, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -3144,23 +3144,23 @@ static PyObject *__pyx_pw_6libhaf_7haf_rpt_complex(PyObject *__pyx_self, PyObjec default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 162, __pyx_L3_error) - __pyx_v_rpt = __Pyx_PyObject_to_MemoryviewSlice_ds_int(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_rpt.memview)) __PYX_ERR(0, 162, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 164, __pyx_L3_error) + __pyx_v_rpt = __Pyx_PyObject_to_MemoryviewSlice_ds_int(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_rpt.memview)) __PYX_ERR(0, 164, __pyx_L3_error) if (values[2]) { - __pyx_v_mu = __Pyx_PyObject_to_MemoryviewSlice_ds___pyx_t_double_complex(values[2], PyBUF_WRITABLE); if (unlikely(!__pyx_v_mu.memview)) __PYX_ERR(0, 162, __pyx_L3_error) + __pyx_v_mu = __Pyx_PyObject_to_MemoryviewSlice_ds___pyx_t_double_complex(values[2], PyBUF_WRITABLE); if (unlikely(!__pyx_v_mu.memview)) __PYX_ERR(0, 164, __pyx_L3_error) } else { __pyx_v_mu = __pyx_k__2; __PYX_INC_MEMVIEW(&__pyx_v_mu, 1); } if (values[3]) { - __pyx_v_loop = __Pyx_PyObject_IsTrue(values[3]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 162, __pyx_L3_error) + __pyx_v_loop = __Pyx_PyObject_IsTrue(values[3]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 164, __pyx_L3_error) } else { __pyx_v_loop = ((int)0); } } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("haf_rpt_complex", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 162, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_rpt_complex", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 164, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.haf_rpt_complex", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -3199,7 +3199,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx PyObject *__pyx_t_15 = NULL; __Pyx_RefNannySetupContext("haf_rpt_complex", 0); - /* "hafnian/hafnian.pyx":179 + /* "hafnian/hafnian.pyx":181 * np.complex128: the hafnian * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -3208,7 +3208,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":183 + /* "hafnian/hafnian.pyx":185 * cdef vector[double complex] mat, d * * for i in range(n): # <<<<<<<<<<<<<< @@ -3220,7 +3220,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":184 + /* "hafnian/hafnian.pyx":186 * * for i in range(n): * nud.push_back(rpt[i]) # <<<<<<<<<<<<<< @@ -3232,10 +3232,10 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_v_nud.push_back((*((int *) ( /* dim=0 */ (__pyx_v_rpt.data + __pyx_t_4 * __pyx_v_rpt.strides[0]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 184, __pyx_L1_error) + __PYX_ERR(0, 186, __pyx_L1_error) } - /* "hafnian/hafnian.pyx":186 + /* "hafnian/hafnian.pyx":188 * nud.push_back(rpt[i]) * * if mu is None: # <<<<<<<<<<<<<< @@ -3245,7 +3245,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_t_5 = ((((PyObject *) __pyx_v_mu.memview) == Py_None) != 0); if (__pyx_t_5) { - /* "hafnian/hafnian.pyx":187 + /* "hafnian/hafnian.pyx":189 * * if mu is None: * d.push_back(A[i, i]) # <<<<<<<<<<<<<< @@ -3258,10 +3258,10 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_v_d.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_6 * __pyx_v_A.strides[0]) ) + __pyx_t_7 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 187, __pyx_L1_error) + __PYX_ERR(0, 189, __pyx_L1_error) } - /* "hafnian/hafnian.pyx":186 + /* "hafnian/hafnian.pyx":188 * nud.push_back(rpt[i]) * * if mu is None: # <<<<<<<<<<<<<< @@ -3271,7 +3271,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx goto __pyx_L5; } - /* "hafnian/hafnian.pyx":189 + /* "hafnian/hafnian.pyx":191 * d.push_back(A[i, i]) * else: * d.push_back(mu[i]) # <<<<<<<<<<<<<< @@ -3284,12 +3284,12 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_v_d.push_back((*((__pyx_t_double_complex *) ( /* dim=0 */ (__pyx_v_mu.data + __pyx_t_8 * __pyx_v_mu.strides[0]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 189, __pyx_L1_error) + __PYX_ERR(0, 191, __pyx_L1_error) } } __pyx_L5:; - /* "hafnian/hafnian.pyx":191 + /* "hafnian/hafnian.pyx":193 * d.push_back(mu[i]) * * for j in range(n): # <<<<<<<<<<<<<< @@ -3301,7 +3301,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx for (__pyx_t_11 = 0; __pyx_t_11 < __pyx_t_10; __pyx_t_11+=1) { __pyx_v_j = __pyx_t_11; - /* "hafnian/hafnian.pyx":192 + /* "hafnian/hafnian.pyx":194 * * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -3314,12 +3314,12 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_v_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_12 * __pyx_v_A.strides[0]) ) + __pyx_t_13 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 192, __pyx_L1_error) + __PYX_ERR(0, 194, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":195 + /* "hafnian/hafnian.pyx":197 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3329,7 +3329,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx __pyx_t_5 = (__pyx_v_loop != 0); if (__pyx_t_5) { - /* "hafnian/hafnian.pyx":196 + /* "hafnian/hafnian.pyx":198 * # Exposes a c function to python * if loop: * return loop_hafnian_rpt_quad(mat, d, nud) # <<<<<<<<<<<<<< @@ -3338,13 +3338,13 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx */ __Pyx_XDECREF(__pyx_r); __pyx_t_14 = hafnian::loop_hafnian_rpt_quad(__pyx_v_mat, __pyx_v_d, __pyx_v_nud); - __pyx_t_15 = __pyx_PyComplex_FromComplex(__pyx_t_14); if (unlikely(!__pyx_t_15)) __PYX_ERR(0, 196, __pyx_L1_error) + __pyx_t_15 = __pyx_PyComplex_FromComplex(__pyx_t_14); if (unlikely(!__pyx_t_15)) __PYX_ERR(0, 198, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_15); __pyx_r = __pyx_t_15; __pyx_t_15 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":195 + /* "hafnian/hafnian.pyx":197 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3353,7 +3353,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx */ } - /* "hafnian/hafnian.pyx":198 + /* "hafnian/hafnian.pyx":200 * return loop_hafnian_rpt_quad(mat, d, nud) * * return hafnian_rpt_quad(mat, nud) # <<<<<<<<<<<<<< @@ -3362,13 +3362,13 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx */ __Pyx_XDECREF(__pyx_r); __pyx_t_14 = hafnian::hafnian_rpt_quad(__pyx_v_mat, __pyx_v_nud); - __pyx_t_15 = __pyx_PyComplex_FromComplex(__pyx_t_14); if (unlikely(!__pyx_t_15)) __PYX_ERR(0, 198, __pyx_L1_error) + __pyx_t_15 = __pyx_PyComplex_FromComplex(__pyx_t_14); if (unlikely(!__pyx_t_15)) __PYX_ERR(0, 200, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_15); __pyx_r = __pyx_t_15; __pyx_t_15 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":162 + /* "hafnian/hafnian.pyx":164 * * * def haf_rpt_complex(double complex[:, :] A, int[:] rpt, double complex[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< @@ -3390,7 +3390,7 @@ static PyObject *__pyx_pf_6libhaf_6haf_rpt_complex(CYTHON_UNUSED PyObject *__pyx return __pyx_r; } -/* "hafnian/hafnian.pyx":205 +/* "hafnian/hafnian.pyx":207 * * * def haf_int(long long[:, :] A): # <<<<<<<<<<<<<< @@ -3408,7 +3408,7 @@ static PyObject *__pyx_pw_6libhaf_9haf_int(PyObject *__pyx_self, PyObject *__pyx __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("haf_int (wrapper)", 0); assert(__pyx_arg_A); { - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_PY_LONG_LONG(__pyx_arg_A, PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 205, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_PY_LONG_LONG(__pyx_arg_A, PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 207, __pyx_L3_error) } goto __pyx_L4_argument_unpacking_done; __pyx_L3_error:; @@ -3441,7 +3441,7 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ PyObject *__pyx_t_9 = NULL; __Pyx_RefNannySetupContext("haf_int", 0); - /* "hafnian/hafnian.pyx":217 + /* "hafnian/hafnian.pyx":219 * np.int64: the hafnian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -3450,7 +3450,7 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":220 + /* "hafnian/hafnian.pyx":222 * cdef vector[long long] mat * * for i in range(n): # <<<<<<<<<<<<<< @@ -3462,7 +3462,7 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":221 + /* "hafnian/hafnian.pyx":223 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -3474,7 +3474,7 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":222 + /* "hafnian/hafnian.pyx":224 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -3487,12 +3487,12 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ __pyx_v_mat.push_back((*((PY_LONG_LONG *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 222, __pyx_L1_error) + __PYX_ERR(0, 224, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":225 + /* "hafnian/hafnian.pyx":227 * * # Exposes a c function to python * return hafnian_recursive(mat) # <<<<<<<<<<<<<< @@ -3500,13 +3500,13 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_9 = __Pyx_PyInt_From_PY_LONG_LONG(hafnian::hafnian_recursive(__pyx_v_mat)); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 225, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyInt_From_PY_LONG_LONG(hafnian::hafnian_recursive(__pyx_v_mat)); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 227, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __pyx_r = __pyx_t_9; __pyx_t_9 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":205 + /* "hafnian/hafnian.pyx":207 * * * def haf_int(long long[:, :] A): # <<<<<<<<<<<<<< @@ -3526,7 +3526,7 @@ static PyObject *__pyx_pf_6libhaf_8haf_int(CYTHON_UNUSED PyObject *__pyx_self, _ return __pyx_r; } -/* "hafnian/hafnian.pyx":232 +/* "hafnian/hafnian.pyx":234 * * * def haf_complex(double complex[:, :] A, bint loop=False, bint recursive=True, quad=True): # <<<<<<<<<<<<<< @@ -3590,7 +3590,7 @@ static PyObject *__pyx_pw_6libhaf_11haf_complex(PyObject *__pyx_self, PyObject * } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_complex") < 0)) __PYX_ERR(0, 232, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_complex") < 0)) __PYX_ERR(0, 234, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -3605,14 +3605,14 @@ static PyObject *__pyx_pw_6libhaf_11haf_complex(PyObject *__pyx_self, PyObject * default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 232, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 234, __pyx_L3_error) if (values[1]) { - __pyx_v_loop = __Pyx_PyObject_IsTrue(values[1]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 232, __pyx_L3_error) + __pyx_v_loop = __Pyx_PyObject_IsTrue(values[1]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 234, __pyx_L3_error) } else { __pyx_v_loop = ((int)0); } if (values[2]) { - __pyx_v_recursive = __Pyx_PyObject_IsTrue(values[2]); if (unlikely((__pyx_v_recursive == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 232, __pyx_L3_error) + __pyx_v_recursive = __Pyx_PyObject_IsTrue(values[2]); if (unlikely((__pyx_v_recursive == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 234, __pyx_L3_error) } else { __pyx_v_recursive = ((int)1); } @@ -3620,7 +3620,7 @@ static PyObject *__pyx_pw_6libhaf_11haf_complex(PyObject *__pyx_self, PyObject * } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("haf_complex", 0, 1, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 232, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_complex", 0, 1, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 234, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.haf_complex", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -3653,7 +3653,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se PyObject *__pyx_t_11 = NULL; __Pyx_RefNannySetupContext("haf_complex", 0); - /* "hafnian/hafnian.pyx":246 + /* "hafnian/hafnian.pyx":248 * np.complex128: the hafnian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -3662,7 +3662,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":249 + /* "hafnian/hafnian.pyx":251 * cdef vector[double complex] mat * * for i in range(n): # <<<<<<<<<<<<<< @@ -3674,7 +3674,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":250 + /* "hafnian/hafnian.pyx":252 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -3686,7 +3686,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":251 + /* "hafnian/hafnian.pyx":253 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -3699,12 +3699,12 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se __pyx_v_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 251, __pyx_L1_error) + __PYX_ERR(0, 253, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":254 + /* "hafnian/hafnian.pyx":256 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3714,7 +3714,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se __pyx_t_9 = (__pyx_v_loop != 0); if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":255 + /* "hafnian/hafnian.pyx":257 * # Exposes a c function to python * if loop: * return loop_hafnian(mat) # <<<<<<<<<<<<<< @@ -3723,13 +3723,13 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::loop_hafnian<__pyx_t_double_complex>(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 255, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 257, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":254 + /* "hafnian/hafnian.pyx":256 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -3738,7 +3738,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ } - /* "hafnian/hafnian.pyx":257 + /* "hafnian/hafnian.pyx":259 * return loop_hafnian(mat) * * if recursive: # <<<<<<<<<<<<<< @@ -3748,17 +3748,17 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se __pyx_t_9 = (__pyx_v_recursive != 0); if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":258 + /* "hafnian/hafnian.pyx":260 * * if recursive: * if quad: # <<<<<<<<<<<<<< * return hafnian_recursive_quad(mat) * return hafnian_recursive(mat) */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 258, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 260, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":259 + /* "hafnian/hafnian.pyx":261 * if recursive: * if quad: * return hafnian_recursive_quad(mat) # <<<<<<<<<<<<<< @@ -3767,13 +3767,13 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::hafnian_recursive_quad(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 259, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 261, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":258 + /* "hafnian/hafnian.pyx":260 * * if recursive: * if quad: # <<<<<<<<<<<<<< @@ -3782,7 +3782,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ } - /* "hafnian/hafnian.pyx":260 + /* "hafnian/hafnian.pyx":262 * if quad: * return hafnian_recursive_quad(mat) * return hafnian_recursive(mat) # <<<<<<<<<<<<<< @@ -3791,13 +3791,13 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::hafnian_recursive<__pyx_t_double_complex>(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 260, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 262, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":257 + /* "hafnian/hafnian.pyx":259 * return loop_hafnian(mat) * * if recursive: # <<<<<<<<<<<<<< @@ -3806,7 +3806,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ } - /* "hafnian/hafnian.pyx":262 + /* "hafnian/hafnian.pyx":264 * return hafnian_recursive(mat) * * return hafnian(mat) # <<<<<<<<<<<<<< @@ -3815,13 +3815,13 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::hafnian<__pyx_t_double_complex>(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 262, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 264, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":232 + /* "hafnian/hafnian.pyx":234 * * * def haf_complex(double complex[:, :] A, bint loop=False, bint recursive=True, quad=True): # <<<<<<<<<<<<<< @@ -3841,7 +3841,7 @@ static PyObject *__pyx_pf_6libhaf_10haf_complex(CYTHON_UNUSED PyObject *__pyx_se return __pyx_r; } -/* "hafnian/hafnian.pyx":265 +/* "hafnian/hafnian.pyx":267 * * * def haf_real(double[:, :] A, bint loop=False, bint recursive=True, quad=True, bint approx=False, nsamples=1000): # <<<<<<<<<<<<<< @@ -3924,7 +3924,7 @@ static PyObject *__pyx_pw_6libhaf_13haf_real(PyObject *__pyx_self, PyObject *__p } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_real") < 0)) __PYX_ERR(0, 265, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "haf_real") < 0)) __PYX_ERR(0, 267, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -3943,20 +3943,20 @@ static PyObject *__pyx_pw_6libhaf_13haf_real(PyObject *__pyx_self, PyObject *__p default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 265, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 267, __pyx_L3_error) if (values[1]) { - __pyx_v_loop = __Pyx_PyObject_IsTrue(values[1]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 265, __pyx_L3_error) + __pyx_v_loop = __Pyx_PyObject_IsTrue(values[1]); if (unlikely((__pyx_v_loop == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 267, __pyx_L3_error) } else { __pyx_v_loop = ((int)0); } if (values[2]) { - __pyx_v_recursive = __Pyx_PyObject_IsTrue(values[2]); if (unlikely((__pyx_v_recursive == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 265, __pyx_L3_error) + __pyx_v_recursive = __Pyx_PyObject_IsTrue(values[2]); if (unlikely((__pyx_v_recursive == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 267, __pyx_L3_error) } else { __pyx_v_recursive = ((int)1); } __pyx_v_quad = values[3]; if (values[4]) { - __pyx_v_approx = __Pyx_PyObject_IsTrue(values[4]); if (unlikely((__pyx_v_approx == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 265, __pyx_L3_error) + __pyx_v_approx = __Pyx_PyObject_IsTrue(values[4]); if (unlikely((__pyx_v_approx == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 267, __pyx_L3_error) } else { __pyx_v_approx = ((int)0); } @@ -3964,7 +3964,7 @@ static PyObject *__pyx_pw_6libhaf_13haf_real(PyObject *__pyx_self, PyObject *__p } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("haf_real", 0, 1, 6, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 265, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("haf_real", 0, 1, 6, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 267, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.haf_real", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -3996,7 +3996,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_t_10 = NULL; __Pyx_RefNannySetupContext("haf_real", 0); - /* "hafnian/hafnian.pyx":283 + /* "hafnian/hafnian.pyx":285 * np.float64: the hafnian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -4005,7 +4005,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":286 + /* "hafnian/hafnian.pyx":288 * cdef vector[double] mat * * for i in range(n): # <<<<<<<<<<<<<< @@ -4017,7 +4017,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":287 + /* "hafnian/hafnian.pyx":289 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -4029,7 +4029,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":288 + /* "hafnian/hafnian.pyx":290 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -4042,12 +4042,12 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, __pyx_v_mat.push_back((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 288, __pyx_L1_error) + __PYX_ERR(0, 290, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":291 + /* "hafnian/hafnian.pyx":293 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -4057,7 +4057,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, __pyx_t_9 = (__pyx_v_loop != 0); if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":292 + /* "hafnian/hafnian.pyx":294 * # Exposes a c function to python * if loop: * return loop_hafnian(mat) # <<<<<<<<<<<<<< @@ -4065,13 +4065,13 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, * if approx: */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::loop_hafnian(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 292, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::loop_hafnian(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":291 + /* "hafnian/hafnian.pyx":293 * * # Exposes a c function to python * if loop: # <<<<<<<<<<<<<< @@ -4080,7 +4080,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, */ } - /* "hafnian/hafnian.pyx":294 + /* "hafnian/hafnian.pyx":296 * return loop_hafnian(mat) * * if approx: # <<<<<<<<<<<<<< @@ -4090,7 +4090,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, __pyx_t_9 = (__pyx_v_approx != 0); if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":295 + /* "hafnian/hafnian.pyx":297 * * if approx: * return hafnian_approx(mat, nsamples) # <<<<<<<<<<<<<< @@ -4098,14 +4098,14 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, * if recursive: */ __Pyx_XDECREF(__pyx_r); - __pyx_t_1 = __Pyx_PyInt_As_int(__pyx_v_nsamples); if (unlikely((__pyx_t_1 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 295, __pyx_L1_error) - __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_approx(__pyx_v_mat, __pyx_t_1)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 295, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyInt_As_int(__pyx_v_nsamples); if (unlikely((__pyx_t_1 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 297, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_approx(__pyx_v_mat, __pyx_t_1)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 297, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":294 + /* "hafnian/hafnian.pyx":296 * return loop_hafnian(mat) * * if approx: # <<<<<<<<<<<<<< @@ -4114,7 +4114,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, */ } - /* "hafnian/hafnian.pyx":297 + /* "hafnian/hafnian.pyx":299 * return hafnian_approx(mat, nsamples) * * if recursive: # <<<<<<<<<<<<<< @@ -4124,17 +4124,17 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, __pyx_t_9 = (__pyx_v_recursive != 0); if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":298 + /* "hafnian/hafnian.pyx":300 * * if recursive: * if quad: # <<<<<<<<<<<<<< * return hafnian_recursive_quad(mat) * return hafnian_recursive(mat) */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 298, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 300, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":299 + /* "hafnian/hafnian.pyx":301 * if recursive: * if quad: * return hafnian_recursive_quad(mat) # <<<<<<<<<<<<<< @@ -4142,13 +4142,13 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_recursive_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 299, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_recursive_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 301, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":298 + /* "hafnian/hafnian.pyx":300 * * if recursive: * if quad: # <<<<<<<<<<<<<< @@ -4157,7 +4157,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, */ } - /* "hafnian/hafnian.pyx":300 + /* "hafnian/hafnian.pyx":302 * if quad: * return hafnian_recursive_quad(mat) * return hafnian_recursive(mat) # <<<<<<<<<<<<<< @@ -4165,13 +4165,13 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, * return hafnian(mat) */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_recursive(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 300, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian_recursive(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":297 + /* "hafnian/hafnian.pyx":299 * return hafnian_approx(mat, nsamples) * * if recursive: # <<<<<<<<<<<<<< @@ -4180,7 +4180,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, */ } - /* "hafnian/hafnian.pyx":302 + /* "hafnian/hafnian.pyx":304 * return hafnian_recursive(mat) * * return hafnian(mat) # <<<<<<<<<<<<<< @@ -4188,13 +4188,13 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 302, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::hafnian(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 304, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":265 + /* "hafnian/hafnian.pyx":267 * * * def haf_real(double[:, :] A, bint loop=False, bint recursive=True, quad=True, bint approx=False, nsamples=1000): # <<<<<<<<<<<<<< @@ -4214,7 +4214,7 @@ static PyObject *__pyx_pf_6libhaf_12haf_real(CYTHON_UNUSED PyObject *__pyx_self, return __pyx_r; } -/* "hafnian/hafnian.pyx":310 +/* "hafnian/hafnian.pyx":312 * * * def perm_complex(double complex[:, :] A, quad=True): # <<<<<<<<<<<<<< @@ -4260,7 +4260,7 @@ static PyObject *__pyx_pw_6libhaf_15perm_complex(PyObject *__pyx_self, PyObject } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "perm_complex") < 0)) __PYX_ERR(0, 310, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "perm_complex") < 0)) __PYX_ERR(0, 312, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -4271,12 +4271,12 @@ static PyObject *__pyx_pw_6libhaf_15perm_complex(PyObject *__pyx_self, PyObject default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 310, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 312, __pyx_L3_error) __pyx_v_quad = values[1]; } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("perm_complex", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 310, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("perm_complex", 0, 1, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 312, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.perm_complex", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -4309,7 +4309,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s PyObject *__pyx_t_11 = NULL; __Pyx_RefNannySetupContext("perm_complex", 0); - /* "hafnian/hafnian.pyx":321 + /* "hafnian/hafnian.pyx":323 * np.complex128: the hafnian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -4318,7 +4318,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":324 + /* "hafnian/hafnian.pyx":326 * cdef vector[double complex] mat * * for i in range(n): # <<<<<<<<<<<<<< @@ -4330,7 +4330,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":325 + /* "hafnian/hafnian.pyx":327 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -4342,7 +4342,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":326 + /* "hafnian/hafnian.pyx":328 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -4355,22 +4355,22 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s __pyx_v_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 326, __pyx_L1_error) + __PYX_ERR(0, 328, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":329 + /* "hafnian/hafnian.pyx":331 * * # Exposes a c function to python * if quad: # <<<<<<<<<<<<<< * return permanent_quad(mat) * */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 329, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 331, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":330 + /* "hafnian/hafnian.pyx":332 * # Exposes a c function to python * if quad: * return permanent_quad(mat) # <<<<<<<<<<<<<< @@ -4379,13 +4379,13 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::permanent_quad(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 330, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 332, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":329 + /* "hafnian/hafnian.pyx":331 * * # Exposes a c function to python * if quad: # <<<<<<<<<<<<<< @@ -4394,7 +4394,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s */ } - /* "hafnian/hafnian.pyx":332 + /* "hafnian/hafnian.pyx":334 * return permanent_quad(mat) * * return permanent(mat) # <<<<<<<<<<<<<< @@ -4403,13 +4403,13 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s */ __Pyx_XDECREF(__pyx_r); __pyx_t_10 = hafnian::permanent<__pyx_t_double_complex>(__pyx_v_mat); - __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 332, __pyx_L1_error) + __pyx_t_11 = __pyx_PyComplex_FromComplex(__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 334, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_r = __pyx_t_11; __pyx_t_11 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":310 + /* "hafnian/hafnian.pyx":312 * * * def perm_complex(double complex[:, :] A, quad=True): # <<<<<<<<<<<<<< @@ -4429,7 +4429,7 @@ static PyObject *__pyx_pf_6libhaf_14perm_complex(CYTHON_UNUSED PyObject *__pyx_s return __pyx_r; } -/* "hafnian/hafnian.pyx":335 +/* "hafnian/hafnian.pyx":337 * * * def perm_real(double [:, :] A, quad=True, fsum=False): # <<<<<<<<<<<<<< @@ -4485,7 +4485,7 @@ static PyObject *__pyx_pw_6libhaf_17perm_real(PyObject *__pyx_self, PyObject *__ } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "perm_real") < 0)) __PYX_ERR(0, 335, __pyx_L3_error) + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "perm_real") < 0)) __PYX_ERR(0, 337, __pyx_L3_error) } } else { switch (PyTuple_GET_SIZE(__pyx_args)) { @@ -4498,13 +4498,13 @@ static PyObject *__pyx_pw_6libhaf_17perm_real(PyObject *__pyx_self, PyObject *__ default: goto __pyx_L5_argtuple_error; } } - __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 335, __pyx_L3_error) + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 337, __pyx_L3_error) __pyx_v_quad = values[1]; __pyx_v_fsum = values[2]; } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("perm_real", 0, 1, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 335, __pyx_L3_error) + __Pyx_RaiseArgtupleInvalid("perm_real", 0, 1, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 337, __pyx_L3_error) __pyx_L3_error:; __Pyx_AddTraceback("libhaf.perm_real", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); @@ -4536,7 +4536,7 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self PyObject *__pyx_t_10 = NULL; __Pyx_RefNannySetupContext("perm_real", 0); - /* "hafnian/hafnian.pyx":348 + /* "hafnian/hafnian.pyx":350 * np.float64: the hafnian of matrix A * """ * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< @@ -4545,7 +4545,7 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self */ __pyx_v_n = (__pyx_v_A.shape[0]); - /* "hafnian/hafnian.pyx":351 + /* "hafnian/hafnian.pyx":353 * cdef vector[double] mat * * for i in range(n): # <<<<<<<<<<<<<< @@ -4557,7 +4557,7 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) { __pyx_v_i = __pyx_t_3; - /* "hafnian/hafnian.pyx":352 + /* "hafnian/hafnian.pyx":354 * * for i in range(n): * for j in range(n): # <<<<<<<<<<<<<< @@ -4569,7 +4569,7 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) { __pyx_v_j = __pyx_t_6; - /* "hafnian/hafnian.pyx":353 + /* "hafnian/hafnian.pyx":355 * for i in range(n): * for j in range(n): * mat.push_back(A[i, j]) # <<<<<<<<<<<<<< @@ -4582,22 +4582,22 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self __pyx_v_mat.push_back((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_7 * __pyx_v_A.strides[0]) ) + __pyx_t_8 * __pyx_v_A.strides[1]) )))); } catch(...) { __Pyx_CppExn2PyErr(); - __PYX_ERR(0, 353, __pyx_L1_error) + __PYX_ERR(0, 355, __pyx_L1_error) } } } - /* "hafnian/hafnian.pyx":355 + /* "hafnian/hafnian.pyx":357 * mat.push_back(A[i, j]) * * if fsum: # <<<<<<<<<<<<<< * return permanent_fsum(mat) * */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 355, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_fsum); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 357, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":356 + /* "hafnian/hafnian.pyx":358 * * if fsum: * return permanent_fsum(mat) # <<<<<<<<<<<<<< @@ -4605,13 +4605,13 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self * # Exposes a c function to python */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent_fsum(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 356, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent_fsum(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 358, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":355 + /* "hafnian/hafnian.pyx":357 * mat.push_back(A[i, j]) * * if fsum: # <<<<<<<<<<<<<< @@ -4620,17 +4620,17 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self */ } - /* "hafnian/hafnian.pyx":359 + /* "hafnian/hafnian.pyx":361 * * # Exposes a c function to python * if quad: # <<<<<<<<<<<<<< * return permanent_quad(mat) * */ - __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 359, __pyx_L1_error) + __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_v_quad); if (unlikely(__pyx_t_9 < 0)) __PYX_ERR(0, 361, __pyx_L1_error) if (__pyx_t_9) { - /* "hafnian/hafnian.pyx":360 + /* "hafnian/hafnian.pyx":362 * # Exposes a c function to python * if quad: * return permanent_quad(mat) # <<<<<<<<<<<<<< @@ -4638,49 +4638,382 @@ static PyObject *__pyx_pf_6libhaf_16perm_real(CYTHON_UNUSED PyObject *__pyx_self * return permanent(mat) */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 360, __pyx_L1_error) + __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent_quad(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 362, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_r = __pyx_t_10; __pyx_t_10 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":359 + /* "hafnian/hafnian.pyx":361 + * + * # Exposes a c function to python + * if quad: # <<<<<<<<<<<<<< + * return permanent_quad(mat) + * + */ + } + + /* "hafnian/hafnian.pyx":364 + * return permanent_quad(mat) + * + * return permanent(mat) # <<<<<<<<<<<<<< + * + * + */ + __Pyx_XDECREF(__pyx_r); + __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 364, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_10); + __pyx_r = __pyx_t_10; + __pyx_t_10 = 0; + goto __pyx_L0; + + /* "hafnian/hafnian.pyx":337 + * + * + * def perm_real(double [:, :] A, quad=True, fsum=False): # <<<<<<<<<<<<<< + * """Returns the hafnian of a real matrix A via the C++ hafnian library. + * + */ + + /* function exit code */ + __pyx_L1_error:; + __Pyx_XDECREF(__pyx_t_10); + __Pyx_AddTraceback("libhaf.perm_real", __pyx_clineno, __pyx_lineno, __pyx_filename); + __pyx_r = NULL; + __pyx_L0:; + __PYX_XDEC_MEMVIEW(&__pyx_v_A, 1); + __Pyx_XGIVEREF(__pyx_r); + __Pyx_RefNannyFinishContext(); + return __pyx_r; +} + +/* "hafnian/hafnian.pyx":370 + * # Batch hafnian + * + * def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): # <<<<<<<<<<<<<< + * cdef int i, j, n = A.shape[0] + * cdef vector[double complex] R_mat, y_mat + */ + +/* Python wrapper */ +static PyObject *__pyx_pw_6libhaf_19hermite_multidimensional(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/ +static char __pyx_doc_6libhaf_18hermite_multidimensional[] = "hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False)"; +static PyMethodDef __pyx_mdef_6libhaf_19hermite_multidimensional = {"hermite_multidimensional", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6libhaf_19hermite_multidimensional, METH_VARARGS|METH_KEYWORDS, __pyx_doc_6libhaf_18hermite_multidimensional}; +static PyObject *__pyx_pw_6libhaf_19hermite_multidimensional(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) { + __Pyx_memviewslice __pyx_v_A = { 0, 0, { 0 }, { 0 }, { 0 } }; + __Pyx_memviewslice __pyx_v_d = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_v_resolution; + PyObject *__pyx_v_ren = 0; + PyObject *__pyx_r = 0; + __Pyx_RefNannyDeclarations + __Pyx_RefNannySetupContext("hermite_multidimensional (wrapper)", 0); + { + static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_A,&__pyx_n_s_d,&__pyx_n_s_resolution,&__pyx_n_s_ren,0}; + PyObject* values[4] = {0,0,0,0}; + values[3] = ((PyObject *)Py_False); + if (unlikely(__pyx_kwds)) { + Py_ssize_t kw_args; + const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args); + switch (pos_args) { + case 4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3); + CYTHON_FALLTHROUGH; + case 3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2); + CYTHON_FALLTHROUGH; + case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1); + CYTHON_FALLTHROUGH; + case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0); + CYTHON_FALLTHROUGH; + case 0: break; + default: goto __pyx_L5_argtuple_error; + } + kw_args = PyDict_Size(__pyx_kwds); + switch (pos_args) { + case 0: + if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_A)) != 0)) kw_args--; + else goto __pyx_L5_argtuple_error; + CYTHON_FALLTHROUGH; + case 1: + if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_d)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("hermite_multidimensional", 0, 3, 4, 1); __PYX_ERR(0, 370, __pyx_L3_error) + } + CYTHON_FALLTHROUGH; + case 2: + if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_resolution)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("hermite_multidimensional", 0, 3, 4, 2); __PYX_ERR(0, 370, __pyx_L3_error) + } + CYTHON_FALLTHROUGH; + case 3: + if (kw_args > 0) { + PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_ren); + if (value) { values[3] = value; kw_args--; } + } + } + if (unlikely(kw_args > 0)) { + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "hermite_multidimensional") < 0)) __PYX_ERR(0, 370, __pyx_L3_error) + } + } else { + switch (PyTuple_GET_SIZE(__pyx_args)) { + case 4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3); + CYTHON_FALLTHROUGH; + case 3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2); + values[1] = PyTuple_GET_ITEM(__pyx_args, 1); + values[0] = PyTuple_GET_ITEM(__pyx_args, 0); + break; + default: goto __pyx_L5_argtuple_error; + } + } + __pyx_v_A = __Pyx_PyObject_to_MemoryviewSlice_dsds___pyx_t_double_complex(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 370, __pyx_L3_error) + __pyx_v_d = __Pyx_PyObject_to_MemoryviewSlice_ds___pyx_t_double_complex(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_d.memview)) __PYX_ERR(0, 370, __pyx_L3_error) + __pyx_v_resolution = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_resolution == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 370, __pyx_L3_error) + __pyx_v_ren = values[3]; + } + goto __pyx_L4_argument_unpacking_done; + __pyx_L5_argtuple_error:; + __Pyx_RaiseArgtupleInvalid("hermite_multidimensional", 0, 3, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 370, __pyx_L3_error) + __pyx_L3_error:; + __Pyx_AddTraceback("libhaf.hermite_multidimensional", __pyx_clineno, __pyx_lineno, __pyx_filename); + __Pyx_RefNannyFinishContext(); + return NULL; + __pyx_L4_argument_unpacking_done:; + __pyx_r = __pyx_pf_6libhaf_18hermite_multidimensional(__pyx_self, __pyx_v_A, __pyx_v_d, __pyx_v_resolution, __pyx_v_ren); + + /* function exit code */ + __Pyx_RefNannyFinishContext(); + return __pyx_r; +} + +static PyObject *__pyx_pf_6libhaf_18hermite_multidimensional(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, __Pyx_memviewslice __pyx_v_d, int __pyx_v_resolution, PyObject *__pyx_v_ren) { + int __pyx_v_i; + int __pyx_v_j; + int __pyx_v_n; + std::vector<__pyx_t_double_complex> __pyx_v_R_mat; + std::vector<__pyx_t_double_complex> __pyx_v_y_mat; + int __pyx_v_renorm; + PyObject *__pyx_r = NULL; + __Pyx_RefNannyDeclarations + int __pyx_t_1; + int __pyx_t_2; + int __pyx_t_3; + int __pyx_t_4; + int __pyx_t_5; + int __pyx_t_6; + int __pyx_t_7; + Py_ssize_t __pyx_t_8; + Py_ssize_t __pyx_t_9; + Py_ssize_t __pyx_t_10; + PyObject *__pyx_t_11 = NULL; + __Pyx_RefNannySetupContext("hermite_multidimensional", 0); + + /* "hafnian/hafnian.pyx":371 + * + * def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): + * cdef int i, j, n = A.shape[0] # <<<<<<<<<<<<<< + * cdef vector[double complex] R_mat, y_mat + * + */ + __pyx_v_n = (__pyx_v_A.shape[0]); + + /* "hafnian/hafnian.pyx":374 + * cdef vector[double complex] R_mat, y_mat + * + * cdef int renorm = 0 # <<<<<<<<<<<<<< + * + * if ren: + */ + __pyx_v_renorm = 0; + + /* "hafnian/hafnian.pyx":376 + * cdef int renorm = 0 + * + * if ren: # <<<<<<<<<<<<<< + * renorm = 1 + * + */ + __pyx_t_1 = __Pyx_PyObject_IsTrue(__pyx_v_ren); if (unlikely(__pyx_t_1 < 0)) __PYX_ERR(0, 376, __pyx_L1_error) + if (__pyx_t_1) { + + /* "hafnian/hafnian.pyx":377 + * + * if ren: + * renorm = 1 # <<<<<<<<<<<<<< + * + * for i in range(n): + */ + __pyx_v_renorm = 1; + + /* "hafnian/hafnian.pyx":376 + * cdef int renorm = 0 + * + * if ren: # <<<<<<<<<<<<<< + * renorm = 1 + * + */ + } + + /* "hafnian/hafnian.pyx":379 + * renorm = 1 + * + * for i in range(n): # <<<<<<<<<<<<<< + * for j in range(n): + * R_mat.push_back(A[i, j]) + */ + __pyx_t_2 = __pyx_v_n; + __pyx_t_3 = __pyx_t_2; + for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) { + __pyx_v_i = __pyx_t_4; + + /* "hafnian/hafnian.pyx":380 + * + * for i in range(n): + * for j in range(n): # <<<<<<<<<<<<<< + * R_mat.push_back(A[i, j]) + * + */ + __pyx_t_5 = __pyx_v_n; + __pyx_t_6 = __pyx_t_5; + for (__pyx_t_7 = 0; __pyx_t_7 < __pyx_t_6; __pyx_t_7+=1) { + __pyx_v_j = __pyx_t_7; + + /* "hafnian/hafnian.pyx":381 + * for i in range(n): + * for j in range(n): + * R_mat.push_back(A[i, j]) # <<<<<<<<<<<<<< + * + * + */ + __pyx_t_8 = __pyx_v_i; + __pyx_t_9 = __pyx_v_j; + try { + __pyx_v_R_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_8 * __pyx_v_A.strides[0]) ) + __pyx_t_9 * __pyx_v_A.strides[1]) )))); + } catch(...) { + __Pyx_CppExn2PyErr(); + __PYX_ERR(0, 381, __pyx_L1_error) + } + } + } + + /* "hafnian/hafnian.pyx":384 + * + * + * for i in range(n): # <<<<<<<<<<<<<< + * y_mat.push_back(d[i]) + * + */ + __pyx_t_2 = __pyx_v_n; + __pyx_t_3 = __pyx_t_2; + for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) { + __pyx_v_i = __pyx_t_4; + + /* "hafnian/hafnian.pyx":385 + * + * for i in range(n): + * y_mat.push_back(d[i]) # <<<<<<<<<<<<<< + * + * return hermite_multidimensional_cpp(R_mat, y_mat, resolution, renorm) + */ + __pyx_t_10 = __pyx_v_i; + try { + __pyx_v_y_mat.push_back((*((__pyx_t_double_complex *) ( /* dim=0 */ (__pyx_v_d.data + __pyx_t_10 * __pyx_v_d.strides[0]) )))); + } catch(...) { + __Pyx_CppExn2PyErr(); + __PYX_ERR(0, 385, __pyx_L1_error) + } + } + + /* "hafnian/hafnian.pyx":387 + * y_mat.push_back(d[i]) + * + * return hermite_multidimensional_cpp(R_mat, y_mat, resolution, renorm) # <<<<<<<<<<<<<< + */ + __Pyx_XDECREF(__pyx_r); + __pyx_t_11 = __pyx_convert_vector_to_py___pyx_t_double_complex(hafnian::hermite_multidimensional_cpp(__pyx_v_R_mat, __pyx_v_y_mat, __pyx_v_resolution, __pyx_v_renorm)); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 387, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_11); + __pyx_r = __pyx_t_11; + __pyx_t_11 = 0; + goto __pyx_L0; + + /* "hafnian/hafnian.pyx":370 + * # Batch hafnian + * + * def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): # <<<<<<<<<<<<<< + * cdef int i, j, n = A.shape[0] + * cdef vector[double complex] R_mat, y_mat + */ + + /* function exit code */ + __pyx_L1_error:; + __Pyx_XDECREF(__pyx_t_11); + __Pyx_AddTraceback("libhaf.hermite_multidimensional", __pyx_clineno, __pyx_lineno, __pyx_filename); + __pyx_r = NULL; + __pyx_L0:; + __PYX_XDEC_MEMVIEW(&__pyx_v_A, 1); + __PYX_XDEC_MEMVIEW(&__pyx_v_d, 1); + __Pyx_XGIVEREF(__pyx_r); + __Pyx_RefNannyFinishContext(); + return __pyx_r; +} + +/* "vector.to_py":60 * - * # Exposes a c function to python - * if quad: # <<<<<<<<<<<<<< - * return permanent_quad(mat) + * @cname("__pyx_convert_vector_to_py___pyx_t_double_complex") + * cdef object __pyx_convert_vector_to_py___pyx_t_double_complex(vector[X]& v): # <<<<<<<<<<<<<< + * return [v[i] for i in range(v.size())] * */ - } - /* "hafnian/hafnian.pyx":362 - * return permanent_quad(mat) +static PyObject *__pyx_convert_vector_to_py___pyx_t_double_complex(const std::vector<__pyx_t_double_complex> &__pyx_v_v) { + size_t __pyx_v_i; + PyObject *__pyx_r = NULL; + __Pyx_RefNannyDeclarations + PyObject *__pyx_t_1 = NULL; + size_t __pyx_t_2; + size_t __pyx_t_3; + size_t __pyx_t_4; + __pyx_t_double_complex __pyx_t_5; + PyObject *__pyx_t_6 = NULL; + __Pyx_RefNannySetupContext("__pyx_convert_vector_to_py___pyx_t_double_complex", 0); + + /* "vector.to_py":61 + * @cname("__pyx_convert_vector_to_py___pyx_t_double_complex") + * cdef object __pyx_convert_vector_to_py___pyx_t_double_complex(vector[X]& v): + * return [v[i] for i in range(v.size())] # <<<<<<<<<<<<<< * - * return permanent(mat) # <<<<<<<<<<<<<< * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_10 = PyFloat_FromDouble(hafnian::permanent(__pyx_v_mat)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 362, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_10); - __pyx_r = __pyx_t_10; - __pyx_t_10 = 0; + __pyx_t_1 = PyList_New(0); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 61, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_1); + __pyx_t_2 = __pyx_v_v.size(); + __pyx_t_3 = __pyx_t_2; + for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) { + __pyx_v_i = __pyx_t_4; + __pyx_t_5 = (__pyx_v_v[__pyx_v_i]); + __pyx_t_6 = __pyx_PyComplex_FromComplex(__pyx_t_5); if (unlikely(!__pyx_t_6)) __PYX_ERR(1, 61, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_6); + if (unlikely(__Pyx_ListComp_Append(__pyx_t_1, (PyObject*)__pyx_t_6))) __PYX_ERR(1, 61, __pyx_L1_error) + __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; + } + __pyx_r = __pyx_t_1; + __pyx_t_1 = 0; goto __pyx_L0; - /* "hafnian/hafnian.pyx":335 - * + /* "vector.to_py":60 * - * def perm_real(double [:, :] A, quad=True, fsum=False): # <<<<<<<<<<<<<< - * """Returns the hafnian of a real matrix A via the C++ hafnian library. + * @cname("__pyx_convert_vector_to_py___pyx_t_double_complex") + * cdef object __pyx_convert_vector_to_py___pyx_t_double_complex(vector[X]& v): # <<<<<<<<<<<<<< + * return [v[i] for i in range(v.size())] * */ /* function exit code */ __pyx_L1_error:; - __Pyx_XDECREF(__pyx_t_10); - __Pyx_AddTraceback("libhaf.perm_real", __pyx_clineno, __pyx_lineno, __pyx_filename); - __pyx_r = NULL; + __Pyx_XDECREF(__pyx_t_1); + __Pyx_XDECREF(__pyx_t_6); + __Pyx_AddTraceback("vector.to_py.__pyx_convert_vector_to_py___pyx_t_double_complex", __pyx_clineno, __pyx_lineno, __pyx_filename); + __pyx_r = 0; __pyx_L0:; - __PYX_XDEC_MEMVIEW(&__pyx_v_A, 1); __Pyx_XGIVEREF(__pyx_r); __Pyx_RefNannyFinishContext(); return __pyx_r; @@ -18197,6 +18530,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_b_O, __pyx_k_O, sizeof(__pyx_k_O), 0, 0, 0, 1}, {&__pyx_kp_s_Out_of_bounds_on_buffer_access_a, __pyx_k_Out_of_bounds_on_buffer_access_a, sizeof(__pyx_k_Out_of_bounds_on_buffer_access_a), 0, 0, 1, 0}, {&__pyx_n_s_PickleError, __pyx_k_PickleError, sizeof(__pyx_k_PickleError), 0, 0, 1, 1}, + {&__pyx_n_s_R_mat, __pyx_k_R_mat, sizeof(__pyx_k_R_mat), 0, 0, 1, 1}, {&__pyx_n_s_TypeError, __pyx_k_TypeError, sizeof(__pyx_k_TypeError), 0, 0, 1, 1}, {&__pyx_kp_s_Unable_to_convert_item_to_object, __pyx_k_Unable_to_convert_item_to_object, sizeof(__pyx_k_Unable_to_convert_item_to_object), 0, 0, 1, 0}, {&__pyx_n_s_ValueError, __pyx_k_ValueError, sizeof(__pyx_k_ValueError), 0, 0, 1, 1}, @@ -18229,6 +18563,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_haf_rpt_complex, __pyx_k_haf_rpt_complex, sizeof(__pyx_k_haf_rpt_complex), 0, 0, 1, 1}, {&__pyx_n_s_haf_rpt_real, __pyx_k_haf_rpt_real, sizeof(__pyx_k_haf_rpt_real), 0, 0, 1, 1}, {&__pyx_kp_s_hafnian_hafnian_pyx, __pyx_k_hafnian_hafnian_pyx, sizeof(__pyx_k_hafnian_hafnian_pyx), 0, 0, 1, 0}, + {&__pyx_n_s_hermite_multidimensional, __pyx_k_hermite_multidimensional, sizeof(__pyx_k_hermite_multidimensional), 0, 0, 1, 1}, {&__pyx_n_s_i, __pyx_k_i, sizeof(__pyx_k_i), 0, 0, 1, 1}, {&__pyx_n_s_id, __pyx_k_id, sizeof(__pyx_k_id), 0, 0, 1, 1}, {&__pyx_n_s_import, __pyx_k_import, sizeof(__pyx_k_import), 0, 0, 1, 1}, @@ -18270,6 +18605,9 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_reduce, __pyx_k_reduce, sizeof(__pyx_k_reduce), 0, 0, 1, 1}, {&__pyx_n_s_reduce_cython, __pyx_k_reduce_cython, sizeof(__pyx_k_reduce_cython), 0, 0, 1, 1}, {&__pyx_n_s_reduce_ex, __pyx_k_reduce_ex, sizeof(__pyx_k_reduce_ex), 0, 0, 1, 1}, + {&__pyx_n_s_ren, __pyx_k_ren, sizeof(__pyx_k_ren), 0, 0, 1, 1}, + {&__pyx_n_s_renorm, __pyx_k_renorm, sizeof(__pyx_k_renorm), 0, 0, 1, 1}, + {&__pyx_n_s_resolution, __pyx_k_resolution, sizeof(__pyx_k_resolution), 0, 0, 1, 1}, {&__pyx_n_s_rpt, __pyx_k_rpt, sizeof(__pyx_k_rpt), 0, 0, 1, 1}, {&__pyx_n_s_setstate, __pyx_k_setstate, sizeof(__pyx_k_setstate), 0, 0, 1, 1}, {&__pyx_n_s_setstate_cython, __pyx_k_setstate_cython, sizeof(__pyx_k_setstate_cython), 0, 0, 1, 1}, @@ -18290,10 +18628,11 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_kp_s_unable_to_allocate_shape_and_str, __pyx_k_unable_to_allocate_shape_and_str, sizeof(__pyx_k_unable_to_allocate_shape_and_str), 0, 0, 1, 0}, {&__pyx_n_s_unpack, __pyx_k_unpack, sizeof(__pyx_k_unpack), 0, 0, 1, 1}, {&__pyx_n_s_update, __pyx_k_update, sizeof(__pyx_k_update), 0, 0, 1, 1}, + {&__pyx_n_s_y_mat, __pyx_k_y_mat, sizeof(__pyx_k_y_mat), 0, 0, 1, 1}, {0, 0, 0, 0, 0, 0, 0} }; static CYTHON_SMALL_CODE int __Pyx_InitCachedBuiltins(void) { - __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 76, __pyx_L1_error) + __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 78, __pyx_L1_error) __pyx_builtin_ValueError = __Pyx_GetBuiltinName(__pyx_n_s_ValueError); if (!__pyx_builtin_ValueError) __PYX_ERR(1, 133, __pyx_L1_error) __pyx_builtin_MemoryError = __Pyx_GetBuiltinName(__pyx_n_s_MemoryError); if (!__pyx_builtin_MemoryError) __PYX_ERR(1, 148, __pyx_L1_error) __pyx_builtin_enumerate = __Pyx_GetBuiltinName(__pyx_n_s_enumerate); if (!__pyx_builtin_enumerate) __PYX_ERR(1, 151, __pyx_L1_error) @@ -18502,113 +18841,125 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { __Pyx_GOTREF(__pyx_tuple__20); __Pyx_GIVEREF(__pyx_tuple__20); - /* "hafnian/hafnian.pyx":54 + /* "hafnian/hafnian.pyx":56 * * * def torontonian_complex(double complex[:, :] A, fsum=False): # <<<<<<<<<<<<<< * """Returns the Torontonian of a complex matrix A via the C++ hafnian library. * */ - __pyx_tuple__21 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat, __pyx_n_s_m); if (unlikely(!__pyx_tuple__21)) __PYX_ERR(0, 54, __pyx_L1_error) + __pyx_tuple__21 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat, __pyx_n_s_m); if (unlikely(!__pyx_tuple__21)) __PYX_ERR(0, 56, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__21); __Pyx_GIVEREF(__pyx_tuple__21); - __pyx_codeobj__22 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__21, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_torontonian_complex, 54, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__22)) __PYX_ERR(0, 54, __pyx_L1_error) + __pyx_codeobj__22 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__21, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_torontonian_complex, 56, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__22)) __PYX_ERR(0, 56, __pyx_L1_error) - /* "hafnian/hafnian.pyx":86 + /* "hafnian/hafnian.pyx":88 * * * def torontonian_real(double[:, :] A, fsum=False): # <<<<<<<<<<<<<< * """Returns the Torontonian of a real matrix A via the C++ hafnian library. * */ - __pyx_tuple__23 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat, __pyx_n_s_m); if (unlikely(!__pyx_tuple__23)) __PYX_ERR(0, 86, __pyx_L1_error) + __pyx_tuple__23 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat, __pyx_n_s_m); if (unlikely(!__pyx_tuple__23)) __PYX_ERR(0, 88, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__23); __Pyx_GIVEREF(__pyx_tuple__23); - __pyx_codeobj__24 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__23, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_torontonian_real, 86, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__24)) __PYX_ERR(0, 86, __pyx_L1_error) + __pyx_codeobj__24 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__23, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_torontonian_real, 88, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__24)) __PYX_ERR(0, 88, __pyx_L1_error) - /* "hafnian/hafnian.pyx":123 + /* "hafnian/hafnian.pyx":125 * * * def haf_rpt_real(double[:, :] A, int[:] rpt, double[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< * r"""Returns the hafnian of a real matrix A via the C++ hafnian library * using the rpt method. This method is more efficient for matrices with */ - __pyx_tuple__25 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_rpt, __pyx_n_s_mu, __pyx_n_s_loop, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_nud, __pyx_n_s_mat, __pyx_n_s_d); if (unlikely(!__pyx_tuple__25)) __PYX_ERR(0, 123, __pyx_L1_error) + __pyx_tuple__25 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_rpt, __pyx_n_s_mu, __pyx_n_s_loop, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_nud, __pyx_n_s_mat, __pyx_n_s_d); if (unlikely(!__pyx_tuple__25)) __PYX_ERR(0, 125, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__25); __Pyx_GIVEREF(__pyx_tuple__25); - __pyx_codeobj__26 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__25, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_rpt_real, 123, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__26)) __PYX_ERR(0, 123, __pyx_L1_error) + __pyx_codeobj__26 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__25, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_rpt_real, 125, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__26)) __PYX_ERR(0, 125, __pyx_L1_error) - /* "hafnian/hafnian.pyx":162 + /* "hafnian/hafnian.pyx":164 * * * def haf_rpt_complex(double complex[:, :] A, int[:] rpt, double complex[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< * r"""Returns the hafnian of a complex matrix A via the C++ hafnian library * using the rpt method. This method is more efficient for matrices with */ - __pyx_tuple__27 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_rpt, __pyx_n_s_mu, __pyx_n_s_loop, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_nud, __pyx_n_s_mat, __pyx_n_s_d); if (unlikely(!__pyx_tuple__27)) __PYX_ERR(0, 162, __pyx_L1_error) + __pyx_tuple__27 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_rpt, __pyx_n_s_mu, __pyx_n_s_loop, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_nud, __pyx_n_s_mat, __pyx_n_s_d); if (unlikely(!__pyx_tuple__27)) __PYX_ERR(0, 164, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__27); __Pyx_GIVEREF(__pyx_tuple__27); - __pyx_codeobj__28 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__27, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_rpt_complex, 162, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__28)) __PYX_ERR(0, 162, __pyx_L1_error) + __pyx_codeobj__28 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__27, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_rpt_complex, 164, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__28)) __PYX_ERR(0, 164, __pyx_L1_error) - /* "hafnian/hafnian.pyx":205 + /* "hafnian/hafnian.pyx":207 * * * def haf_int(long long[:, :] A): # <<<<<<<<<<<<<< * """Returns the hafnian of an integer matrix A via the C++ hafnian library. * Modified with permission from https://github.com/eklotek/Hafnian. */ - __pyx_tuple__29 = PyTuple_Pack(6, __pyx_n_s_A, __pyx_n_s_A, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__29)) __PYX_ERR(0, 205, __pyx_L1_error) + __pyx_tuple__29 = PyTuple_Pack(6, __pyx_n_s_A, __pyx_n_s_A, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__29)) __PYX_ERR(0, 207, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__29); __Pyx_GIVEREF(__pyx_tuple__29); - __pyx_codeobj__30 = (PyObject*)__Pyx_PyCode_New(1, 0, 6, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__29, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_int, 205, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__30)) __PYX_ERR(0, 205, __pyx_L1_error) + __pyx_codeobj__30 = (PyObject*)__Pyx_PyCode_New(1, 0, 6, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__29, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_int, 207, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__30)) __PYX_ERR(0, 207, __pyx_L1_error) - /* "hafnian/hafnian.pyx":232 + /* "hafnian/hafnian.pyx":234 * * * def haf_complex(double complex[:, :] A, bint loop=False, bint recursive=True, quad=True): # <<<<<<<<<<<<<< * """Returns the hafnian of a complex matrix A via the C++ hafnian library. * */ - __pyx_tuple__31 = PyTuple_Pack(8, __pyx_n_s_A, __pyx_n_s_loop, __pyx_n_s_recursive, __pyx_n_s_quad, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__31)) __PYX_ERR(0, 232, __pyx_L1_error) + __pyx_tuple__31 = PyTuple_Pack(8, __pyx_n_s_A, __pyx_n_s_loop, __pyx_n_s_recursive, __pyx_n_s_quad, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__31)) __PYX_ERR(0, 234, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__31); __Pyx_GIVEREF(__pyx_tuple__31); - __pyx_codeobj__32 = (PyObject*)__Pyx_PyCode_New(4, 0, 8, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__31, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_complex, 232, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__32)) __PYX_ERR(0, 232, __pyx_L1_error) + __pyx_codeobj__32 = (PyObject*)__Pyx_PyCode_New(4, 0, 8, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__31, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_complex, 234, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__32)) __PYX_ERR(0, 234, __pyx_L1_error) - /* "hafnian/hafnian.pyx":265 + /* "hafnian/hafnian.pyx":267 * * * def haf_real(double[:, :] A, bint loop=False, bint recursive=True, quad=True, bint approx=False, nsamples=1000): # <<<<<<<<<<<<<< * """Returns the hafnian of a real matrix A via the C++ hafnian library. * */ - __pyx_tuple__33 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_loop, __pyx_n_s_recursive, __pyx_n_s_quad, __pyx_n_s_approx, __pyx_n_s_nsamples, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__33)) __PYX_ERR(0, 265, __pyx_L1_error) + __pyx_tuple__33 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_loop, __pyx_n_s_recursive, __pyx_n_s_quad, __pyx_n_s_approx, __pyx_n_s_nsamples, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__33)) __PYX_ERR(0, 267, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__33); __Pyx_GIVEREF(__pyx_tuple__33); - __pyx_codeobj__34 = (PyObject*)__Pyx_PyCode_New(6, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__33, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_real, 265, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__34)) __PYX_ERR(0, 265, __pyx_L1_error) + __pyx_codeobj__34 = (PyObject*)__Pyx_PyCode_New(6, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__33, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_haf_real, 267, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__34)) __PYX_ERR(0, 267, __pyx_L1_error) - /* "hafnian/hafnian.pyx":310 + /* "hafnian/hafnian.pyx":312 * * * def perm_complex(double complex[:, :] A, quad=True): # <<<<<<<<<<<<<< * """Returns the hafnian of a complex matrix A via the C++ hafnian library. * */ - __pyx_tuple__35 = PyTuple_Pack(6, __pyx_n_s_A, __pyx_n_s_quad, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__35)) __PYX_ERR(0, 310, __pyx_L1_error) + __pyx_tuple__35 = PyTuple_Pack(6, __pyx_n_s_A, __pyx_n_s_quad, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__35)) __PYX_ERR(0, 312, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__35); __Pyx_GIVEREF(__pyx_tuple__35); - __pyx_codeobj__36 = (PyObject*)__Pyx_PyCode_New(2, 0, 6, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__35, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_perm_complex, 310, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__36)) __PYX_ERR(0, 310, __pyx_L1_error) + __pyx_codeobj__36 = (PyObject*)__Pyx_PyCode_New(2, 0, 6, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__35, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_perm_complex, 312, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__36)) __PYX_ERR(0, 312, __pyx_L1_error) - /* "hafnian/hafnian.pyx":335 + /* "hafnian/hafnian.pyx":337 * * * def perm_real(double [:, :] A, quad=True, fsum=False): # <<<<<<<<<<<<<< * """Returns the hafnian of a real matrix A via the C++ hafnian library. * */ - __pyx_tuple__37 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_quad, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__37)) __PYX_ERR(0, 335, __pyx_L1_error) + __pyx_tuple__37 = PyTuple_Pack(7, __pyx_n_s_A, __pyx_n_s_quad, __pyx_n_s_fsum, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_mat); if (unlikely(!__pyx_tuple__37)) __PYX_ERR(0, 337, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__37); __Pyx_GIVEREF(__pyx_tuple__37); - __pyx_codeobj__38 = (PyObject*)__Pyx_PyCode_New(3, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__37, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_perm_real, 335, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__38)) __PYX_ERR(0, 335, __pyx_L1_error) + __pyx_codeobj__38 = (PyObject*)__Pyx_PyCode_New(3, 0, 7, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__37, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_perm_real, 337, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__38)) __PYX_ERR(0, 337, __pyx_L1_error) + + /* "hafnian/hafnian.pyx":370 + * # Batch hafnian + * + * def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): # <<<<<<<<<<<<<< + * cdef int i, j, n = A.shape[0] + * cdef vector[double complex] R_mat, y_mat + */ + __pyx_tuple__39 = PyTuple_Pack(10, __pyx_n_s_A, __pyx_n_s_d, __pyx_n_s_resolution, __pyx_n_s_ren, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_R_mat, __pyx_n_s_y_mat, __pyx_n_s_renorm); if (unlikely(!__pyx_tuple__39)) __PYX_ERR(0, 370, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__39); + __Pyx_GIVEREF(__pyx_tuple__39); + __pyx_codeobj__40 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__39, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_hafnian_hafnian_pyx, __pyx_n_s_hermite_multidimensional, 370, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__40)) __PYX_ERR(0, 370, __pyx_L1_error) /* "View.MemoryView":286 * return self.name @@ -18617,9 +18968,9 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { * cdef strided = Enum("") # default * cdef indirect = Enum("") */ - __pyx_tuple__39 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct_or_indirect); if (unlikely(!__pyx_tuple__39)) __PYX_ERR(1, 286, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__39); - __Pyx_GIVEREF(__pyx_tuple__39); + __pyx_tuple__41 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct_or_indirect); if (unlikely(!__pyx_tuple__41)) __PYX_ERR(1, 286, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__41); + __Pyx_GIVEREF(__pyx_tuple__41); /* "View.MemoryView":287 * @@ -18628,9 +18979,9 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { * cdef indirect = Enum("") * */ - __pyx_tuple__40 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct); if (unlikely(!__pyx_tuple__40)) __PYX_ERR(1, 287, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__40); - __Pyx_GIVEREF(__pyx_tuple__40); + __pyx_tuple__42 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct); if (unlikely(!__pyx_tuple__42)) __PYX_ERR(1, 287, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__42); + __Pyx_GIVEREF(__pyx_tuple__42); /* "View.MemoryView":288 * cdef generic = Enum("") @@ -18639,9 +18990,9 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { * * */ - __pyx_tuple__41 = PyTuple_Pack(1, __pyx_kp_s_strided_and_indirect); if (unlikely(!__pyx_tuple__41)) __PYX_ERR(1, 288, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__41); - __Pyx_GIVEREF(__pyx_tuple__41); + __pyx_tuple__43 = PyTuple_Pack(1, __pyx_kp_s_strided_and_indirect); if (unlikely(!__pyx_tuple__43)) __PYX_ERR(1, 288, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__43); + __Pyx_GIVEREF(__pyx_tuple__43); /* "View.MemoryView":291 * @@ -18650,9 +19001,9 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { * cdef indirect_contiguous = Enum("") * */ - __pyx_tuple__42 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_direct); if (unlikely(!__pyx_tuple__42)) __PYX_ERR(1, 291, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__42); - __Pyx_GIVEREF(__pyx_tuple__42); + __pyx_tuple__44 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_direct); if (unlikely(!__pyx_tuple__44)) __PYX_ERR(1, 291, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__44); + __Pyx_GIVEREF(__pyx_tuple__44); /* "View.MemoryView":292 * @@ -18661,19 +19012,19 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) { * * */ - __pyx_tuple__43 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_indirect); if (unlikely(!__pyx_tuple__43)) __PYX_ERR(1, 292, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__43); - __Pyx_GIVEREF(__pyx_tuple__43); + __pyx_tuple__45 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_indirect); if (unlikely(!__pyx_tuple__45)) __PYX_ERR(1, 292, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__45); + __Pyx_GIVEREF(__pyx_tuple__45); /* "(tree fragment)":1 * def __pyx_unpickle_Enum(__pyx_type, long __pyx_checksum, __pyx_state): # <<<<<<<<<<<<<< * cdef object __pyx_PickleError * cdef object __pyx_result */ - __pyx_tuple__44 = PyTuple_Pack(5, __pyx_n_s_pyx_type, __pyx_n_s_pyx_checksum, __pyx_n_s_pyx_state, __pyx_n_s_pyx_PickleError, __pyx_n_s_pyx_result); if (unlikely(!__pyx_tuple__44)) __PYX_ERR(1, 1, __pyx_L1_error) - __Pyx_GOTREF(__pyx_tuple__44); - __Pyx_GIVEREF(__pyx_tuple__44); - __pyx_codeobj__45 = (PyObject*)__Pyx_PyCode_New(3, 0, 5, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__44, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_stringsource, __pyx_n_s_pyx_unpickle_Enum, 1, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__45)) __PYX_ERR(1, 1, __pyx_L1_error) + __pyx_tuple__46 = PyTuple_Pack(5, __pyx_n_s_pyx_type, __pyx_n_s_pyx_checksum, __pyx_n_s_pyx_state, __pyx_n_s_pyx_PickleError, __pyx_n_s_pyx_result); if (unlikely(!__pyx_tuple__46)) __PYX_ERR(1, 1, __pyx_L1_error) + __Pyx_GOTREF(__pyx_tuple__46); + __Pyx_GIVEREF(__pyx_tuple__46); + __pyx_codeobj__47 = (PyObject*)__Pyx_PyCode_New(3, 0, 5, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__46, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_stringsource, __pyx_n_s_pyx_unpickle_Enum, 1, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__47)) __PYX_ERR(1, 1, __pyx_L1_error) __Pyx_RefNannyFinishContext(); return 0; __pyx_L1_error:; @@ -18971,9 +19322,10 @@ if (!__Pyx_RefNanny) { __pyx_d = PyModule_GetDict(__pyx_m); if (unlikely(!__pyx_d)) __PYX_ERR(0, 1, __pyx_L1_error) Py_INCREF(__pyx_d); __pyx_b = PyImport_AddModule(__Pyx_BUILTIN_MODULE_NAME); if (unlikely(!__pyx_b)) __PYX_ERR(0, 1, __pyx_L1_error) - Py_INCREF(__pyx_b); __pyx_cython_runtime = PyImport_AddModule((char *) "cython_runtime"); if (unlikely(!__pyx_cython_runtime)) __PYX_ERR(0, 1, __pyx_L1_error) - Py_INCREF(__pyx_cython_runtime); + #if CYTHON_COMPILING_IN_PYPY + Py_INCREF(__pyx_b); + #endif if (PyObject_SetAttrString(__pyx_m, "__builtins__", __pyx_b) < 0) __PYX_ERR(0, 1, __pyx_L1_error); /*--- Initialize various global constants etc. ---*/ if (__Pyx_InitGlobals() < 0) __PYX_ERR(0, 1, __pyx_L1_error) @@ -19008,120 +19360,132 @@ if (!__Pyx_RefNanny) { if (__Pyx_patch_abc() < 0) __PYX_ERR(0, 1, __pyx_L1_error) #endif - /* "hafnian/hafnian.pyx":54 + /* "hafnian/hafnian.pyx":56 * * * def torontonian_complex(double complex[:, :] A, fsum=False): # <<<<<<<<<<<<<< * """Returns the Torontonian of a complex matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_1torontonian_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 54, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_1torontonian_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 56, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_torontonian_complex, __pyx_t_1) < 0) __PYX_ERR(0, 54, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_torontonian_complex, __pyx_t_1) < 0) __PYX_ERR(0, 56, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":86 + /* "hafnian/hafnian.pyx":88 * * * def torontonian_real(double[:, :] A, fsum=False): # <<<<<<<<<<<<<< * """Returns the Torontonian of a real matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_3torontonian_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 86, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_3torontonian_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 88, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_torontonian_real, __pyx_t_1) < 0) __PYX_ERR(0, 86, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_torontonian_real, __pyx_t_1) < 0) __PYX_ERR(0, 88, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":123 + /* "hafnian/hafnian.pyx":125 * * * def haf_rpt_real(double[:, :] A, int[:] rpt, double[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< * r"""Returns the hafnian of a real matrix A via the C++ hafnian library * using the rpt method. This method is more efficient for matrices with */ - __pyx_t_2 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(Py_None, PyBUF_WRITABLE); if (unlikely(!__pyx_t_2.memview)) __PYX_ERR(0, 123, __pyx_L1_error) + __pyx_t_2 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(Py_None, PyBUF_WRITABLE); if (unlikely(!__pyx_t_2.memview)) __PYX_ERR(0, 125, __pyx_L1_error) __pyx_k_ = __pyx_t_2; __pyx_t_2.memview = NULL; __pyx_t_2.data = NULL; - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_5haf_rpt_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 123, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_5haf_rpt_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 125, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_rpt_real, __pyx_t_1) < 0) __PYX_ERR(0, 123, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_rpt_real, __pyx_t_1) < 0) __PYX_ERR(0, 125, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":162 + /* "hafnian/hafnian.pyx":164 * * * def haf_rpt_complex(double complex[:, :] A, int[:] rpt, double complex[:] mu=None, bint loop=False): # <<<<<<<<<<<<<< * r"""Returns the hafnian of a complex matrix A via the C++ hafnian library * using the rpt method. This method is more efficient for matrices with */ - __pyx_t_3 = __Pyx_PyObject_to_MemoryviewSlice_ds___pyx_t_double_complex(Py_None, PyBUF_WRITABLE); if (unlikely(!__pyx_t_3.memview)) __PYX_ERR(0, 162, __pyx_L1_error) + __pyx_t_3 = __Pyx_PyObject_to_MemoryviewSlice_ds___pyx_t_double_complex(Py_None, PyBUF_WRITABLE); if (unlikely(!__pyx_t_3.memview)) __PYX_ERR(0, 164, __pyx_L1_error) __pyx_k__2 = __pyx_t_3; __pyx_t_3.memview = NULL; __pyx_t_3.data = NULL; - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_7haf_rpt_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 162, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_7haf_rpt_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 164, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_rpt_complex, __pyx_t_1) < 0) __PYX_ERR(0, 162, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_rpt_complex, __pyx_t_1) < 0) __PYX_ERR(0, 164, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":205 + /* "hafnian/hafnian.pyx":207 * * * def haf_int(long long[:, :] A): # <<<<<<<<<<<<<< * """Returns the hafnian of an integer matrix A via the C++ hafnian library. * Modified with permission from https://github.com/eklotek/Hafnian. */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_9haf_int, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 205, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_9haf_int, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 207, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_int, __pyx_t_1) < 0) __PYX_ERR(0, 205, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_int, __pyx_t_1) < 0) __PYX_ERR(0, 207, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":232 + /* "hafnian/hafnian.pyx":234 * * * def haf_complex(double complex[:, :] A, bint loop=False, bint recursive=True, quad=True): # <<<<<<<<<<<<<< * """Returns the hafnian of a complex matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_11haf_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 232, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_11haf_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 234, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_complex, __pyx_t_1) < 0) __PYX_ERR(0, 232, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_complex, __pyx_t_1) < 0) __PYX_ERR(0, 234, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":265 + /* "hafnian/hafnian.pyx":267 * * * def haf_real(double[:, :] A, bint loop=False, bint recursive=True, quad=True, bint approx=False, nsamples=1000): # <<<<<<<<<<<<<< * """Returns the hafnian of a real matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_13haf_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 265, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_13haf_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 267, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_real, __pyx_t_1) < 0) __PYX_ERR(0, 265, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_haf_real, __pyx_t_1) < 0) __PYX_ERR(0, 267, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":310 + /* "hafnian/hafnian.pyx":312 * * * def perm_complex(double complex[:, :] A, quad=True): # <<<<<<<<<<<<<< * """Returns the hafnian of a complex matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_15perm_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 310, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_15perm_complex, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 312, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_perm_complex, __pyx_t_1) < 0) __PYX_ERR(0, 310, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_perm_complex, __pyx_t_1) < 0) __PYX_ERR(0, 312, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "hafnian/hafnian.pyx":335 + /* "hafnian/hafnian.pyx":337 * * * def perm_real(double [:, :] A, quad=True, fsum=False): # <<<<<<<<<<<<<< * """Returns the hafnian of a real matrix A via the C++ hafnian library. * */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_17perm_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 335, __pyx_L1_error) + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_17perm_real, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 337, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_1); + if (PyDict_SetItem(__pyx_d, __pyx_n_s_perm_real, __pyx_t_1) < 0) __PYX_ERR(0, 337, __pyx_L1_error) + __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; + + /* "hafnian/hafnian.pyx":370 + * # Batch hafnian + * + * def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): # <<<<<<<<<<<<<< + * cdef int i, j, n = A.shape[0] + * cdef vector[double complex] R_mat, y_mat + */ + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6libhaf_19hermite_multidimensional, NULL, __pyx_n_s_libhaf); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 370, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_perm_real, __pyx_t_1) < 0) __PYX_ERR(0, 335, __pyx_L1_error) + if (PyDict_SetItem(__pyx_d, __pyx_n_s_hermite_multidimensional, __pyx_t_1) < 0) __PYX_ERR(0, 370, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; /* "hafnian/hafnian.pyx":1 @@ -19154,7 +19518,7 @@ if (!__Pyx_RefNanny) { * cdef strided = Enum("") # default * cdef indirect = Enum("") */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__39, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 286, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__41, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 286, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(generic); __Pyx_DECREF_SET(generic, __pyx_t_1); @@ -19168,7 +19532,7 @@ if (!__Pyx_RefNanny) { * cdef indirect = Enum("") * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__40, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 287, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__42, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 287, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(strided); __Pyx_DECREF_SET(strided, __pyx_t_1); @@ -19182,7 +19546,7 @@ if (!__Pyx_RefNanny) { * * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__41, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 288, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__43, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 288, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(indirect); __Pyx_DECREF_SET(indirect, __pyx_t_1); @@ -19196,7 +19560,7 @@ if (!__Pyx_RefNanny) { * cdef indirect_contiguous = Enum("") * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__42, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 291, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__44, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 291, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(contiguous); __Pyx_DECREF_SET(contiguous, __pyx_t_1); @@ -19210,7 +19574,7 @@ if (!__Pyx_RefNanny) { * * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__43, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 292, __pyx_L1_error) + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_MemviewEnum_type), __pyx_tuple__45, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(1, 292, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(indirect_contiguous); __Pyx_DECREF_SET(indirect_contiguous, __pyx_t_1); @@ -19520,7 +19884,11 @@ __Pyx_init_memviewslice(struct __pyx_memoryview_obj *memview, int i, retval=-1; Py_buffer *buf = &memview->view; __Pyx_RefNannySetupContext("init_memviewslice", 0); - if (memviewslice->memview || memviewslice->data) { + if (!buf) { + PyErr_SetString(PyExc_ValueError, + "buf is NULL."); + goto fail; + } else if (memviewslice->memview || memviewslice->data) { PyErr_SetString(PyExc_ValueError, "memviewslice is already initialized!"); goto fail; @@ -20457,32 +20825,6 @@ static CYTHON_INLINE PyObject *__Pyx_GetAttr3(PyObject *o, PyObject *n, PyObject return (likely(r)) ? r : __Pyx_GetAttr3Default(d); } -/* PyDictVersioning */ -#if CYTHON_USE_DICT_VERSIONS && CYTHON_USE_TYPE_SLOTS -static CYTHON_INLINE PY_UINT64_T __Pyx_get_tp_dict_version(PyObject *obj) { - PyObject *dict = Py_TYPE(obj)->tp_dict; - return likely(dict) ? __PYX_GET_DICT_VERSION(dict) : 0; -} -static CYTHON_INLINE PY_UINT64_T __Pyx_get_object_dict_version(PyObject *obj) { - PyObject **dictptr = NULL; - Py_ssize_t offset = Py_TYPE(obj)->tp_dictoffset; - if (offset) { -#if CYTHON_COMPILING_IN_CPYTHON - dictptr = (likely(offset > 0)) ? (PyObject **) ((char *)obj + offset) : _PyObject_GetDictPtr(obj); -#else - dictptr = _PyObject_GetDictPtr(obj); -#endif - } - return (dictptr && *dictptr) ? __PYX_GET_DICT_VERSION(*dictptr) : 0; -} -static CYTHON_INLINE int __Pyx_object_dict_version_matches(PyObject* obj, PY_UINT64_T tp_dict_version, PY_UINT64_T obj_dict_version) { - PyObject *dict = Py_TYPE(obj)->tp_dict; - if (unlikely(!dict) || unlikely(tp_dict_version != __PYX_GET_DICT_VERSION(dict))) - return 0; - return obj_dict_version == __Pyx_get_object_dict_version(obj); -} -#endif - /* GetModuleGlobalName */ #if CYTHON_USE_DICT_VERSIONS static PyObject *__Pyx__GetModuleGlobalName(PyObject *name, PY_UINT64_T *dict_version, PyObject **dict_cached_value) @@ -21486,13 +21828,13 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, return __pyx_t_double_complex_from_parts(a.real / b.real, a.imag / b.imag); } else { double r = b.imag / b.real; - double s = (double)(1.0) / (b.real + b.imag * r); + double s = 1.0 / (b.real + b.imag * r); return __pyx_t_double_complex_from_parts( (a.real + a.imag * r) * s, (a.imag - a.real * r) * s); } } else { double r = b.real / b.imag; - double s = (double)(1.0) / (b.imag + b.real * r); + double s = 1.0 / (b.imag + b.real * r); return __pyx_t_double_complex_from_parts( (a.real * r + a.imag) * s, (a.imag * r - a.real) * s); } @@ -21727,7 +22069,7 @@ static int __Pyx_BufFmt_ParseNumber(const char** ts) { return -1; } else { count = *t++ - '0'; - while (*t >= '0' && *t <= '9') { + while (*t >= '0' && *t < '9') { count *= 10; count += *t++ - '0'; } @@ -22301,7 +22643,7 @@ __pyx_check_suboffsets(Py_buffer *buf, int dim, CYTHON_UNUSED int ndim, int spec } } if (spec & __Pyx_MEMVIEW_PTR) { - if (!buf->suboffsets || (buf->suboffsets[dim] < 0)) { + if (!buf->suboffsets || (buf->suboffsets && buf->suboffsets[dim] < 0)) { PyErr_Format(PyExc_ValueError, "Buffer is not indirectly accessible " "in dimension %d.", dim); @@ -22555,6 +22897,28 @@ static int __Pyx_ValidateAndInit_memviewslice( return result; } +/* CIntFromPyVerify */ + #define __PYX_VERIFY_RETURN_INT(target_type, func_type, func_value)\ + __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, 0) +#define __PYX_VERIFY_RETURN_INT_EXC(target_type, func_type, func_value)\ + __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, 1) +#define __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, exc)\ + {\ + func_type value = func_value;\ + if (sizeof(target_type) < sizeof(func_type)) {\ + if (unlikely(value != (func_type) (target_type) value)) {\ + func_type zero = 0;\ + if (exc && unlikely(value == (func_type)-1 && PyErr_Occurred()))\ + return (target_type) -1;\ + if (is_unsigned && unlikely(value < zero))\ + goto raise_neg_overflow;\ + else\ + goto raise_overflow;\ + }\ + }\ + return (target_type) value;\ + } + /* CIntToPy */ static CYTHON_INLINE PyObject* __Pyx_PyInt_From_int(int value) { const int neg_one = (int) ((int) 0 - (int) 1), const_zero = (int) 0; @@ -22586,28 +22950,6 @@ static int __Pyx_ValidateAndInit_memviewslice( } } -/* CIntFromPyVerify */ - #define __PYX_VERIFY_RETURN_INT(target_type, func_type, func_value)\ - __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, 0) -#define __PYX_VERIFY_RETURN_INT_EXC(target_type, func_type, func_value)\ - __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, 1) -#define __PYX__VERIFY_RETURN_INT(target_type, func_type, func_value, exc)\ - {\ - func_type value = func_value;\ - if (sizeof(target_type) < sizeof(func_type)) {\ - if (unlikely(value != (func_type) (target_type) value)) {\ - func_type zero = 0;\ - if (exc && unlikely(value == (func_type)-1 && PyErr_Occurred()))\ - return (target_type) -1;\ - if (is_unsigned && unlikely(value < zero))\ - goto raise_neg_overflow;\ - else\ - goto raise_overflow;\ - }\ - }\ - return (target_type) value;\ - } - /* CIntToPy */ static CYTHON_INLINE PyObject* __Pyx_PyInt_From_PY_LONG_LONG(PY_LONG_LONG value) { const PY_LONG_LONG neg_one = (PY_LONG_LONG) ((PY_LONG_LONG) 0 - (PY_LONG_LONG) 1), const_zero = (PY_LONG_LONG) 0; @@ -22895,6 +23237,195 @@ __pyx_memoryview_copy_new_contig(const __Pyx_memviewslice *from_mvs, return (int) -1; } +/* CIntFromPy */ + static CYTHON_INLINE size_t __Pyx_PyInt_As_size_t(PyObject *x) { + const size_t neg_one = (size_t) ((size_t) 0 - (size_t) 1), const_zero = (size_t) 0; + const int is_unsigned = neg_one > const_zero; +#if PY_MAJOR_VERSION < 3 + if (likely(PyInt_Check(x))) { + if (sizeof(size_t) < sizeof(long)) { + __PYX_VERIFY_RETURN_INT(size_t, long, PyInt_AS_LONG(x)) + } else { + long val = PyInt_AS_LONG(x); + if (is_unsigned && unlikely(val < 0)) { + goto raise_neg_overflow; + } + return (size_t) val; + } + } else +#endif + if (likely(PyLong_Check(x))) { + if (is_unsigned) { +#if CYTHON_USE_PYLONG_INTERNALS + const digit* digits = ((PyLongObject*)x)->ob_digit; + switch (Py_SIZE(x)) { + case 0: return (size_t) 0; + case 1: __PYX_VERIFY_RETURN_INT(size_t, digit, digits[0]) + case 2: + if (8 * sizeof(size_t) > 1 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 2 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) >= 2 * PyLong_SHIFT) { + return (size_t) (((((size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0])); + } + } + break; + case 3: + if (8 * sizeof(size_t) > 2 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 3 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((((unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) >= 3 * PyLong_SHIFT) { + return (size_t) (((((((size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0])); + } + } + break; + case 4: + if (8 * sizeof(size_t) > 3 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 4 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((((((unsigned long)digits[3]) << PyLong_SHIFT) | (unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) >= 4 * PyLong_SHIFT) { + return (size_t) (((((((((size_t)digits[3]) << PyLong_SHIFT) | (size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0])); + } + } + break; + } +#endif +#if CYTHON_COMPILING_IN_CPYTHON + if (unlikely(Py_SIZE(x) < 0)) { + goto raise_neg_overflow; + } +#else + { + int result = PyObject_RichCompareBool(x, Py_False, Py_LT); + if (unlikely(result < 0)) + return (size_t) -1; + if (unlikely(result == 1)) + goto raise_neg_overflow; + } +#endif + if (sizeof(size_t) <= sizeof(unsigned long)) { + __PYX_VERIFY_RETURN_INT_EXC(size_t, unsigned long, PyLong_AsUnsignedLong(x)) +#ifdef HAVE_LONG_LONG + } else if (sizeof(size_t) <= sizeof(unsigned PY_LONG_LONG)) { + __PYX_VERIFY_RETURN_INT_EXC(size_t, unsigned PY_LONG_LONG, PyLong_AsUnsignedLongLong(x)) +#endif + } + } else { +#if CYTHON_USE_PYLONG_INTERNALS + const digit* digits = ((PyLongObject*)x)->ob_digit; + switch (Py_SIZE(x)) { + case 0: return (size_t) 0; + case -1: __PYX_VERIFY_RETURN_INT(size_t, sdigit, (sdigit) (-(sdigit)digits[0])) + case 1: __PYX_VERIFY_RETURN_INT(size_t, digit, +digits[0]) + case -2: + if (8 * sizeof(size_t) - 1 > 1 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 2 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, long, -(long) (((((unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 2 * PyLong_SHIFT) { + return (size_t) (((size_t)-1)*(((((size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + case 2: + if (8 * sizeof(size_t) > 1 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 2 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 2 * PyLong_SHIFT) { + return (size_t) ((((((size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + case -3: + if (8 * sizeof(size_t) - 1 > 2 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 3 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, long, -(long) (((((((unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 3 * PyLong_SHIFT) { + return (size_t) (((size_t)-1)*(((((((size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + case 3: + if (8 * sizeof(size_t) > 2 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 3 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((((unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 3 * PyLong_SHIFT) { + return (size_t) ((((((((size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + case -4: + if (8 * sizeof(size_t) - 1 > 3 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 4 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, long, -(long) (((((((((unsigned long)digits[3]) << PyLong_SHIFT) | (unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 4 * PyLong_SHIFT) { + return (size_t) (((size_t)-1)*(((((((((size_t)digits[3]) << PyLong_SHIFT) | (size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + case 4: + if (8 * sizeof(size_t) > 3 * PyLong_SHIFT) { + if (8 * sizeof(unsigned long) > 4 * PyLong_SHIFT) { + __PYX_VERIFY_RETURN_INT(size_t, unsigned long, (((((((((unsigned long)digits[3]) << PyLong_SHIFT) | (unsigned long)digits[2]) << PyLong_SHIFT) | (unsigned long)digits[1]) << PyLong_SHIFT) | (unsigned long)digits[0]))) + } else if (8 * sizeof(size_t) - 1 > 4 * PyLong_SHIFT) { + return (size_t) ((((((((((size_t)digits[3]) << PyLong_SHIFT) | (size_t)digits[2]) << PyLong_SHIFT) | (size_t)digits[1]) << PyLong_SHIFT) | (size_t)digits[0]))); + } + } + break; + } +#endif + if (sizeof(size_t) <= sizeof(long)) { + __PYX_VERIFY_RETURN_INT_EXC(size_t, long, PyLong_AsLong(x)) +#ifdef HAVE_LONG_LONG + } else if (sizeof(size_t) <= sizeof(PY_LONG_LONG)) { + __PYX_VERIFY_RETURN_INT_EXC(size_t, PY_LONG_LONG, PyLong_AsLongLong(x)) +#endif + } + } + { +#if CYTHON_COMPILING_IN_PYPY && !defined(_PyLong_AsByteArray) + PyErr_SetString(PyExc_RuntimeError, + "_PyLong_AsByteArray() not available in PyPy, cannot convert large numbers"); +#else + size_t val; + PyObject *v = __Pyx_PyNumber_IntOrLong(x); + #if PY_MAJOR_VERSION < 3 + if (likely(v) && !PyLong_Check(v)) { + PyObject *tmp = v; + v = PyNumber_Long(tmp); + Py_DECREF(tmp); + } + #endif + if (likely(v)) { + int one = 1; int is_little = (int)*(unsigned char *)&one; + unsigned char *bytes = (unsigned char *)&val; + int ret = _PyLong_AsByteArray((PyLongObject *)v, + bytes, sizeof(val), + is_little, !is_unsigned); + Py_DECREF(v); + if (likely(!ret)) + return val; + } +#endif + return (size_t) -1; + } + } else { + size_t val; + PyObject *tmp = __Pyx_PyNumber_IntOrLong(x); + if (!tmp) return (size_t) -1; + val = __Pyx_PyInt_As_size_t(tmp); + Py_DECREF(tmp); + return val; + } +raise_overflow: + PyErr_SetString(PyExc_OverflowError, + "value too large to convert to size_t"); + return (size_t) -1; +raise_neg_overflow: + PyErr_SetString(PyExc_OverflowError, + "can't convert negative value to size_t"); + return (size_t) -1; +} + /* CIntFromPy */ static CYTHON_INLINE long __Pyx_PyInt_As_long(PyObject *x) { const long neg_one = (long) ((long) 0 - (long) 1), const_zero = (long) 0; diff --git a/hafnian/hafnian.pyx b/hafnian/hafnian.pyx index 6006ae123..32ceb9465 100644 --- a/hafnian/hafnian.pyx +++ b/hafnian/hafnian.pyx @@ -46,6 +46,8 @@ cdef extern from "../src/hafnian.hpp" namespace "hafnian": double complex torontonian_quad(vector[double complex] &mat) double torontonian_fsum[T](vector[T] &mat) + vector[double complex] hermite_multidimensional_cpp(vector[double complex] &mat, vector[double complex] &d, int &resolution, bint &renorm) + # ============================================================================== # Torontonian @@ -361,3 +363,25 @@ def perm_real(double [:, :] A, quad=True, fsum=False): return permanent(mat) + +# ============================================================================== +# Batch hafnian + +def hermite_multidimensional(double complex[:, :] A, double complex[:] d, int resolution, ren=False): + cdef int i, j, n = A.shape[0] + cdef vector[double complex] R_mat, y_mat + + cdef int renorm = 0 + + if ren: + renorm = 1 + + for i in range(n): + for j in range(n): + R_mat.push_back(A[i, j]) + + + for i in range(n): + y_mat.push_back(d[i]) + + return hermite_multidimensional_cpp(R_mat, y_mat, resolution, renorm) diff --git a/hafnian/quantum.py b/hafnian/quantum.py index 7729b708e..a38488658 100644 --- a/hafnian/quantum.py +++ b/hafnian/quantum.py @@ -89,6 +89,8 @@ from scipy.special import factorial as fac from ._hafnian import hafnian, hafnian_repeated, reduction +from ._hermite_multidimensional import hermite_multidimensional, hafnian_batched + def reduced_gaussian(mu, cov, modes): @@ -311,7 +313,6 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba complex: the density matrix element """ rpt = i + j - beta = Beta(mu, hbar=hbar) A = Amat(cov, hbar=hbar) if np.linalg.norm(beta) < tol: @@ -325,7 +326,6 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba # replace the diagonal of A with gamma # gamma = X @ np.linalg.inv(Q).conj() @ beta gamma = beta.conj() - A @ beta - if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: A_rpt = reduction(A, rpt) np.fill_diagonal(A_rpt, reduction(gamma, rpt)) @@ -339,6 +339,81 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba return haf / np.sqrt(np.prod(fac(rpt))) +def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2): + r"""Returns the density matrix of a (PNR post-selected) Gaussian state. + + The resulting density matrix will have shape + + .. math:: \underbrace{D\times D \times \cdots \times D}_{2M} + + where :math:`D` is the Fock space cutoff, and :math:`M` is the + number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. + + Note that we use the Strawberry Fields convention for indexing the density + matrix; the first two dimensions correspond to subsystem 1, the second two + dimensions correspond to subsystem 2, etc. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + post_select (dict): dictionary containing the post-selected modes, of + the form ``{mode: value}``. If post_select is None the whole non post-selected density matrix + is calculated directly using (multidimensional) Hermite polynomials, which is significantly faster + than calculating one hafnian at a time. + normalize (bool): If ``True``, a post-selected density matrix is re-normalized. + cutoff (dim): the final length (i.e., Hilbert space dimension) of each + mode in the density matrix. + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + np.array[complex]: the density matrix of the Gaussian state + """ + N = len(mu) // 2 + pref = prefactor(mu, cov, hbar=hbar) + + if post_select is None: + A = Amat(cov, hbar=hbar) + if np.allclose(mu, np.zeros_like(mu)): + return pref * hermite_multidimensional(-A, cutoff, renorm=True) + try: + y = np.linalg.inv(A) @ np.linalg.inv(Qmat(cov)) @ Beta(mu, hbar=hbar) + return pref * hermite_multidimensional(-A, cutoff, y=-y, renorm=True) + except np.linalg.LinAlgError: + pass + post_select = {} + + M = N - len(post_select) + rho = np.zeros([cutoff] * (2 * M), dtype=np.complex128) + + for idx in product(range(cutoff), repeat=2 * M): + el = [] + + counter = count(0) + modes = (np.arange(2 * N) % N).tolist() + el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] + + el = np.array(el).reshape(2, -1) + el0 = el[0].tolist() + el1 = el[1].tolist() + + sf_idx = np.array(idx).reshape(2, -1) + sf_el = tuple(sf_idx[::-1].T.flatten()) + + rho[sf_el] = density_matrix_element(mu, cov, el0, el1, include_prefactor=False, hbar=hbar) + + rho *= pref + + if normalize: + # construct the standard 2D density matrix, and take the trace + new_ax = np.arange(2 * M).reshape([M, 2]).T.flatten() + tr = np.trace(rho.transpose(new_ax).reshape([cutoff ** M, cutoff ** M])).real + # renormalize + rho /= tr + + return rho + + def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2, check_purity=True): r"""Returns the :math:`\langle i | \psi\rangle` element of the state ket of a Gaussian state defined by covariance matrix cov. @@ -428,100 +503,44 @@ def state_vector(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2, c if not is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): raise ValueError("The covariance matrix does not correspond to a pure state") - if post_select is None: - post_select = {} - beta = Beta(mu, hbar=hbar) A = Amat(cov, hbar=hbar) + Q = Qmat(cov, hbar=hbar) (n, _) = cov.shape N = n // 2 B = A[0:N, 0:N] alpha = beta[0:N] - M = N - len(post_select) - psi = np.zeros([cutoff] * (M), dtype=np.complex128) - - for idx in product(range(cutoff), repeat=M): - el = [] - - counter = count(0) - modes = (np.arange(N)).tolist() - el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] - psi[idx] = pure_state_amplitude(mu, cov, el, check_purity=False, include_prefactor=False) - pref = np.exp(-0.5 * (np.linalg.norm(alpha) ** 2 - alpha.conj() @ B @ alpha.conj())) - psi = psi * pref - - if normalize: - norm = np.sqrt(np.sum(np.abs(psi) ** 2)) - psi = psi / norm - - return psi - - -def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2): - r"""Returns the density matrix of a (PNR post-selected) Gaussian state. - - The resulting density matrix will have shape - - .. math:: \underbrace{D\times D \times \cdots \times D}_{2M} - - where :math:`D` is the Fock space cutoff, and :math:`M` is the - number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. - - Note that we use the Strawberry Fields convention for indexing the density - matrix; the first two dimensions correspond to subsystem 1, the second two - dimensions correspond to subsystem 2, etc. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - post_select (dict): dictionary containing the post-selected modes, of - the form ``{mode: value}``. - normalize (bool): If ``True``, a post-selected density matrix is re-normalized. - cutoff (dim): the final length (i.e., Hilbert space dimension) of each - mode in the density matrix. - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - Returns: - np.array[complex]: the density matrix of the Gaussian state - """ if post_select is None: - post_select = {} - - N = len(mu) // 2 - M = N - len(post_select) - - rho = np.zeros([cutoff] * (2 * M), dtype=np.complex128) - - for idx in product(range(cutoff), repeat=2 * M): - el = [] - - counter = count(0) - modes = (np.arange(2 * N) % N).tolist() - el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] + psi = ( + pref + * hafnian_batched(B, cutoff, mu=alpha, renorm=True) + / np.sqrt(np.sqrt(np.linalg.det(Q).real)) + ) + else: + M = N - len(post_select) + psi = np.zeros([cutoff] * (M), dtype=np.complex128) - el = np.array(el).reshape(2, -1) - el0 = el[0].tolist() - el1 = el[1].tolist() + for idx in product(range(cutoff), repeat=M): + el = [] - sf_idx = np.array(idx).reshape(2, -1) - sf_el = tuple(sf_idx[::-1].T.flatten()) - - rho[sf_el] = density_matrix_element(mu, cov, el0, el1, include_prefactor=False, hbar=hbar) + counter = count(0) + modes = (np.arange(N)).tolist() + el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] + psi[idx] = pure_state_amplitude( + mu, cov, el, check_purity=False, include_prefactor=False + ) - rho *= prefactor(mu, cov, hbar=hbar) + psi = psi * pref if normalize: - # construct the standard 2D density matrix, and take the trace - new_ax = np.arange(2 * M).reshape([M, 2]).T.flatten() - tr = np.trace(rho.transpose(new_ax).reshape([cutoff ** M, cutoff ** M])).real - # renormalize - rho /= tr + norm = np.sqrt(np.sum(np.abs(psi) ** 2)) + psi = psi / norm - return rho + return psi def find_scaling_adjacency_matrix(A, n_mean): diff --git a/hafnian/tests/test_hafnian_approx.py b/hafnian/tests/test_hafnian_approx.py index f37d5ac43..7e608cf83 100644 --- a/hafnian/tests/test_hafnian_approx.py +++ b/hafnian/tests/test_hafnian_approx.py @@ -30,7 +30,7 @@ def test_rank_one(n): A = np.outer(x, x) exact = factorial2(n - 1) * np.prod(x) approx = haf_real(A, approx=True, nsamples=10000) - assert np.allclose(approx, exact, rtol=1e-1, atol=0) + assert np.allclose(approx, exact, rtol=2e-1, atol=0) def test_approx_complex_error(): diff --git a/hafnian/tests/test_hermite_multidimensional.py b/hafnian/tests/test_hermite_multidimensional.py new file mode 100644 index 000000000..cddfc80f1 --- /dev/null +++ b/hafnian/tests/test_hermite_multidimensional.py @@ -0,0 +1,123 @@ +# Copyright 2019 Xanadu Quantum Technologies Inc. + +# 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. +"""Tests for the batch hafnian wrapper function""" +# pylint: disable=no-self-use,redefined-outer-name +from itertools import product + +import numpy as np + +from scipy.special import eval_hermitenorm, eval_hermite + +from hafnian import hermite_multidimensional, hafnian_batched, hafnian_repeated + + +def test_hermite_multidimensional_renorm(): + """ This tests the renormalized batchhafnian wrapper function to compute photon number statistics for a fixed gaussian state. + """ + B = np.sqrt(0.5) * np.array([[0, 1], [1, 0]]) + 0 * 1j + res = 10 + expected = np.diag(0.5 ** (np.arange(0, res) / 2)) + array = hermite_multidimensional(-B, res, renorm=True) + + assert np.allclose(array, expected) + + +def test_reduction_to_physicists_polys(): + """Tests that the multidimensional hermite polynomials reduce to the regular physicists' hermite polynomials in the appropriate limit""" + x = np.arange(-1, 1, 0.1) + init = 1 + n_max = 5 + A = np.ones([init, init], dtype=complex) + vals = np.array( + [hermite_multidimensional(2 * A, n_max, y=np.array([x0], dtype=complex)) for x0 in x] + ).T + expected = np.array([eval_hermite(i, x) for i in range(len(vals))]) + assert np.allclose(vals, expected) + + +def test_reduction_to_probabilist_polys(): + """Tests that the multidimensional hermite polynomials reduce to the regular probabilist' hermite polynomials in the appropriate limit""" + x = np.arange(-1, 1, 0.1) + init = 1 + n_max = 5 + A = np.ones([init, init], dtype=complex) + vals = np.array( + [hermite_multidimensional(A, n_max, y=np.array([x0], dtype=complex)) for x0 in x] + ).T + expected = np.array([eval_hermitenorm(i, x) for i in range(len(vals))]) + assert np.allclose(vals, expected) + + +def test_hafnian_batched(): + """Test hafnian_batched against hafnian_repeated for a random symmetric matrix""" + n_modes = 4 + A = np.random.rand(n_modes, n_modes) + 1j * np.random.rand(n_modes, n_modes) + A += A.T + n_photon = 5 + v1 = np.array([hafnian_repeated(A, q) for q in product(np.arange(n_photon), repeat=n_modes)]) + assert np.allclose(hafnian_batched(A, n_photon, make_tensor=False), v1) + + +def test_hafnian_batched_loops(): + """Test hafnian_batched with loops against hafnian_repeated with loops for a random symmetric matrix + and a random vector of loops + """ + n_modes = 4 + A = np.random.rand(n_modes, n_modes) + 1j * np.random.rand(n_modes, n_modes) + A += A.T + mu = np.random.rand(n_modes) + 1j * np.random.rand(n_modes) + n_photon = 5 + v1 = np.array( + [ + hafnian_repeated(A, q, mu=mu, loop=True) + for q in product(np.arange(n_photon), repeat=n_modes) + ] + ) + expected = hafnian_batched(A, n_photon, mu=mu, make_tensor=False) + + assert np.allclose(expected, v1) + + +def test_hafnian_batched_loops_no_edges(): + """Test hafnian_batched with loops against hafnian_repeated with loops for a random symmetric matrix + and a random vector of loops + """ + n_modes = 4 + A = np.zeros([n_modes, n_modes], dtype=complex) + mu = np.random.rand(n_modes) + 1j * np.random.rand(n_modes) + n_photon = 5 + v1 = np.array( + [ + hafnian_repeated(A, q, mu=mu, loop=True) + for q in product(np.arange(n_photon), repeat=n_modes) + ] + ) + expected = hafnian_batched(A, n_photon, mu=mu, make_tensor=False) + + assert np.allclose(expected, v1) + + +def test_hafnian_batched_zero_loops_no_edges(): + """Test hafnian_batched with loops against hafnian_repeated with loops for a the zero matrix + and a loops + """ + n_modes = 4 + A = np.zeros([n_modes, n_modes], dtype=complex) + n_photon = 5 + v1 = np.array( + [hafnian_repeated(A, q, loop=True) for q in product(np.arange(n_photon), repeat=n_modes)] + ) + expected = hafnian_batched(A, n_photon, make_tensor=False) + + assert np.allclose(expected, v1) diff --git a/hafnian/tests/test_quantum.py b/hafnian/tests/test_quantum.py index 345e92221..2d9f53730 100644 --- a/hafnian/tests/test_quantum.py +++ b/hafnian/tests/test_quantum.py @@ -354,8 +354,8 @@ def test_density_matrix_squeezed(): assert np.allclose(res, expected) -def test_density_matrix_displaced_squeezed(): - """Test density matrix for a squeezed state""" +def test_coherent_squeezed(): + """Test density matrix for a squeezed displaced state""" r = 0.43 mu = np.array([0.24, -0.2]) diff --git a/hafnian/tests/test_samples.py b/hafnian/tests/test_samples.py index b617498a8..f75063cd4 100644 --- a/hafnian/tests/test_samples.py +++ b/hafnian/tests/test_samples.py @@ -225,7 +225,7 @@ def test_hafnian_sample_graph(self): mean_n = 0.5 samples = hafnian_sample_graph(A, mean_n, samples=n_samples) approx_mean_n = np.sum(samples) / n_samples - assert np.allclose(mean_n, approx_mean_n, rtol=1e-1) + assert np.allclose(mean_n, approx_mean_n, rtol=2e-1) def test_single_pm_graphs(self): """Tests that the number of photons is the same for modes i and n-i @@ -259,8 +259,6 @@ def test_probability_vacuum(self): A = np.triu(A) A = A + np.transpose(A) - for i in range(n): - A[i, i] = 0 n_mean = 1.0 Q = gen_Qmat_from_graph(A, n_mean) prob0 = np.abs(1 / (np.sqrt(np.linalg.det(Q)))) diff --git a/setup.py b/setup.py index 7bf9b79a7..5b3bd0a7b 100644 --- a/setup.py +++ b/setup.py @@ -112,6 +112,7 @@ def cythonize(x, compile_time_env=None): "src/hafnian_approx.hpp", "src/torontonian.hpp", "src/permanent.hpp", + "src/hermite_multidimensional.hpp", "src/stdafx.h", "src/fsum.hpp"], include_dirs=C_INCLUDE_PATH, diff --git a/src/hafnian.hpp b/src/hafnian.hpp index 8bdac0d6c..a0eaa23b9 100644 --- a/src/hafnian.hpp +++ b/src/hafnian.hpp @@ -19,6 +19,7 @@ #include #include #include +#include /** * @namespace hafnian diff --git a/src/hermite_multidimensional.hpp b/src/hermite_multidimensional.hpp new file mode 100644 index 000000000..8b69e4c22 --- /dev/null +++ b/src/hermite_multidimensional.hpp @@ -0,0 +1,205 @@ +// Copyright 2019 Xanadu Quantum Technologies Inc. + + +// 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. + +#pragma once +#include + +#include +#include + +typedef unsigned long long int ullint; + + +ullint vec2index(std::vector &pos, int resolution) { + int dim = pos.size(); + ullint nextCoordinate = 0; + + nextCoordinate = pos[0]-1; + for(int ii = 0; ii < dim-1; ii++) { + nextCoordinate = nextCoordinate*resolution + (pos[ii+1]-1); + } + + return nextCoordinate; + +} + +std::vector find_rep(int val, int base, int n) { + std::vector x(n, 0); + int local_val = val; + + x[0] = 1; + + for (int i = 1; i < n; i++) + x[i] = x[i-1]*base; + + std::vector digits(n, 0); + + for (int i = 0; i < n; i++) { + digits[i] = local_val/x[n-i-1]; + local_val = local_val - digits[i] * x[n-i-1]; + } + + return digits; + +} + + +long double sqrtfactorial(int nn) +{ + long double n = static_cast(nn); + + if(n > 1) + return std::sqrt(n) * sqrtfactorial(n - 1); + else + return 1; +} + +/** + * Renormalizes and unnormalized photon number statistics of a Gaussian state. + * Based on the MATLAB code available at: https://github.com/clementsw/gaussian-optics + * + * @param tn unnormalized flattened vector of size \f$res**nmodes$ representing unnormalized photon number statistics + * \f$2n\times 2n\f$ row-ordered symmetric matrix. + * @param nmodes number of modes + * @param res highest number of photons to be resolved. + * + * $return Renormalized photon number statistics + */ +template +inline std::vector renormalization(std::vector tn, int nmodes, int res) { + std::vector invsqfacts(res, 0); + std::vector digits(nmodes, 0); + + ullint Hdim = pow(res, nmodes); + + for (int i = 0; i < res; i++) + invsqfacts[i] = sqrtfactorial(i); + + for (ullint i = 0; i < Hdim; i++) { + digits = find_rep(i, res, nmodes); + long double pref = 1; + for (int j = 0; j < nmodes; j++) + pref *= 1.0L/invsqfacts[digits[j]]; + tn[i] = tn[i]*static_cast(pref); + } + + return tn; + +} + + + +namespace hafnian { + +/** + * Returns photon number statistics of a Gaussian state for a given covariance matrix `mat`. + * Based on the MATLAB code available at: https://github.com/clementsw/gaussian-optics + * + * @param mat a flattened vector of size \f$2n^2\f$, representing an + * \f$2n\times 2n\f$ row-ordered symmetric matrix. + * @param d a flattened vector of size \f$2n\f$, representing the first order moments. + * @param resolution highest number of photons to be resolved. + * + */ +template +inline std::vector hermite_multidimensional_cpp(std::vector &R_mat, std::vector &y_mat, int &resolution, int &renorm) { + int dim = std::sqrt(static_cast(R_mat.size())); + + namespace eg = Eigen; + + eg::Matrix R = eg::Map, eg::Unaligned>(R_mat.data(), dim, dim); + eg::Matrix y = eg::Map, eg::Unaligned>(y_mat.data(), dim, dim); + + ullint Hdim = pow(resolution, dim); + std::vector H(Hdim, 0); + std::vector ren_factor(Hdim, 0); + + H[0] = 1; + ren_factor[0] = 1; + + std::vector nextPos(dim, 1); + std::vector jumpFrom(dim, 1); + std::vector ek(dim, 0); + std::vector factors(resolution+1, 0); + int jump = 0; + + + for (ullint jj = 0; jj < Hdim-1; jj++) { + + if (jump > 0) { + jumpFrom[jump] += 1; + jump = 0; + } + + + for (int ii = 0; ii < dim; ii++) { + std::vector forwardStep(dim, 0); + forwardStep[ii] = 1; + + if ( forwardStep[ii] + nextPos[ii] > resolution) { + nextPos[ii] = 1; + jumpFrom[ii] = 1; + jump = ii+1; + } + else { + jumpFrom[ii] = nextPos[ii]; + nextPos[ii] = nextPos[ii] + 1; + break; + } + } + + for (int ii = 0; ii < dim; ii++) + ek[ii] = nextPos[ii] - jumpFrom[ii]; + + int k = 0; + for(; k < static_cast(ek.size()); k++) { + if(ek[k]) break; + } + + ullint nextCoordinate = vec2index(nextPos, resolution); + ullint fromCoordinate = vec2index(jumpFrom, resolution); + + + for (int ii = 0; ii < dim; ii++) { + H[nextCoordinate] = H[nextCoordinate] + R(k, ii) * y(ii, 0); + } + H[nextCoordinate] = H[nextCoordinate] * H[fromCoordinate]; + + std::vector tmpjump(dim, 0); + + for (int ii = 0; ii < dim; ii++) { + if (jumpFrom[ii] > 1) { + std::vector prevJump(dim, 0); + prevJump[ii] = 1; + std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus()); + ullint prevCoordinate = vec2index(tmpjump, resolution); + H[nextCoordinate] = H[nextCoordinate] - (static_cast(jumpFrom[ii]-1))*static_cast(R(k,ii))*H[prevCoordinate]; + + } + } + + } + + if (renorm) { + H = renormalization(H, dim, resolution); + } + + return H; + +} + + + +} diff --git a/src/permanent.hpp b/src/permanent.hpp index f9dcd6777..6a463bd4f 100644 --- a/src/permanent.hpp +++ b/src/permanent.hpp @@ -22,10 +22,10 @@ typedef long long int llint; //typedef long double qp; #if (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__) || defined(__INTEL_COMPILER)) - typedef __float128 qp; - //#include +typedef __float128 qp; +//#include #else - typedef long double qp; +typedef long double qp; #endif /** @@ -210,7 +210,7 @@ inline double perm_fsum(std::vector &mat) { #endif std::vector tot(nthreads, static_cast(0)); - + std::vector threadbound_low(nthreads); std::vector threadbound_hi(nthreads); @@ -224,7 +224,7 @@ inline double perm_fsum(std::vector &mat) { #pragma omp parallel for shared(tot) for (int ii = 0; ii < nthreads; ii++) { - fsum::sc_partials permtmp; // = 0; + fsum::sc_partials permtmp; // = 0; int cntr = 0; std::vector chitmp(n, 0); @@ -272,9 +272,8 @@ inline double perm_fsum(std::vector &mat) { permtmp -= rowsumprod; } - tot[ii] = permtmp; + tot[ii] = permtmp; } - return static_cast(std::accumulate(tot.begin(), tot.end(), static_cast(0))); } diff --git a/src/tests/Makefile b/src/tests/Makefile index 9ef5a13df..07786f9ab 100644 --- a/src/tests/Makefile +++ b/src/tests/Makefile @@ -18,7 +18,8 @@ EIGEN_INCLUDE_DIR=/usr/include/eigen3 CC=g++ CFLAGS=-std=c++11 -O3 -g -Wall -fopenmp -fPIC -I/usr/include \ -I.. -I$(EIGEN_INCLUDE_DIR) -I$(GOOGLETEST_DIR)/include \ - -march=native + -march=native + LFLAGS=-fopenmp -L$(GOOGLETEST_DIR)/lib -lgtest all: test diff --git a/src/tests/hafnian_unittest.cpp b/src/tests/hafnian_unittest.cpp index ae94a2afb..74784654b 100644 --- a/src/tests/hafnian_unittest.cpp +++ b/src/tests/hafnian_unittest.cpp @@ -19,6 +19,8 @@ #include const double tol = 1.0e-10f; +const double tol2 = 1.0e-7f; + namespace permanent { TEST(PermanentRealFsum, CompleteGraph) { @@ -41,14 +43,14 @@ TEST(PermanentFsum, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 3; i++) { - for(int j = 0; j < 3; j++) { + for (int j = 0; j < 3; j++) { double randnum = distribution(generator); - mat[i*3+j] = randnum; + mat[i * 3 + j] = randnum; } } - double expected = mat[2]*mat[4]*mat[6] + mat[1]*mat[5]*mat[6] + mat[2]*mat[3]*mat[7] - + mat[0]*mat[5]*mat[7] + mat[1]*mat[3]*mat[8] + mat[0]*mat[4]*mat[8]; + double expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + mat[2] * mat[3] * mat[7] + + mat[0] * mat[5] * mat[7] + mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; EXPECT_NEAR(expected, hafnian::permanent_fsum(mat), tol); @@ -74,14 +76,14 @@ TEST(PermanentReal, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 3; i++) { - for(int j = 0; j < 3; j++) { + for (int j = 0; j < 3; j++) { double randnum = distribution(generator); - mat[i*3+j] = randnum; + mat[i * 3 + j] = randnum; } } - double expected = mat[2]*mat[4]*mat[6] + mat[1]*mat[5]*mat[6] + mat[2]*mat[3]*mat[7] - + mat[0]*mat[5]*mat[7] + mat[1]*mat[3]*mat[8] + mat[0]*mat[4]*mat[8]; + double expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + mat[2] * mat[3] * mat[7] + + mat[0] * mat[5] * mat[7] + mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; EXPECT_NEAR(expected, hafnian::permanent_quad(mat), tol); @@ -97,15 +99,15 @@ TEST(PermanentComplex, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 3; i++) { - for(int j = 0; j < 3; j++) { + for (int j = 0; j < 3; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat[i*3+j] = std::complex(randnum1, randnum2); + mat[i * 3 + j] = std::complex(randnum1, randnum2); } } - std::complex expected = mat[2]*mat[4]*mat[6] + mat[1]*mat[5]*mat[6] + mat[2]*mat[3]*mat[7] - + mat[0]*mat[5]*mat[7] + mat[1]*mat[3]*mat[8] + mat[0]*mat[4]*mat[8]; + std::complex expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + mat[2] * mat[3] * mat[7] + + mat[0] * mat[5] * mat[7] + mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; std::complex perm = hafnian::permanent_quad(mat); @@ -139,13 +141,13 @@ TEST(HafianRecursiveDouble, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum = distribution(generator); - mat[i*4+j] = randnum; + mat[i * 4 + j] = randnum; } } - double expected = mat[1]*mat[11] + mat[2]*mat[7] + mat[3]*mat[6]; + double expected = mat[1] * mat[11] + mat[2] * mat[7] + mat[3] * mat[6]; EXPECT_NEAR(expected, hafnian::hafnian_recursive_quad(mat), tol); } @@ -220,14 +222,14 @@ TEST(HafianRecursiveDoubleComplex, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat[i*4+j] = std::complex(randnum1, randnum2); + mat[i * 4 + j] = std::complex(randnum1, randnum2); } } - std::complex expected = mat[1]*mat[11] + mat[2]*mat[7] + mat[3]*mat[6]; + std::complex expected = mat[1] * mat[11] + mat[2] * mat[7] + mat[3] * mat[6]; std::complex haf = hafnian::hafnian_recursive_quad(mat); @@ -294,13 +296,13 @@ TEST(HafianEigenDouble, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum = distribution(generator); - mat[i*4+j] = randnum; + mat[i * 4 + j] = randnum; } } - double expected = mat[1]*mat[11] + mat[2]*mat[7] + mat[3]*mat[6]; + double expected = mat[1] * mat[11] + mat[2] * mat[7] + mat[3] * mat[6]; EXPECT_NEAR(expected, hafnian::hafnian_eigen(mat), tol); } @@ -375,14 +377,14 @@ TEST(HafianEigenDoubleComplex, Random) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat[i*4+j] = std::complex(randnum1, randnum2); + mat[i * 4 + j] = std::complex(randnum1, randnum2); } } - std::complex expected = mat[1]*mat[11] + mat[2]*mat[7] + mat[3]*mat[6]; + std::complex expected = mat[1] * mat[11] + mat[2] * mat[7] + mat[3] * mat[6]; std::complex haf = hafnian::hafnian_eigen(mat); @@ -451,8 +453,8 @@ TEST(HafnianApproxNonngeative, Random) { } for (int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - mat4[i*n+j] = x4[i]*x4[j]; + for (int j = 0; j < n; j++) { + mat4[i * n + j] = x4[i] * x4[j]; } } @@ -466,8 +468,8 @@ TEST(HafnianApproxNonngeative, Random) { } for (int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - mat6[i*n+j] = x6[i]*x6[j]; + for (int j = 0; j < n; j++) { + mat6[i * n + j] = x6[i] * x6[j]; } } @@ -480,8 +482,8 @@ TEST(HafnianApproxNonngeative, Random) { } for (int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - mat8[i*n+j] = x8[i]*x8[j]; + for (int j = 0; j < n; j++) { + mat8[i * n + j] = x8[i] * x8[j]; } } @@ -493,9 +495,9 @@ TEST(HafnianApproxNonngeative, Random) { double haf6 = hafnian::hafnian_approx(mat6, nsamples); double haf8 = hafnian::hafnian_approx(mat8, nsamples); - EXPECT_NEAR(expected4, haf4, haf4/15.0); - EXPECT_NEAR(expected6, haf6, haf6/15.0); - EXPECT_NEAR(expected8, haf8, haf8/15.0); + EXPECT_NEAR(expected4, haf4, haf4 / 15.0); + EXPECT_NEAR(expected6, haf6, haf6 / 15.0); + EXPECT_NEAR(expected8, haf8, haf8 / 15.0); } @@ -530,16 +532,16 @@ TEST(HafnianRepeatedDouble, AllOneRpt) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum = distribution(generator); - mat2rand[i*2+j] = randnum; + mat2rand[i * 2 + j] = randnum; } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum = distribution(generator); - mat4rand[i*4+j] = randnum; + mat4rand[i * 4 + j] = randnum; } } @@ -548,7 +550,7 @@ TEST(HafnianRepeatedDouble, AllOneRpt) { double expected2rand = mat2rand[1]; double expected4 = 3; - double expected4rand = mat4rand[1]*mat4rand[11] + mat4rand[2]*mat4rand[7] + mat4rand[3]*mat4rand[6];; + double expected4rand = mat4rand[1] * mat4rand[11] + mat4rand[2] * mat4rand[7] + mat4rand[3] * mat4rand[6];; double haf2 = hafnian::hafnian_rpt_quad(mat2, rpt2); double haf2rand = hafnian::hafnian_rpt_quad(mat2rand, rpt2); @@ -592,18 +594,18 @@ TEST(HafnianRepeatedComplex, AllOneRpt) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat2rand[i*2+j] = std::complex(randnum1, randnum2); + mat2rand[i * 2 + j] = std::complex(randnum1, randnum2); } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat4rand[i*4+j] = std::complex(randnum1, randnum2); + mat4rand[i * 4 + j] = std::complex(randnum1, randnum2); } } @@ -615,7 +617,7 @@ TEST(HafnianRepeatedComplex, AllOneRpt) { double expected4_re = 3; double expected4_im = 0; - std::complex expected4rand = mat4rand[1]*mat4rand[11] + mat4rand[2]*mat4rand[7] + mat4rand[3]*mat4rand[6];; + std::complex expected4rand = mat4rand[1] * mat4rand[11] + mat4rand[2] * mat4rand[7] + mat4rand[3] * mat4rand[6];; double expected4rand_re = std::real(expected4rand); double expected4rand_im = std::imag(expected4rand); @@ -675,16 +677,16 @@ TEST(LoopHafnianEigenDouble, EvenRandom) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum1 = distribution(generator); - mat2[i*2+j] = randnum1; + mat2[i * 2 + j] = randnum1; } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); - mat4[i*4+j] = randnum1; + mat4[i * 4 + j] = randnum1; } } @@ -692,11 +694,11 @@ TEST(LoopHafnianEigenDouble, EvenRandom) { double haf4 = hafnian::loop_hafnian_eigen(mat4); double expected2 = mat2[1] + mat2[0] * mat2[3]; - double expected4 = mat4[1]*mat4[11] + mat4[2]*mat4[7] + mat4[3]*mat4[6] + - mat4[0]*mat4[5]*mat4[11] + mat4[1]*mat4[10]*mat4[15] + - mat4[2]*mat4[5]*mat4[15] + mat4[0]*mat4[10]*mat4[7] + - mat4[0]*mat4[15]*mat4[6] + mat4[3]*mat4[5]*mat4[10] + - mat4[0]*mat4[5]*mat4[10]*mat4[15]; + double expected4 = mat4[1] * mat4[11] + mat4[2] * mat4[7] + mat4[3] * mat4[6] + + mat4[0] * mat4[5] * mat4[11] + mat4[1] * mat4[10] * mat4[15] + + mat4[2] * mat4[5] * mat4[15] + mat4[0] * mat4[10] * mat4[7] + + mat4[0] * mat4[15] * mat4[6] + mat4[3] * mat4[5] * mat4[10] + + mat4[0] * mat4[5] * mat4[10] * mat4[15]; EXPECT_NEAR(expected2, haf2, tol); EXPECT_NEAR(expected4, haf4, tol); @@ -742,18 +744,18 @@ TEST(LoopHafnianEigenComplex, EvenRandom) { std::normal_distribution distribution(1.0, 0.0); for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat2[i*2+j] = std::complex(randnum1, randnum2); + mat2[i * 2 + j] = std::complex(randnum1, randnum2); } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat4[i*4+j] = std::complex(randnum1, randnum2); + mat4[i * 4 + j] = std::complex(randnum1, randnum2); } } @@ -761,11 +763,11 @@ TEST(LoopHafnianEigenComplex, EvenRandom) { std::complex haf4 = hafnian::loop_hafnian_eigen(mat4); std::complex expected2 = mat2[1] + mat2[0] * mat2[3]; - std::complex expected4 = mat4[1]*mat4[11] + mat4[2]*mat4[7] + mat4[3]*mat4[6] + - mat4[0]*mat4[5]*mat4[11] + mat4[1]*mat4[10]*mat4[15] + - mat4[2]*mat4[5]*mat4[15] + mat4[0]*mat4[10]*mat4[7] + - mat4[0]*mat4[15]*mat4[6] + mat4[3]*mat4[5]*mat4[10] + - mat4[0]*mat4[5]*mat4[10]*mat4[15]; + std::complex expected4 = mat4[1] * mat4[11] + mat4[2] * mat4[7] + mat4[3] * mat4[6] + + mat4[0] * mat4[5] * mat4[11] + mat4[1] * mat4[10] * mat4[15] + + mat4[2] * mat4[5] * mat4[15] + mat4[0] * mat4[10] * mat4[7] + + mat4[0] * mat4[15] * mat4[6] + mat4[3] * mat4[5] * mat4[10] + + mat4[0] * mat4[5] * mat4[10] * mat4[15]; EXPECT_NEAR(std::real(expected2), std::real(haf2), tol); EXPECT_NEAR(std::imag(expected2), std::imag(haf2), tol); @@ -824,10 +826,10 @@ TEST(LoopHafnianRepeatedDouble, EvenOnes) { std::vector rpt6(6, 1); for (int i = 0; i < 4; i++) - mu4[i] = mat4[i*4+i]; + mu4[i] = mat4[i * 4 + i]; for (int i = 0; i < 6; i++) - mu6[i] = mat6[i*6+i]; + mu6[i] = mat6[i * 6 + i]; double haf4 = hafnian::loop_hafnian_rpt_quad(mat4, mu4, rpt4); double haf6 = hafnian::loop_hafnian_rpt_quad(mat6, mu6, rpt6); @@ -857,34 +859,34 @@ TEST(LoopHafnianRepeatedDouble, EvenRandom) { for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum1 = distribution(generator); - mat2[i*2+j] = randnum1; + mat2[i * 2 + j] = randnum1; } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); - mat4[i*4+j] = randnum1; + mat4[i * 4 + j] = randnum1; } } for (int i = 0; i < 2; i++) - mu2[i] = mat2[i*2+i]; + mu2[i] = mat2[i * 2 + i]; for (int i = 0; i < 4; i++) - mu4[i] = mat4[i*4+i]; + mu4[i] = mat4[i * 4 + i]; double haf2 = hafnian::loop_hafnian_rpt_quad(mat2, mu2, rpt2); double haf4 = hafnian::loop_hafnian_rpt_quad(mat4, mu4, rpt4); double expected2 = mat2[1] + mat2[0] * mat2[3]; - double expected4 = mat4[1]*mat4[11] + mat4[2]*mat4[7] + mat4[3]*mat4[6] + - mat4[0]*mat4[5]*mat4[11] + mat4[1]*mat4[10]*mat4[15] + - mat4[2]*mat4[5]*mat4[15] + mat4[0]*mat4[10]*mat4[7] + - mat4[0]*mat4[15]*mat4[6] + mat4[3]*mat4[5]*mat4[10] + - mat4[0]*mat4[5]*mat4[10]*mat4[15]; + double expected4 = mat4[1] * mat4[11] + mat4[2] * mat4[7] + mat4[3] * mat4[6] + + mat4[0] * mat4[5] * mat4[11] + mat4[1] * mat4[10] * mat4[15] + + mat4[2] * mat4[5] * mat4[15] + mat4[0] * mat4[10] * mat4[7] + + mat4[0] * mat4[15] * mat4[6] + mat4[3] * mat4[5] * mat4[10] + + mat4[0] * mat4[5] * mat4[10] * mat4[15]; EXPECT_NEAR(expected2, haf2, tol); EXPECT_NEAR(expected4, haf4, tol); @@ -903,10 +905,10 @@ TEST(LoopHafnianRepeatedDouble, Odd) { std::vector rpt5(5, 1); for (int i = 0; i < 3; i++) - mu3[i] = mat3[i*3+i]; + mu3[i] = mat3[i * 3 + i]; for (int i = 0; i < 5; i++) - mu5[i] = mat5[i*5+i]; + mu5[i] = mat5[i * 5 + i]; double haf3 = hafnian::loop_hafnian_rpt_quad(mat3, mu3, rpt3); double haf5 = hafnian::loop_hafnian_rpt_quad(mat5, mu5, rpt5); @@ -954,13 +956,13 @@ TEST(LoopHafnianRepeatedComplex, EvenOnes) { std::vector rpt6(6, 1); for (int i = 0; i < 4; i++) - mu4[i] = mat4[i*4+i]; + mu4[i] = mat4[i * 4 + i]; std::complex haf4 = hafnian::loop_hafnian_rpt_quad(mat4, mu4, rpt4); EXPECT_NEAR(10, std::real(haf4), tol); EXPECT_NEAR(0, std::imag(haf4), tol); for (int i = 0; i < 6; i++) - mu6[i] = mat6[i*6+i]; + mu6[i] = mat6[i * 6 + i]; std::complex haf6 = hafnian::loop_hafnian_rpt_quad(mat6, mu6, rpt6); @@ -988,36 +990,36 @@ TEST(LoopHafnianRepeatedComplex, EvenRandom) { for (int i = 0; i < 2; i++) { - for(int j = 0; j < 2; j++) { + for (int j = 0; j < 2; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat2[i*2+j] = std::complex(randnum1, randnum2); + mat2[i * 2 + j] = std::complex(randnum1, randnum2); } } for (int i = 0; i < 4; i++) { - for(int j = 0; j < 4; j++) { + for (int j = 0; j < 4; j++) { double randnum1 = distribution(generator); double randnum2 = distribution(generator); - mat4[i*4+j] = std::complex(randnum1, randnum2); + mat4[i * 4 + j] = std::complex(randnum1, randnum2); } } for (int i = 0; i < 2; i++) - mu2[i] = mat2[i*2+i]; + mu2[i] = mat2[i * 2 + i]; for (int i = 0; i < 4; i++) - mu4[i] = mat4[i*4+i]; + mu4[i] = mat4[i * 4 + i]; std::complex haf2 = hafnian::loop_hafnian_rpt_quad(mat2, mu2, rpt2); std::complex haf4 = hafnian::loop_hafnian_rpt_quad(mat4, mu4, rpt4); std::complex expected2 = mat2[1] + mat2[0] * mat2[3]; - std::complex expected4 = mat4[1]*mat4[11] + mat4[2]*mat4[7] + mat4[3]*mat4[6] + - mat4[0]*mat4[5]*mat4[11] + mat4[1]*mat4[10]*mat4[15] + - mat4[2]*mat4[5]*mat4[15] + mat4[0]*mat4[10]*mat4[7] + - mat4[0]*mat4[15]*mat4[6] + mat4[3]*mat4[5]*mat4[10] + - mat4[0]*mat4[5]*mat4[10]*mat4[15]; + std::complex expected4 = mat4[1] * mat4[11] + mat4[2] * mat4[7] + mat4[3] * mat4[6] + + mat4[0] * mat4[5] * mat4[11] + mat4[1] * mat4[10] * mat4[15] + + mat4[2] * mat4[5] * mat4[15] + mat4[0] * mat4[10] * mat4[7] + + mat4[0] * mat4[15] * mat4[6] + mat4[3] * mat4[5] * mat4[10] + + mat4[0] * mat4[5] * mat4[10] * mat4[15]; EXPECT_NEAR(std::real(expected2), std::real(haf2), tol); EXPECT_NEAR(std::imag(expected2), std::imag(haf2), tol); @@ -1039,10 +1041,10 @@ TEST(LoopHafnianRepeatedComplex, Odd) { std::vector rpt5(5, 1); for (int i = 0; i < 3; i++) - mu3[i] = mat3[i*3+i]; + mu3[i] = mat3[i * 3 + i]; for (int i = 0; i < 5; i++) - mu5[i] = mat5[i*5+i]; + mu5[i] = mat5[i * 5 + i]; std::complex haf3 = hafnian::loop_hafnian_rpt_quad(mat3, mu3, rpt3); std::complex haf5 = hafnian::loop_hafnian_rpt_quad(mat5, mu5, rpt5); @@ -1068,16 +1070,16 @@ TEST(TorontonianDouble, TMSV) { double r = asinh(std::sqrt(mean_n)); int n = 4; - for(int i = 0; i < n; i++) - mat4[i*n+n-i-1] = tanh(r)*1.0; + for (int i = 0; i < n; i++) + mat4[i * n + n - i - 1] = tanh(r) * 1.0; n = 8; - for(int i = 0; i < n; i++) - mat8[i*n+n-i-1] = tanh(r)*1.0; + for (int i = 0; i < n; i++) + mat8[i * n + n - i - 1] = tanh(r) * 1.0; n = 16; - for(int i = 0; i < n; i++) - mat16[i*n+n-i-1] = tanh(r)*1.0; + for (int i = 0; i < n; i++) + mat16[i * n + n - i - 1] = tanh(r) * 1.0; double tor4 = hafnian::torontonian_quad(mat4); double tor8 = hafnian::torontonian_quad(mat8); @@ -1093,7 +1095,7 @@ TEST(TorontonianDouble, TMSV) { TEST(TorontonianDouble, Vacuum) { int n_modes = 5; - std::vector mat(2*n_modes*2*n_modes, 0.0); + std::vector mat(2 * n_modes * 2 * n_modes, 0.0); double tor_val = hafnian::torontonian_quad(mat); @@ -1103,11 +1105,11 @@ TEST(TorontonianDouble, Vacuum) { TEST(TorontonianDouble, Analytical) { int n = 1; double nbar = 0.25; - std::vector mat1(2*n*2*n, 0.0); - for(int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - mat1[i*2*n+j] = nbar/(static_cast(n)*(1.0+nbar)); - mat1[(i+n)*2*n+(j+n)] = nbar/(static_cast(n)*(1.0+nbar)); + std::vector mat1(2 * n * 2 * n, 0.0); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + mat1[i * 2 * n + j] = nbar / (static_cast(n) * (1.0 + nbar)); + mat1[(i + n) * 2 * n + (j + n)] = nbar / (static_cast(n) * (1.0 + nbar)); } } @@ -1117,11 +1119,11 @@ TEST(TorontonianDouble, Analytical) { n = 2; nbar = 0.25; - std::vector mat2(2*n*2*n, 0.0); - for(int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - mat2[i*2*n+j] = nbar/(static_cast(n)*(1.0+nbar)); - mat2[(i+n)*2*n+(j+n)] = nbar/(static_cast(n)*(1.0+nbar)); + std::vector mat2(2 * n * 2 * n, 0.0); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + mat2[i * 2 * n + j] = nbar / (static_cast(n) * (1.0 + nbar)); + mat2[(i + n) * 2 * n + (j + n)] = nbar / (static_cast(n) * (1.0 + nbar)); } } @@ -1135,3 +1137,177 @@ TEST(TorontonianDouble, Analytical) { } + +namespace batchhafnian { +TEST(BatchHafnian, CompleteGraph) { + std::vector> mat4{std::complex(-0.28264629150778969, 0.39867701584672210), std::complex(-0.06086128222348247, -0.12220227033305252), std::complex(-0.22959477315790058, 0.00000000000000008), std::complex(-0.00660678867199307, -0.09884501458235322), std::complex(-0.06086128222348247, -0.12220227033305252), std::complex(0.38245649793510783, -0.41413300040003126), std::complex(-0.00660678867199307, 0.09884501458235322), std::complex(-0.13684045954832844, 0.00000000000000006), std::complex(-0.22959477315790058, -0.00000000000000008), std::complex(-0.00660678867199307, 0.09884501458235322), std::complex(-0.28264629150778969, -0.39867701584672210), std::complex(-0.06086128222348247, 0.12220227033305252), std::complex(-0.00660678867199307, -0.09884501458235322), std::complex(-0.13684045954832844, -0.00000000000000006), std::complex(-0.06086128222348247, +0.12220227033305252), std::complex(0.38245649793510783, 0.41413300040003126)}; + std::vector> d4{std::complex(0.66917130190858, -1.52776303400764), std::complex(-2.95847055822102, -1.29582519437023), std::complex(0.66917130190858, 1.52776303400764), std::complex(-2.95847055822102, 1.29582519437023)}; + std::vector> out(256, 0.0); + std::vector expected_re{1.00000000e+00, -1.64614736e+00, 1.94351456e+00, -1.44618627e+00, + 4.35642368e-01, -1.32047906e+00, 2.23766490e+00, -1.86917564e+00, + -6.76966967e-01, 5.73670333e-01, -7.33188149e-02, -1.21997190e-01, + 2.32161778e-01, -5.57198229e-01, 1.18563164e+00, -1.79235874e+00, + -1.64614736e+00, 3.23047167e+00, -4.65694018e+00, 4.44401287e+00, + -4.63159381e-02, 1.31073870e+00, -3.22177207e+00, 3.63237405e+00, + 1.23991893e+00, -1.44213928e+00, 8.01092161e-01, 2.28567603e-01, + -6.99231782e-01, 1.51022665e+00, -2.91603997e+00, 4.30125549e+00, + 1.94351456e+00, -4.65694018e+00, 8.15053238e+00, -9.76981613e+00, + -7.95376620e-01, -3.06685257e-01, 2.99900529e+00, -5.19576276e+00, + -1.40243198e+00, 2.03208134e+00, -1.62929470e+00, -2.50514870e-01, + 1.12880996e+00, -2.69285454e+00, 5.33966585e+00, -8.00210813e+00, + -1.44618627e+00, 4.44401287e+00, -9.76981613e+00, 1.52284285e+01, + 1.11641813e+00, -8.67834158e-01, -1.21356826e+00, 4.85970544e+00, + 4.82389499e-01, -9.16927653e-01, 9.19681681e-01, 5.20655922e-01, + -8.34195372e-01, 2.92229672e+00, -6.78809904e+00, 1.10996783e+01, + 4.35642368e-01, -4.63159381e-02, -7.95376620e-01, 1.11641813e+00, + 1.56877658e+00, -2.31606077e+00, 2.25637651e+00, -1.05357461e+00, + 4.79157988e-01, -1.44851326e+00, 2.00611116e+00, -5.18263640e-01, + -6.15836536e-01, 5.36430564e-01, -2.88165936e-01, 4.59048908e-01, + -1.32047906e+00, 1.31073870e+00, -3.06685257e-01, -8.67834158e-01, + -2.31606077e+00, 4.15988484e+00, -5.13015068e+00, 3.41102851e+00, + 1.55349628e-01, 1.12665346e+00, -2.55195698e+00, 1.37496136e+00, + 5.89607772e-01, -4.94801663e-01, -7.51210237e-02, 5.26517347e-01, + 2.23766490e+00, -3.22177207e+00, 2.99900529e+00, -1.21356826e+00, + 2.25637651e+00, -5.13015068e+00, 8.00781690e+00, -7.39736288e+00, + -1.05510937e+00, 1.44993971e-01, 1.83226945e+00, -2.23251713e+00, + 2.06934935e-01, -6.05111407e-01, 1.55663126e+00, -2.69619939e+00, + -1.86917564e+00, 3.63237405e+00, -5.19576276e+00, 4.85970544e+00, + -1.05357461e+00, 3.41102851e+00, -7.39736288e+00, 1.03896013e+01, + 6.44552723e-01, -6.96168471e-01, -2.62839607e-01, 1.95309353e+00, + -1.42746493e+00, 2.66108892e+00, -4.01938103e+00, 4.99610368e+00, + -6.76966967e-01, 1.23991893e+00, -1.40243198e+00, 4.82389499e-01, + 4.79157988e-01, 1.55349628e-01, -1.05510937e+00, 6.44552723e-01, + 2.08027043e+00, -2.71815189e+00, 2.08749438e+00, -3.59011119e-01, + 2.39920807e-01, -1.06932525e+00, 1.14339407e+00, 9.25081052e-01, + 5.73670333e-01, -1.44213928e+00, 2.03208134e+00, -9.16927653e-01, + -1.44851326e+00, 1.12665346e+00, 1.44993971e-01, -6.96168471e-01, + -2.71815189e+00, 4.45234470e+00, -4.61154260e+00, 1.74127108e+00, + 6.25361244e-01, 3.75915531e-01, -1.21876790e+00, -5.13059479e-01, + -7.33188149e-02, 8.01092161e-01, -1.62929470e+00, 9.19681681e-01, + 2.00611116e+00, -2.55195698e+00, 1.83226945e+00, -2.62839607e-01, + 2.08749438e+00, -4.61154260e+00, 6.55944518e+00, -4.66420500e+00, + -1.28536521e+00, 7.81945188e-01, 5.63969365e-01, -7.03167365e-01, + -1.21997190e-01, 2.28567603e-01, -2.50514870e-01, 5.20655922e-01, + -5.18263640e-01, 1.37496136e+00, -2.23251713e+00, 1.95309353e+00, + -3.59011119e-01, 1.74127108e+00, -4.66420500e+00, 7.36229153e+00, + 1.42777111e-02, -3.18771206e-01, -5.24341968e-02, 1.68507163e+00, + 2.32161778e-01, -6.99231782e-01, 1.12880996e+00, -8.34195372e-01, + -6.15836536e-01, 5.89607772e-01, 2.06934935e-01, -1.42746493e+00, + 2.39920807e-01, 6.25361244e-01, -1.28536521e+00, 1.42777111e-02, + 2.99123653e+00, -3.75195474e+00, 2.82213119e+00, -7.81358057e-01, + -5.57198229e-01, 1.51022665e+00, -2.69285454e+00, 2.92229672e+00, + 5.36430564e-01, -4.94801663e-01, -6.05111407e-01, 2.66108892e+00, + -1.06932525e+00, 3.75915531e-01, 7.81945188e-01, -3.18771206e-01, + -3.75195474e+00, 6.10393167e+00, -6.48641622e+00, 3.30150978e+00, + 1.18563164e+00, -2.91603997e+00, 5.33966585e+00, -6.78809904e+00, + -2.88165936e-01, -7.51210237e-02, 1.55663126e+00, -4.01938103e+00, + 1.14339407e+00, -1.21876790e+00, 5.63969365e-01, -5.24341968e-02, + 2.82213119e+00, -6.48641622e+00, 9.98950443e+00, -9.09570553e+00, + -1.79235874e+00, 4.30125549e+00, -8.00210813e+00, 1.10996783e+01, + 4.59048908e-01, 5.26517347e-01, -2.69619939e+00, 4.99610368e+00, + 9.25081052e-01, -5.13059479e-01, -7.03167365e-01, 1.68507163e+00, + -7.81358057e-01, 3.30150978e+00, -9.09570553e+00, 1.56159162e+01}; + std::vector expected_im{0.00000000e+00, -6.19540212e-01, 1.62557597e+00, -2.04268077e+00, + -1.07209959e+00, 1.37273368e+00, -1.04855753e+00, 2.44875659e-01, + -5.35426993e-01, 1.06382831e+00, -1.13167436e+00, -5.50895821e-02, + 2.33832577e-01, -3.78336056e-01, 7.66353380e-01, -1.42492084e+00, + 6.19540212e-01, -8.32667268e-17, -1.64140851e+00, 3.13391670e+00, + 1.93588686e+00, -3.06589811e+00, 3.30672773e+00, -1.86213343e+00, + 3.61695050e-01, -1.18989013e+00, 1.65240909e+00, -7.67200444e-02, + -5.09572526e-02, 1.60560058e-01, -6.31216079e-01, 1.58488036e+00, + -1.62557597e+00, 1.64140851e+00, -2.00406005e-15, -2.89715424e+00, + -2.45819767e+00, 4.73270217e+00, -6.50126689e+00, 5.38692527e+00, + 1.31935650e-01, 6.45842049e-01, -1.46639188e+00, 1.74105305e-01, + -7.03114482e-01, 9.82371295e-01, -7.30475300e-01, -4.69877164e-01, + 2.04268077e+00, -3.13391670e+00, 2.89715424e+00, -9.52005210e-15, + 1.83179431e+00, -4.41516458e+00, 7.85181522e+00, -9.37606710e+00, + -2.10974632e-01, 5.27110853e-02, 3.02027970e-01, 3.02537735e-01, + 1.73508640e+00, -2.84340169e+00, 3.30111115e+00, -2.05599107e+00, + 1.07209959e+00, -1.93588686e+00, 2.45819767e+00, -1.83179431e+00, + -5.55111512e-17, -9.23929364e-01, 2.07252046e+00, -1.72348979e+00, + -1.45132762e+00, 1.63837310e+00, -9.25631147e-01, -8.65199492e-02, + -1.80257925e-02, -4.95003549e-03, 7.10339932e-01, -2.21351707e+00, + -1.37273368e+00, 3.06589811e+00, -4.73270217e+00, 4.41516458e+00, + 9.23929364e-01, -6.97358837e-16, -1.96212896e+00, 2.69923158e+00, + 2.26051161e+00, -3.21251263e+00, 2.71847065e+00, -3.91405380e-01, + -4.80113320e-01, 7.08913327e-01, -1.77090980e+00, 4.06599993e+00, + 1.04855753e+00, -3.30672773e+00, 6.50126689e+00, -7.85181522e+00, + -2.07252046e+00, 1.96212896e+00, -2.15452656e-15, -2.46721598e+00, + -2.22778040e+00, 3.99810349e+00, -4.56372028e+00, 1.81520753e+00, + 8.17667372e-01, -1.51861788e+00, 3.08750334e+00, -5.85075242e+00, + -2.44875659e-01, 1.86213343e+00, -5.38692527e+00, 9.37606710e+00, + 1.72348979e+00, -2.69923158e+00, 2.46721598e+00, -8.65973959e-15, + 7.76349098e-01, -2.07828888e+00, 3.67741704e+00, -3.39674623e+00, + 1.67268754e-03, 1.05174705e+00, -3.26705990e+00, 6.24576761e+00, + 5.35426993e-01, -3.61695050e-01, -1.31935650e-01, 2.10974632e-01, + 1.45132762e+00, -2.26051161e+00, 2.22778040e+00, -7.76349098e-01, + -4.49459509e-16, -1.15371273e+00, 2.14384564e+00, -7.75127965e-01, + -1.69420825e+00, 1.75565162e+00, -7.87144208e-01, -2.91297712e-01, + -1.06382831e+00, 1.18989013e+00, -6.45842049e-01, -5.27110853e-02, + -1.63837310e+00, 3.21251263e+00, -3.99810349e+00, 2.07828888e+00, + 1.15371273e+00, -2.03309591e-15, -1.93589426e+00, 1.65877008e+00, + 2.16797656e+00, -2.87072370e+00, 1.92413422e+00, 6.63359067e-01, + 1.13167436e+00, -1.65240909e+00, 1.46639188e+00, -3.02027970e-01, + 9.25631147e-01, -2.71847065e+00, 4.56372028e+00, -3.67741704e+00, + -2.14384564e+00, 1.93589426e+00, -3.84414722e-15, -1.84525302e+00, + -1.51949589e+00, 2.78597545e+00, -2.83838661e+00, 1.47024726e-02, + 5.50895821e-02, 7.67200444e-02, -1.74105305e-01, -3.02537735e-01, + 8.65199492e-02, 3.91405380e-01, -1.81520753e+00, 3.39674623e+00, + 7.75127965e-01, -1.65877008e+00, 1.84525302e+00, -8.71004657e-15, + 5.32280868e-02, -7.29284151e-01, 2.11417457e+00, -2.58202829e+00, + -2.33832577e-01, 5.09572526e-02, 7.03114482e-01, -1.73508640e+00, + 1.80257925e-02, 4.80113320e-01, -8.17667372e-01, -1.67268754e-03, + 1.69420825e+00, -2.16797656e+00, 1.51949589e+00, -5.32280868e-02, + -1.26663794e-15, -1.59424776e+00, 2.80268987e+00, -9.62872951e-01, + 3.78336056e-01, -1.60560058e-01, -9.82371295e-01, 2.84340169e+00, + 4.95003549e-03, -7.08913327e-01, 1.51861788e+00, -1.05174705e+00, + -1.75565162e+00, 2.87072370e+00, -2.78597545e+00, 7.29284151e-01, + 1.59424776e+00, -3.49720253e-15, -2.61480295e+00, 2.53217806e+00, + -7.66353380e-01, 6.31216079e-01, 7.30475300e-01, -3.30111115e+00, + -7.10339932e-01, 1.77090980e+00, -3.08750334e+00, 3.26705990e+00, + 7.87144208e-01, -1.92413422e+00, 2.83838661e+00, -2.11417457e+00, + -2.80268987e+00, 2.61480295e+00, -6.24500451e-15, -3.10967275e+00, + 1.42492084e+00, -1.58488036e+00, 4.69877164e-01, 2.05599107e+00, + 2.21351707e+00, -4.06599993e+00, 5.85075242e+00, -6.24576761e+00, + 2.91297712e-01, -6.63359067e-01, -1.47024726e-02, 2.58202829e+00, + 9.62872951e-01, -2.53217806e+00, 3.10967275e+00, -1.42559575e-14}; + int res = 4; + int renorm = 0; + + out = hafnian::hermite_multidimensional_cpp(mat4, d4, res, renorm); + + for (int i = 0; i < 256; i++) { + EXPECT_NEAR(expected_re[i], std::real(out[i]), tol2); + EXPECT_NEAR(expected_im[i], std::imag(out[i]), tol2); + } + +} + + + +TEST(BatchHafnian, UnitRenormalization) { + std::vector> B = {std::complex(0, 0), std::complex(-0.70710678, 0), std::complex(-0.70710678, 0), std::complex(0, 0)}; + std::vector> d(4, std::complex(0.0, 0.0)); + + int res = 10; + + std::vector expected_re(res*res, 0); + std::vector expected_im(res*res, 0); + + std::vector> out(res*res, 0.0); + + int renorm = 1; + + for (int i = 0; i < res; i++) + expected_re[i*res+i] = pow(0.5, static_cast(i)/2.0); + + out = hafnian::hermite_multidimensional_cpp(B, d, res, renorm); + + for (int i = 0; i < res*res; i++) { + EXPECT_NEAR(expected_re[i], std::real(out[i]), tol2); + EXPECT_NEAR(expected_im[i], std::imag(out[i]), tol2); + } + +} + +}