From 1da7e239e457b9220b4330045a02f28cfefce776 Mon Sep 17 00:00:00 2001
From: vera-rubin <322284620+vera-rubin@users.noreply.github.com>
Date: Thu, 15 Oct 2026 06:22:21 +0000
Subject: [PATCH] Fix zeroth powers of singular diagonal matrices

---
 quest/src/api/calculations.cpp |  9 +++++----
 tests/unit/calculations.cpp    | 29 +++++++++++++++++++++++++++++
 2 files changed, 34 insertions(+), 4 deletions(-)

diff --git a/quest/src/api/calculations.cpp b/quest/src/api/calculations.cpp
index 47e5d8a63..8e1717128 100644
--- a/quest/src/api/calculations.cpp
+++ b/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):
diff --git a/tests/unit/calculations.cpp b/tests/unit/calculations.cpp
index 4b4db284e..f7c389587 100644
--- a/tests/unit/calculations.cpp
+++ b/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
 }