From eeadf261ad03e8d1996a2aac00b6ea4960eb73e0 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 13:51:15 +0100
Subject: [PATCH 1/9] Clip probabilities

---
 qiskit/primitives/sampler.py                  |  2 +-
 qiskit/quantum_info/states/densitymatrix.py   |  4 ++++
 qiskit/quantum_info/states/quantum_state.py   | 12 +++++++++---
 qiskit/quantum_info/states/stabilizerstate.py | 15 +++++++++++++--
 qiskit/quantum_info/states/statevector.py     |  6 +++++-
 5 files changed, 32 insertions(+), 7 deletions(-)

diff --git a/qiskit/primitives/sampler.py b/qiskit/primitives/sampler.py
index ede3747d0111..c606a85d1498 100644
--- a/qiskit/primitives/sampler.py
+++ b/qiskit/primitives/sampler.py
@@ -119,7 +119,7 @@ def _call(
             qargs_list.append(self._qargs_list[i])
         probabilities = [
             Statevector(bound_circuit_to_instruction(circ)).probabilities_dict(
-                qargs=qargs, decimals=16
+                qargs=qargs, decimals=16, clip=True
             )
             for circ, qargs in zip(bound_circuits, qargs_list)
         ]
diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index 435c02e535cf..71e2a4c7281a 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -420,6 +420,8 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -482,6 +484,8 @@ def probabilities(self, qargs=None, decimals=None):
         )
         if decimals is not None:
             probs = probs.round(decimals=decimals)
+        if clip:
+            probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):
diff --git a/qiskit/quantum_info/states/quantum_state.py b/qiskit/quantum_info/states/quantum_state.py
index 640f7fec3c99..0f8719019eed 100644
--- a/qiskit/quantum_info/states/quantum_state.py
+++ b/qiskit/quantum_info/states/quantum_state.py
@@ -196,7 +196,7 @@ def expectation_value(self, oper, qargs=None):
         pass
 
     @abstractmethod
-    def probabilities(self, qargs=None, decimals=None):
+    def probabilities(self, qargs=None, decimals=None, clip=False):
         """Return the subsystem measurement probability vector.
 
         Measurement probabilities are with respect to measurement in the
@@ -207,13 +207,15 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
         """
         pass
 
-    def probabilities_dict(self, qargs=None, decimals=None):
+    def probabilities_dict(self, qargs=None, decimals=None, clip=False):
         """Return the subsystem measurement probability dictionary.
 
         Measurement probabilities are with respect to measurement in the
@@ -229,12 +231,16 @@ def probabilities_dict(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             dict: The measurement probabilities in dict (ket) form.
         """
         return self._vector_to_dict(
-            self.probabilities(qargs=qargs, decimals=decimals), self.dims(qargs), string_labels=True
+            self.probabilities(qargs=qargs, decimals=decimals, clip=clip),
+            self.dims(qargs),
+            string_labels=True,
         )
 
     def sample_memory(self, shots, qargs=None):
diff --git a/qiskit/quantum_info/states/stabilizerstate.py b/qiskit/quantum_info/states/stabilizerstate.py
index 18aa1d541ae1..be51300f8ada 100644
--- a/qiskit/quantum_info/states/stabilizerstate.py
+++ b/qiskit/quantum_info/states/stabilizerstate.py
@@ -318,6 +318,8 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -331,9 +333,12 @@ def probabilities(self, qargs=None, decimals=None):
             place = int(key, 2)
             probs[place] = value
 
+        if clip:
+            probs = np.clip(probs, a_min=0, a_max=1)
+
         return probs
 
