Merge pull request #98 from Budapest-Quantum-Computing-Group/GregoryM…

…orse-patch-4

Fix Hafnian accuracy issue

piquassoboost/sampling/Boson_Sampling_Utilities.py

@@ -1,5 +1,5 @@
 #
-# Copyright 2021-2022 Budapest Quantum Computing Group
+# Copyright 2021 Budapest Quantum Computing Group
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -361,6 +361,59 @@ def calculate(self):
 
 
 
+class PowerTraceHafnianLongDouble(PowerTraceHafnian_wrapper):
+    """
+        This class is designed to calculate the hafnian of a symetrix matrix using the power trace method.
+    """
+
+
+    def __init__(self, matrix):
+
+        # call the constructor of the wrapper class
+        super(PowerTraceHafnianLongDouble, self).__init__(lib=2, matrix=matrix)
+        pass
+
+
+    def calculate(self):
+        """
+            ?????????????????.
+            :return: The hafnian of the matrix.
+        """
+
+
+
+        # call the permanent calculator of the parent class
+        return super(PowerTraceHafnianLongDouble, self).calculate()
+
+
+
+
+class PowerTraceHafnianInf(PowerTraceHafnian_wrapper):
+    """
+        This class is designed to calculate the hafnian of a symetrix matrix using the power trace method.
+    """
+
+
+    def __init__(self, matrix):
+
+        # call the constructor of the wrapper class
+        super(PowerTraceHafnianInf, self).__init__(lib=3, matrix=matrix)
+        pass
+
+
+    def calculate(self):
+        """
+            ?????????????????.
+            :return: The hafnian of the matrix.
+        """
+
+
+
+        # call the permanent calculator of the parent class
+        return super(PowerTraceHafnianInf, self).calculate()
+
+
+
 
 
 class PowerTraceHafnianRecursive(PowerTraceHafnianRecursive_wrapper):
@@ -391,21 +444,21 @@ def calculate(self):
 
 class PowerTraceLoopHafnian(PowerTraceLoopHafnian_wrapper):
     """
-        This class is designed to calculate the loop hafnian of a symetrix matrix using the power trace method.
+        This class is designed to calculate the loop hafnian of a symetrix matrix using the power trace method calculated with double precision.
     """
 
 
     def __init__(self, matrix):
 
         # call the constructor of the wrapper class
-        super(PowerTraceLoopHafnian, self).__init__(matrix=matrix)
+        super(PowerTraceLoopHafnian, self).__init__(matrix=matrix, lib=1)
         pass
 
 
     def calculate(self):
         """
             ?????????????????.
-            :return: The hafnian of the matrix.
+            :return: The loop hafnian of the matrix.
         """
 
 
@@ -415,9 +468,61 @@ def calculate(self):
 
 
 
+class PowerTraceLoopHafnianLongDouble(PowerTraceLoopHafnian_wrapper):
+    """
+        This class is designed to calculate the loop hafnian of a symetrix matrix using the power trace method calculated with extended precision.
+    """
+
+
+    def __init__(self, matrix):
+
+        # call the constructor of the wrapper class
+        super(PowerTraceLoopHafnianLongDouble, self).__init__(lib=2, matrix=matrix )
+        pass
+
+
+    def calculate(self):
+        """
+            ?????????????????.
+            :return: The loop hafnian of the matrix.
+        """
+
+
+
+        # call the permanent calculator of the parent class
+        return super(PowerTraceLoopHafnianLongDouble, self).calculate()
+
+
+
+class PowerTraceLoopHafnianInf(PowerTraceLoopHafnian_wrapper):
+    """
+        This class is designed to calculate the loop hafnian of a symetrix matrix using the power trace method calculated with infinite precision (using the MPFR library).
+    """
+
+
+    def __init__(self, matrix):
+
+        # call the constructor of the wrapper class
+        super(PowerTraceLoopHafnianInf, self).__init__(lib=3, matrix=matrix )
+        pass
+
+
+    def calculate(self):
+        """
+            ?????????????????.
+            :return: The loop hafnian of the matrix.
+        """
+
+
+
+        # call the permanent calculator of the parent class
+        return super(PowerTraceLoopHafnianInf, self).calculate()        
+
+
+
 class PowerTraceLoopHafnianRecursive(PowerTraceLoopHafnianRecursive_wrapper):
     """
-        This class is designed to calculate the hafnian of a symetrix matrix using the power trace method.
+        This class is designed to calculate the loop hafnian of a symetrix matrix using the power trace method.
     """
 
 
