Fix zeroth powers of singular diagonal matrices

quest/src/api/calculations.cpp

@@ -101,10 +101,11 @@ qcomp calcExpecNonHermitianFullStateDiagMatrPower(Qureg qureg, FullStateDiagMatr
     validate_matrixFields(matrix, __func__);
     validate_matrixAndQuregAreCompatible(matrix, qureg, true, __func__);
 
-    // this function never uses the qreal-pow overload (because we make
-    // no assumption matrix is Hermitian i.e. real), and instead always
-    // uses the relatively numerically inaccurate qcomp overload
-    const bool useRealPow = false;
+    // real pow(0,0) returns 1, whereas complex pow(0,0) can return NaN.
+    // For exponent zero, every element becomes 1 even when complex,
+    // so we can safely discard imaginary components and use real pow.
+    // Otherwise, we make no assumption matrix is Hermitian and use complex pow.
+    const bool useRealPow = exponent == qcomp(0,0);
 
     return (qureg.isDensityMatrix)?
         localiser_densmatr_calcExpecFullStateDiagMatr(qureg, matrix, exponent, useRealPow):

tests/unit/calculations.cpp

@@ -700,6 +700,35 @@ TEST_CASE( "calcExpecNonHermitianFullStateDiagMatrPower", TEST_CATEGORY LABEL_MI
         }
     }
 
+    SECTION( "zeroth power of singular matrix" ) {
+
+        qcomp elem = GENERATE( qcomp(1,0), qcomp(-1,0), qcomp(0,1) );
+        qmatrix refMatr = getDiagonalMatrix(getConstantVector(getPow2(getNumCachedQubits()), elem));
+        refMatr[0][0] = 0;
+
+        auto apiFunc = [](Qureg qureg, FullStateDiagMatr matr) {
+            return calcExpecNonHermitianFullStateDiagMatrPower(qureg, matr, 0);
+        };
+
+        SECTION( LABEL_STATEVEC ) {
+
+            auto refFunc = [](qvector state, qmatrix) {
+                return getInnerProduct(state, state);
+            };
+
+            TEST_ON_CACHED_QUREG_AND_MATRIX( getCachedStatevecs(), getCachedFullStateDiagMatrs(), apiFunc, getRefStatevec(), refMatr, refFunc);
+        }
+
+        SECTION( LABEL_DENSMATR ) {
+
+            auto refFunc = [](qmatrix state, qmatrix) {
+                return getTrace(state);
+            };
+
+            TEST_ON_CACHED_QUREG_AND_MATRIX( getCachedDensmatrs(), getCachedFullStateDiagMatrs(), apiFunc, getRefDensmatr(), refMatr, refFunc);
+        }
+    }
+
     /// @todo input validation
 }