-    def probabilities_dict(self, qargs=None, decimals=None):
+    def probabilities_dict(self, qargs=None, decimals=None, clip=False):
         """Return the subsystem measurement probability dictionary.
 
         Measurement probabilities are with respect to measurement in the
@@ -349,6 +354,8 @@ def probabilities_dict(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             dict: The measurement probabilities in dict (ket) form.
@@ -366,7 +373,11 @@ def probabilities_dict(self, qargs=None, decimals=None):
 
         if decimals is not None:
             for key, value in probs.items():
-                probs[key] = round(value, decimals)
+                prob = round(value, decimals)
+                if clip:
+                    prob = np.clip(prob, a_min=0, a_max=1)
+
+                probs[key] = prob
 
         return probs
 
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index e4f7736c5e98..313bfec246d4 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -504,7 +504,7 @@ def expectation_value(self, oper, qargs=None):
         conj = self.conjugate()
         return np.dot(conj.data, val.data)
 
-    def probabilities(self, qargs=None, decimals=None):
+    def probabilities(self, qargs=None, decimals=None, clip=False):
         """Return the subsystem measurement probability vector.
 
         Measurement probabilities are with respect to measurement in the
@@ -515,6 +515,8 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
+            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
+                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -577,6 +579,8 @@ def probabilities(self, qargs=None, decimals=None):
         )
         if decimals is not None:
             probs = probs.round(decimals=decimals)
+        if clip:
+            probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):

From 0b74ba4081909bc66e43cac1a36c61b54eac1e09 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 16:32:42 +0100
Subject: [PATCH 2/9] clip implictly

---
 qiskit/quantum_info/states/densitymatrix.py   |  6 ++----
 qiskit/quantum_info/states/quantum_state.py   | 10 +++-------
 qiskit/quantum_info/states/stabilizerstate.py | 15 +++------------
 qiskit/quantum_info/states/statevector.py     |  8 +++-----
 4 files changed, 11 insertions(+), 28 deletions(-)

diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index 71e2a4c7281a..9c0d1e6964c1 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -420,8 +420,6 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -484,8 +482,8 @@ def probabilities(self, qargs=None, decimals=None):
         )
         if decimals is not None:
             probs = probs.round(decimals=decimals)
-        if clip:
-            probs = np.clip(probs, a_min=0, a_max=1)
+
+        probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):
diff --git a/qiskit/quantum_info/states/quantum_state.py b/qiskit/quantum_info/states/quantum_state.py
index 0f8719019eed..609c0808a9a1 100644
--- a/qiskit/quantum_info/states/quantum_state.py
+++ b/qiskit/quantum_info/states/quantum_state.py
@@ -196,7 +196,7 @@ def expectation_value(self, oper, qargs=None):
         pass
 
     @abstractmethod
-    def probabilities(self, qargs=None, decimals=None, clip=False):
+    def probabilities(self, qargs=None, decimals=None):
         """Return the subsystem measurement probability vector.
 
         Measurement probabilities are with respect to measurement in the
@@ -207,15 +207,13 @@ def probabilities(self, qargs=None, decimals=None, clip=False):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
         """
         pass
 
-    def probabilities_dict(self, qargs=None, decimals=None, clip=False):
+    def probabilities_dict(self, qargs=None, decimals=None):
         """Return the subsystem measurement probability dictionary.
 
         Measurement probabilities are with respect to measurement in the
@@ -231,14 +229,12 @@ def probabilities_dict(self, qargs=None, decimals=None, clip=False):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             dict: The measurement probabilities in dict (ket) form.
         """
         return self._vector_to_dict(
-            self.probabilities(qargs=qargs, decimals=decimals, clip=clip),
+            self.probabilities(qargs=qargs, decimals=decimals),
             self.dims(qargs),
             string_labels=True,
         )
diff --git a/qiskit/quantum_info/states/stabilizerstate.py b/qiskit/quantum_info/states/stabilizerstate.py
index be51300f8ada..8d155ae20219 100644
--- a/qiskit/quantum_info/states/stabilizerstate.py
+++ b/qiskit/quantum_info/states/stabilizerstate.py
@@ -318,8 +318,6 @@ def probabilities(self, qargs=None, decimals=None):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -333,12 +331,10 @@ def probabilities(self, qargs=None, decimals=None):
             place = int(key, 2)
             probs[place] = value
 
