Added C++ permanent source codes and removed Fortran ones

hafnian/__init__.py

@@ -81,11 +81,13 @@
     haf_rpt_real,
     haf_rpt_complex,
     reduction,
+    perm_complex,
+    perm_real,
 )
 
 from ._torontonian import tor
 
-from ._permanent import perm, perm_real, perm_complex, permanent_repeated
+from ._permanent import perm, permanent_repeated
 
 __all__ = [
     "hafnian",

hafnian/_hafnian.py

@@ -16,7 +16,7 @@
 """
 import numpy as np
 
-from .lib.libhaf import haf_complex, haf_real, haf_int, haf_rpt_real, haf_rpt_complex
+from .lib.libhaf import haf_complex, haf_real, haf_int, haf_rpt_real, haf_rpt_complex, perm_complex, perm_real
 
 
 def input_validation(A, tol=1e-12):

hafnian/_permanent.py

@@ -16,14 +16,12 @@
 """
 import numpy as np
 
-from .lib.libperm import perm as libperm
-from ._hafnian import hafnian_repeated
+from .lib.libhaf import perm_complex, perm_real
 
-perm_real = libperm.re
-perm_complex = libperm.comp
+from ._hafnian import hafnian_repeated
 
 
-def perm(A):
+def perm(A, quad=True, fsum=False):
     """Returns the permanent of a matrix via the
     `Ryser formula <https://en.wikipedia.org/wiki/Computing_the_permanent#Ryser_formula>`_.
 
@@ -32,6 +30,9 @@ def perm(A):
 
     Args:
         A (array): a square array.
+        quad (bool): quad (bool): If ``True``, the input matrix is cast to a ``long double``
+            matrix internally for a quadruple precision hafnian computation.
+        fsum (bool): If ``True``, ``fsum`` method is used for summation.
 
     Returns:
         np.float64 or np.complex128: the permanent of matrix A.
@@ -63,10 +64,10 @@ def perm(A):
 
     if A.dtype == np.complex:
         if np.any(np.iscomplex(A)):
-            return perm_complex(A)
-        return perm_real(np.float64(A.real))
+            return perm_complex(A, quad=quad)
+        return perm_real(np.float64(A.real), quad=quad, fsum=fsum)
 
-    return perm_real(A)
+    return perm_real(A, quad=quad, fsum=fsum)
 
 
 def permanent_repeated(A, rpt):