@@ -431,7 +536,7 @@ def __init__(self, matrix, occupancy):
     def calculate(self):
         """
             ?????????????????.
-            :return: The hafnian of the matrix.
+            :return: The loop hafnian of the matrix.
         """
 
 

piquassoboost/sampling/source/PowerTraceHafnian.cpp

@@ -438,9 +438,9 @@ PowerTraceHafnian<small_scalar_type, scalar_type>::ScaleMatrix() {
         // determine the scale factor
         scale_factor = 0.0;
         for (size_t idx=0; idx<mtx_orig.size(); idx++) {
-            scale_factor = scale_factor + std::sqrt( mtx_orig[idx].real()*mtx_orig[idx].real() + mtx_orig[idx].imag()*mtx_orig[idx].imag() );
+            scale_factor = scale_factor + sqrt( mtx_orig[idx].real()*mtx_orig[idx].real() + mtx_orig[idx].imag()*mtx_orig[idx].imag() );
         }
-        scale_factor = scale_factor/mtx_orig.size()/std::sqrt(2);
+        scale_factor = scale_factor/mtx_orig.size()/sqrt(2.0);
         //scale_factor = scale_factor*mtx_orig.rows;
 
         mtx = mtx_orig.copy();

piquassoboost/sampling/source/PowerTraceHafnianUtilities.cpp

@@ -22,7 +22,7 @@
 #include "loop_correction_AVX.h"
 
 
-
+long double sqrt(long double arg) { return sqrtl(arg); }
 
 namespace pic {
 

piquassoboost/sampling/source/PowerTraceHafnianUtilities.hpp

@@ -38,6 +38,8 @@ static double time_szamlalo = 0.0;
 static double time_nevezo = 0.0;
 */
 
+long double sqrt(long double arg);
+
 namespace pic {
 
 
@@ -63,7 +65,7 @@ get_reflection_vector(matrix_type &input, small_scalar_type &norm_v_sqr) {
   sigma = sqrt(norm_v_sqr);
 
 
-  small_scalar_type abs_val = std::sqrt( reflect_vector[0].real()*reflect_vector[0].real() + reflect_vector[0].imag()*reflect_vector[0].imag() );
+  small_scalar_type abs_val = std::abs(reflect_vector[0]);
   norm_v_sqr = 2*(norm_v_sqr + abs_val*sigma);
   if (abs_val != 0.0){
       //small_scalar_type angle = std::arg(reflect_vector[0]); // sigma *= (reflect_vector[0] / std::abs(reflect_vector[0]));
@@ -79,7 +81,7 @@ get_reflection_vector(matrix_type &input, small_scalar_type &norm_v_sqr) {
       return reflect_vector;
 
   // normalize the reflection matrix
-  small_scalar_type norm_v = std::sqrt(norm_v_sqr);
+  small_scalar_type norm_v = sqrt(norm_v_sqr);
   for (size_t idx=0; idx<reflect_vector.size(); idx++) {
       reflect_vector[idx] = reflect_vector[idx]/norm_v;
   }

piquassoboost/sampling/source/PowerTraceLoopHafnian.cpp

@@ -438,9 +438,9 @@ PowerTraceLoopHafnian<small_scalar_type, scalar_type>::ScaleMatrix() {
         // determine the scale factor
         this->scale_factor = 0.0;
         for (size_t idx=0; idx<this->mtx_orig.size(); idx++) {
-            this->scale_factor = this->scale_factor + std::sqrt( this->mtx_orig[idx].real()*this->mtx_orig[idx].real() + this->mtx_orig[idx].imag()*this->mtx_orig[idx].imag() );
+            this->scale_factor = this->scale_factor + sqrt( this->mtx_orig[idx].real()*this->mtx_orig[idx].real() + this->mtx_orig[idx].imag()*this->mtx_orig[idx].imag() );
         }
-        this->scale_factor = this->scale_factor/this->mtx_orig.size()/std::sqrt(2);
+        this->scale_factor = this->scale_factor/this->mtx_orig.size()/sqrt(2.0);
 
         this->mtx = this->mtx_orig.copy();