-        if clip:
-            probs = np.clip(probs, a_min=0, a_max=1)
-
+        probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
-    def probabilities_dict(self, qargs=None, decimals=None, clip=False):
+    def probabilities_dict(self, qargs=None, decimals=None):
         """Return the subsystem measurement probability dictionary.
 
         Measurement probabilities are with respect to measurement in the
@@ -354,8 +350,6 @@ def probabilities_dict(self, qargs=None, decimals=None, clip=False):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             dict: The measurement probabilities in dict (ket) form.
@@ -373,10 +367,7 @@ def probabilities_dict(self, qargs=None, decimals=None, clip=False):
 
         if decimals is not None:
             for key, value in probs.items():
-                prob = round(value, decimals)
-                if clip:
-                    prob = np.clip(prob, a_min=0, a_max=1)
-
+                prob = np.clip(round(value, decimals), a_min=0, a_max=1)
                 probs[key] = prob
 
         return probs
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index 313bfec246d4..bcffe4428b21 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -504,7 +504,7 @@ def expectation_value(self, oper, qargs=None):
         conj = self.conjugate()
         return np.dot(conj.data, val.data)
 
-    def probabilities(self, qargs=None, decimals=None, clip=False):
+    def probabilities(self, qargs=None, decimals=None):
         """Return the subsystem measurement probability vector.
 
         Measurement probabilities are with respect to measurement in the
@@ -515,8 +515,6 @@ def probabilities(self, qargs=None, decimals=None, clip=False):
                 if None return for all subsystems (Default: None).
             decimals (None or int): the number of decimal places to round
                 values. If None no rounding is done (Default: None).
-            clip (bool): if ``True``, clip the probabilities to ``[0, 1]``
-                to account for possible roundoff errors (Default: ``True``).
 
         Returns:
             np.array: The Numpy vector array of probabilities.
@@ -579,8 +577,8 @@ def probabilities(self, qargs=None, decimals=None, clip=False):
         )
         if decimals is not None:
             probs = probs.round(decimals=decimals)
-        if clip:
-            probs = np.clip(probs, a_min=0, a_max=1)
+
+        probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):

From 568db1627afc54b1c3e3b5580276d4d1ba4a5609 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 16:36:59 +0100
Subject: [PATCH 3/9] add reno

---
 .../clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml | 8 ++++++++
 1 file changed, 8 insertions(+)
 create mode 100644 releasenotes/notes/clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml

diff --git a/releasenotes/notes/clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml b/releasenotes/notes/clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml
new file mode 100644
index 000000000000..0ba63168064a
--- /dev/null
+++ b/releasenotes/notes/clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml
@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    Clip probabilities in the :meth:`.QuantumState.probabilities` and 
+    :meth:`.QuantumState.probabilities_dict` methods to the interval ``[0, 1]``.
+    This fixes roundoff errors where probabilities could e.g. be larger than 1, leading 
+    to errors in the shot emulation of the :class:`.Sampler`.
+    Fixed `#9761 <https://github.com/Qiskit/qiskit-terra/issues/9761>`__.

From 9c08d776dd7d23ac5423a1eb292109ce8f928662 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 17:09:23 +0100
Subject: [PATCH 4/9] add tests

---
 qiskit/primitives/sampler.py                          | 2 +-
 qiskit/quantum_info/states/stabilizerstate.py         | 4 +---
 test/python/quantum_info/states/test_densitymatrix.py | 7 +++++++
 test/python/quantum_info/states/test_statevector.py   | 7 +++++++
 4 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/qiskit/primitives/sampler.py b/qiskit/primitives/sampler.py
index c606a85d1498..ede3747d0111 100644
--- a/qiskit/primitives/sampler.py
+++ b/qiskit/primitives/sampler.py
@@ -119,7 +119,7 @@ def _call(
             qargs_list.append(self._qargs_list[i])
         probabilities = [
             Statevector(bound_circuit_to_instruction(circ)).probabilities_dict(
-                qargs=qargs, decimals=16, clip=True
+                qargs=qargs, decimals=16
             )
             for circ, qargs in zip(bound_circuits, qargs_list)
         ]
diff --git a/qiskit/quantum_info/states/stabilizerstate.py b/qiskit/quantum_info/states/stabilizerstate.py
index 8d155ae20219..18aa1d541ae1 100644
--- a/qiskit/quantum_info/states/stabilizerstate.py
+++ b/qiskit/quantum_info/states/stabilizerstate.py
@@ -331,7 +331,6 @@ def probabilities(self, qargs=None, decimals=None):
             place = int(key, 2)
             probs[place] = value
 
-        probs = np.clip(probs, a_min=0, a_max=1)
         return probs
 
     def probabilities_dict(self, qargs=None, decimals=None):
@@ -367,8 +366,7 @@ def probabilities_dict(self, qargs=None, decimals=None):
 
         if decimals is not None:
             for key, value in probs.items():
-                prob = np.clip(round(value, decimals), a_min=0, a_max=1)
-                probs[key] = prob
+                probs[key] = round(value, decimals)
 
         return probs
 
diff --git a/test/python/quantum_info/states/test_densitymatrix.py b/test/python/quantum_info/states/test_densitymatrix.py
index 6a619094ce3e..fcd43979e442 100644
--- a/test/python/quantum_info/states/test_densitymatrix.py
+++ b/test/python/quantum_info/states/test_densitymatrix.py
@@ -1202,6 +1202,13 @@ def test_drawings(self):
             with self.subTest(msg=f"draw('{drawtype}')"):
                 dm.draw(drawtype)
 
+    def test_clip_probabilities(self):
+        """Test probabilities are clipped to [0, 1]."""
+        dm = DensityMatrix([[1.1, 0], [0, 0]])
+
+        self.assertTrue(np.allclose(dm.probabilities(), [1, 0], atol=0))
+        self.assertDictAlmostEqual(dm.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index 050c4b09da07..b54ebc915689 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -1308,6 +1308,13 @@ def test_statevector_len(self):
         self.assertEqual(len(empty_vector), len(empty_sv))
         self.assertEqual(len(dummy_vector), len(sv))
 
+    def test_clip_probabilities(self):
+        """Test probabilities are clipped to [0, 1]."""
+        sv = Statevector([1.1, 0])
+
+        self.assertTrue(np.allclose(sv.probabilities(), [1, 0], atol=0))
+        self.assertDictAlmostEqual(sv.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
+
 
 if __name__ == "__main__":
     unittest.main()

From 086e9d7d38c7fd16ea4648f55a6135558c18dcd4 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 17:12:17 +0100
Subject: [PATCH 5/9] add renormalization

---
 qiskit/quantum_info/states/densitymatrix.py | 3 ++-
 qiskit/quantum_info/states/statevector.py   | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index 9c0d1e6964c1..cd172ee5a3bf 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -483,7 +483,8 @@ def probabilities(self, qargs=None, decimals=None):
         if decimals is not None:
             probs = probs.round(decimals=decimals)
 
-        probs = np.clip(probs, a_min=0, a_max=1)
+        # to account for roundoff errors, we renormalize and clip
+        probs = np.clip(probs / np.linalg.norm(probs), a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index bcffe4428b21..4a020fb488eb 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -578,7 +578,8 @@ def probabilities(self, qargs=None, decimals=None):
         if decimals is not None:
             probs = probs.round(decimals=decimals)
 
-        probs = np.clip(probs, a_min=0, a_max=1)
+        # to account for roundoff errors, we renormalize and clip
+        probs = np.clip(probs / np.linalg.norm(probs), a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):

From b22df7cdbb99ba448840cee539ffb521b29a3efe Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Thu, 9 Mar 2023 18:31:00 +0100
Subject: [PATCH 6/9] use sum, not norm

---
 qiskit/quantum_info/states/densitymatrix.py | 2 +-
 qiskit/quantum_info/states/statevector.py   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index cd172ee5a3bf..fdc3344dad5b 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -484,7 +484,7 @@ def probabilities(self, qargs=None, decimals=None):
             probs = probs.round(decimals=decimals)
 
         # to account for roundoff errors, we renormalize and clip
-        probs = np.clip(probs / np.linalg.norm(probs), a_min=0, a_max=1)
+        probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index 4a020fb488eb..f251eb3d21e4 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -579,7 +579,7 @@ def probabilities(self, qargs=None, decimals=None):
             probs = probs.round(decimals=decimals)
 
         # to account for roundoff errors, we renormalize and clip
-        probs = np.clip(probs / np.linalg.norm(probs), a_min=0, a_max=1)
+        probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
         return probs
 
     def reset(self, qargs=None):

From 6b3c30a268ac0132a1db448092f65d92fc935142 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Fri, 10 Mar 2023 11:28:38 +0100
Subject: [PATCH 7/9] ensure round still works

---
 qiskit/quantum_info/states/densitymatrix.py          |  6 ++++--
 qiskit/quantum_info/states/statevector.py            |  6 ++++--
 .../python/quantum_info/states/test_densitymatrix.py | 12 ++++++++++++
 test/python/quantum_info/states/test_statevector.py  | 11 +++++++++++
 4 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index fdc3344dad5b..02bb245e5ed3 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -480,11 +480,13 @@ def probabilities(self, qargs=None, decimals=None):
         probs = self._subsystem_probabilities(
             np.abs(self.data.diagonal()), self._op_shape.dims_l(), qargs=qargs
         )
-        if decimals is not None:
-            probs = probs.round(decimals=decimals)
 
         # to account for roundoff errors, we renormalize and clip
         probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
+
+        if decimals is not None:
+            probs = probs.round(decimals=decimals)
+
         return probs
 
     def reset(self, qargs=None):
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index f251eb3d21e4..507da14760aa 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -575,11 +575,13 @@ def probabilities(self, qargs=None, decimals=None):
         probs = self._subsystem_probabilities(
             np.abs(self.data) ** 2, self._op_shape.dims_l(), qargs=qargs
         )
-        if decimals is not None:
-            probs = probs.round(decimals=decimals)
 
         # to account for roundoff errors, we renormalize and clip
         probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
+
+        if decimals is not None:
+            probs = probs.round(decimals=decimals)
+
         return probs
 
     def reset(self, qargs=None):
diff --git a/test/python/quantum_info/states/test_densitymatrix.py b/test/python/quantum_info/states/test_densitymatrix.py
index fcd43979e442..7ab89565d12f 100644
--- a/test/python/quantum_info/states/test_densitymatrix.py
+++ b/test/python/quantum_info/states/test_densitymatrix.py
@@ -1209,6 +1209,18 @@ def test_clip_probabilities(self):
         self.assertTrue(np.allclose(dm.probabilities(), [1, 0], atol=0))
         self.assertDictAlmostEqual(dm.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
 
+    def test_round_probabilities(self):
+        """Test probabilities are correctly rounded.
+
+        This is good to test to ensure clipping, renormalizing and rounding work together.
+        """
+        p = np.sqrt(1 / 3)
+        amplitudes = [p, p, p, 0]
+        dm = DensityMatrix(np.outer(amplitudes, amplitudes))
+        expected = [0.33, 0.33, 0.33, 0]
+
+        self.assertTrue(np.allclose(dm.probabilities(decimals=2), expected))
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index b54ebc915689..37877fca6ce9 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -1315,6 +1315,17 @@ def test_clip_probabilities(self):
         self.assertTrue(np.allclose(sv.probabilities(), [1, 0], atol=0))
         self.assertDictAlmostEqual(sv.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
 
+    def test_round_probabilities(self):
+        """Test probabilities are correctly rounded.
+
+        This is good to test to ensure clipping, renormalizing and rounding work together.
+        """
+        p = np.sqrt(1 / 3)
+        sv = Statevector([p, p, p, 0])
+        expected = [0.33, 0.33, 0.33, 0]
+
+        self.assertTrue(np.allclose(sv.probabilities(decimals=2), expected))
+
 
 if __name__ == "__main__":
     unittest.main()

From a03b4c672c544076fd1e9d93a511f735e06bc471 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Fri, 10 Mar 2023 17:14:38 +0100
Subject: [PATCH 8/9] normalizing seems to re-introduce errors

---
 qiskit/quantum_info/states/densitymatrix.py | 4 ++--
 qiskit/quantum_info/states/statevector.py   | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index 02bb245e5ed3..ac7b672b77ef 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -481,8 +481,8 @@ def probabilities(self, qargs=None, decimals=None):
             np.abs(self.data.diagonal()), self._op_shape.dims_l(), qargs=qargs
         )
 
-        # to account for roundoff errors, we renormalize and clip
-        probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
+        # to account for roundoff errors, we clip
+        probs = np.clip(probs, a_min=0, a_max=1)
 
         if decimals is not None:
             probs = probs.round(decimals=decimals)
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index 507da14760aa..5c63a6628379 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -576,8 +576,8 @@ def probabilities(self, qargs=None, decimals=None):
             np.abs(self.data) ** 2, self._op_shape.dims_l(), qargs=qargs
         )
 
-        # to account for roundoff errors, we renormalize and clip
-        probs = np.clip(probs / np.sum(probs), a_min=0, a_max=1)
+        # to account for roundoff errors, we clip
+        probs = np.clip(probs, a_min=0, a_max=1)
 
         if decimals is not None:
             probs = probs.round(decimals=decimals)

From d624e85a3bd99b266899c60743960e22e181145b Mon Sep 17 00:00:00 2001
From: Jake Lishman <jake.lishman@ibm.com>
Date: Thu, 16 Mar 2023 18:26:44 +0000
Subject: [PATCH 9/9] Tighten floating-point tests

In cases of rounding and clipping, it's important that the
floating-point output is bit-for-bit correct, so the fuzzy tests weren't
ideal (some of these were strict, but it was inconsistent).  It's better
to use `assertEqual` rather than `assertTrue` where possible so we get
better errors on failure.
---
 test/python/quantum_info/states/test_densitymatrix.py | 8 +++++---
 test/python/quantum_info/states/test_statevector.py   | 8 ++++----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/test/python/quantum_info/states/test_densitymatrix.py b/test/python/quantum_info/states/test_densitymatrix.py
index 7ab89565d12f..17f2b22f174e 100644
--- a/test/python/quantum_info/states/test_densitymatrix.py
+++ b/test/python/quantum_info/states/test_densitymatrix.py
@@ -1206,8 +1206,9 @@ def test_clip_probabilities(self):
         """Test probabilities are clipped to [0, 1]."""
         dm = DensityMatrix([[1.1, 0], [0, 0]])
 
-        self.assertTrue(np.allclose(dm.probabilities(), [1, 0], atol=0))
-        self.assertDictAlmostEqual(dm.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
+        self.assertEqual(list(dm.probabilities()), [1.0, 0.0])
+        # The "1" key should be exactly zero and therefore omitted.
+        self.assertEqual(dm.probabilities_dict(), {"0": 1.0})
 
     def test_round_probabilities(self):
         """Test probabilities are correctly rounded.
@@ -1219,7 +1220,8 @@ def test_round_probabilities(self):
         dm = DensityMatrix(np.outer(amplitudes, amplitudes))
         expected = [0.33, 0.33, 0.33, 0]
 
-        self.assertTrue(np.allclose(dm.probabilities(decimals=2), expected))
+        # Exact floating-point check because fixing the rounding should ensure this is exact.
+        self.assertEqual(list(dm.probabilities(decimals=2)), expected)
 
 
 if __name__ == "__main__":
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index 37877fca6ce9..d81cd33d75f8 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -1312,8 +1312,9 @@ def test_clip_probabilities(self):
         """Test probabilities are clipped to [0, 1]."""
         sv = Statevector([1.1, 0])
 
-        self.assertTrue(np.allclose(sv.probabilities(), [1, 0], atol=0))
-        self.assertDictAlmostEqual(sv.probabilities_dict(), {"0": 1, "1": 0}, delta=0)
+        self.assertEqual(list(sv.probabilities()), [1.0, 0.0])
+        # The "1" key should be zero and therefore omitted.
+        self.assertEqual(sv.probabilities_dict(), {"0": 1.0})
 
     def test_round_probabilities(self):
         """Test probabilities are correctly rounded.
@@ -1323,8 +1324,7 @@ def test_round_probabilities(self):
         p = np.sqrt(1 / 3)
         sv = Statevector([p, p, p, 0])
         expected = [0.33, 0.33, 0.33, 0]
-
-        self.assertTrue(np.allclose(sv.probabilities(decimals=2), expected))
+        self.assertEqual(list(sv.probabilities(decimals=2)), expected)
 
 
 if __name__ == "__main__":