Add SparseObservable.compose

This adds matrix multiplication between `SparseObservable`s (or one
`SparseObservable` and anything that can be coerced to an observable),
using the same signature as similar methods in `quantum_info`.

The current implementation attempts to be relatively memory efficient,
though there are specialisations that may be implemented if we wanted to
be faster:

- if all the contained terms are Paulis, so each pairwise multiplication
  produces either zero or one alphabet letters, we currently have a
  redundant copy-in and copy-out of the Cartesian-product generator that
  effectively doubles the overhead.  We could check whether each term
  contains only Paulis once, and then avoid that overhead in those
  cases.

- alternatively, the Cartesian-product iterator could be taught to use
  the buffers of the output `SparseObservable` directly as its output
  space, rather than requiring copy-out of its buffers, which would
  remove the overhead from all calls, not just the single-Pauli ones,
  but requires a bit more bookkeeping, since the references would keep
  need to be updated.

- the case of `qargs` being not `None` and `front=True` currently
  involves an extra copy step to simplify the logic.  If this is an
  important case, we could add a right-matmul specialisation of
  `compose_map`.

This commit currently does not contain tests, because we currently don't
have suitable methods to test the validity of the output without
coupling to the particular choice of factorisation of the matrix
multiplication.  We need either a way to convert to a Pauli-only
representation (with the exponential explosion that entails) or to a
matrix; either of these can be uniquely canonicalised.

crates/accelerate/src/sparse_observable/lookup.rs

@@ -0,0 +1,289 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2025
+//
+// This code is licensed under the Apache License, Version 2.0. You may
+// obtain a copy of this license in the LICENSE.txt file in the root directory
+// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+//
+// Any modifications or derivative works of this code must retain this
+// copyright notice, and modified files need to carry a notice indicating
+// that they have been altered from the originals.
+
+use num_complex::Complex64;
+// The lookup tables are pretty illegible if we have all the syntactic noise of `BitTerm`.
+use super::BitTerm::{self, *};
+
+/// Short-hand alias for [Complex64::new] that retains its ability to be used in `const` contexts.
+const fn c64(re: f64, im: f64) -> Complex64 {
+    Complex64::new(re, im)
+}
+
+/// The allowable items of `BitTerm`.  This is used by the lookup expansion; we need a const-safe
+/// way of iterating through all of the variants.
+static BIT_TERMS: [BitTerm; 9] = [X, Plus, Minus, Y, Right, Left, Z, Zero, One];
+
+/// Create a full binary lookup table for all valid combinations of `BitTerm`.
+macro_rules! expand_lookup_binary {
+    (static $table:ident: [$ty:ty], $generate:ident, $dummy:expr) => {
+        static $table: [$ty; 256] = {
+            let mut out = [const { $dummy }; 256];
+            let mut i = 0;
+            while i < 9 {
+                let left = BIT_TERMS[i];
+                let mut j = 0;
+                while j < 9 {
+                    let right = BIT_TERMS[j];
+                    out[((left as usize) << 4) | (right as usize)] = $generate(left, right);
+                    j += 1;
+                }
+                i += 1;
+            }
+            out
+        };
+    };
+}
+
+expand_lookup_binary!(
+    static MATMUL: [Option<&'static [(Complex64, BitTerm)]>],
+    matmul_generate,
+    None
+);
+/// Calculate the matrix multiplication of two [BitTerm]s.
+///
+/// This is done by a static lookup table.
+///
+/// In the output, `None` means that the result was zero.  Beyond that, the slice should be
+/// interpreted as a sum of the coefficient multiplied by the [BitTerm] pairs.  An empty slice means
+/// the identity.
+#[inline(always)]
+pub const fn matmul(left: BitTerm, right: BitTerm) -> Option<&'static [(Complex64, BitTerm)]> {
+    MATMUL[((left as usize) << 4) | (right as usize)]
+}
+
+const fn matmul_generate(left: BitTerm, right: BitTerm) -> Option<&'static [(Complex64, BitTerm)]> {
+    // This typically gets compiled to a non-inlinable (because the code size is too big) two
+    // separate sets of operate-and-jump instructions, but we use it in the `const` context to build
+    // a static lookup table which can easily be inlined and is a couple of bitwise operations
+    // followed by an offset load of the result.
+    //
+    // These rules were generated by a brute-force search over increasing-n-tuples of (coeff, term)
+    // pairs, where `coeff` was restricted to magnitudes `(0.25, 0.5, 1.0)` at complex phases
+    // `(1, i, -1, -i)`.  Under those conditions, the rules are minimal in number of sum terms, but
+    // don't have any further logic to "curate" the choice of output factorisation.
+    match (left, right) {
+        (X, X) => Some(const { &[] }),
+        (X, Plus) => Some(const { &[(c64(1., 0.), Plus)] }),
+        (X, Minus) => Some(const { &[(c64(-1., 0.), Minus)] }),
+        (X, Y) => Some(const { &[(c64(0., 1.), Z)] }),
+        (X, Right) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Z)] }),
+        (X, Left) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Z)] }),
+        (X, Z) => Some(const { &[(c64(0., -1.), Y)] }),
+        (X, Zero) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Y)] }),
+        (X, One) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Y)] }),
+        (Plus, X) => Some(const { &[(c64(1., 0.), Plus)] }),
+        (Plus, Plus) => Some(const { &[(c64(1., 0.), Plus)] }),
+        (Plus, Minus) => None,
+        (Plus, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), Z)] }),
+        (Plus, Right) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), X),
+                    (c64(0.5, 0.), Right),
+                    (c64(0., 0.25), Z),
+                ]
+            },
+        ),
+        (Plus, Left) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), X),
+                    (c64(0.5, 0.), Left),
+                    (c64(0., -0.25), Z),
+                ]
+            },
+        ),
+        (Plus, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), Y)] }),
+        (Plus, Zero) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), X),
+                    (c64(0.5, 0.), Zero),
+                    (c64(0., -0.25), Y),
+                ]
+            },
+        ),
+        (Plus, One) => {
+            Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., 0.25), Y)] })
+        }
+        (Minus, X) => Some(const { &[(c64(-1., 0.), Minus)] }),
+        (Minus, Plus) => None,
+        (Minus, Minus) => Some(const { &[(c64(1., 0.), Minus)] }),
+        (Minus, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), Z)] }),
+        (Minus, Right) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Y),
+                    (c64(0.5, 0.), Minus),
+                    (c64(0., -0.25), Z),
+                ]
+            },
+        ),
+        (Minus, Left) => Some(
+            const {
+                &[
+                    (c64(-0.25, 0.), X),
+                    (c64(0.5, 0.), Left),
+                    (c64(0., 0.25), Z),
+                ]
+            },
+        ),
+        (Minus, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), Y)] }),
+        (Minus, Zero) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Z),
+                    (c64(0.5, 0.), Minus),
+                    (c64(0., 0.25), Y),
+                ]
+            },
+        ),
+        (Minus, One) => Some(
+            const {
+                &[
+                    (c64(-0.25, 0.), X),
+                    (c64(0.5, 0.), One),
+                    (c64(0., -0.25), Y),
+                ]
+            },
+        ),
+        (Y, X) => Some(const { &[(c64(0., -1.), Z)] }),
+        (Y, Plus) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), Z)] }),
+        (Y, Minus) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), Z)] }),
+        (Y, Y) => Some(const { &[] }),
+        (Y, Right) => Some(const { &[(c64(1., 0.), Right)] }),
+        (Y, Left) => Some(const { &[(c64(-1., 0.), Left)] }),
+        (Y, Z) => Some(const { &[(c64(0., 1.), X)] }),
+        (Y, Zero) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), X)] }),
+        (Y, One) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), X)] }),
+        (Right, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Z)] }),
+        (Right, Plus) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), X),
+                    (c64(0.5, 0.), Right),
+                    (c64(0., -0.25), Z),
+                ]
+            },
+        ),
+        (Right, Minus) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Y),
+                    (c64(0.5, 0.), Minus),
+                    (c64(0., 0.25), Z),
+                ]
+            },
+        ),
+        (Right, Y) => Some(const { &[(c64(1., 0.), Right)] }),
+        (Right, Right) => Some(const { &[(c64(1., 0.), Right)] }),
+        (Right, Left) => None,
+        (Right, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), X)] }),
+        (Right, Zero) => {
+            Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), Zero), (c64(0., 0.25), X)] })
+        }
+        (Right, One) => {
+            Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., -0.25), X)] })
+        }
+        (Left, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Z)] }),
+        (Left, Plus) => {
+            Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), Left), (c64(0., 0.25), Z)] })
+        }
+        (Left, Minus) => Some(
+            const {
+                &[
+                    (c64(-0.25, 0.), X),
+                    (c64(0.5, 0.), Left),
+                    (c64(0., -0.25), Z),
+                ]
+            },
+        ),
+        (Left, Y) => Some(const { &[(c64(-1., 0.), Left)] }),
+        (Left, Right) => None,
+        (Left, Left) => Some(const { &[(c64(1., 0.), Left)] }),
+        (Left, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), X)] }),
+        (Left, Zero) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Z),
+                    (c64(0.5, 0.), Left),
+                    (c64(0., -0.25), X),
+                ]
+            },
+        ),
+        (Left, One) => {
+            Some(const { &[(c64(-0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., 0.25), X)] })
+        }
+        (Z, X) => Some(const { &[(c64(0., 1.), Y)] }),
+        (Z, Plus) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), Y)] }),
+        (Z, Minus) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), Y)] }),
+        (Z, Y) => Some(const { &[(c64(0., -1.), X)] }),
+        (Z, Right) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), X)] }),
+        (Z, Left) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), X)] }),
+        (Z, Z) => Some(const { &[] }),
+        (Z, Zero) => Some(const { &[(c64(1., 0.), Zero)] }),
+        (Z, One) => Some(const { &[(c64(-1., 0.), One)] }),
+        (Zero, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Y)] }),
+        (Zero, Plus) => {
+            Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), Zero), (c64(0., 0.25), Y)] })
+        }
+        (Zero, Minus) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Z),
+                    (c64(0.5, 0.), Minus),
+                    (c64(0., -0.25), Y),
+                ]
+            },
+        ),
+        (Zero, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), X)] }),
+        (Zero, Right) => Some(
+            const {
+                &[
+                    (c64(0.25, 0.), Y),
+                    (c64(0.5, 0.), Zero),
+                    (c64(0., -0.25), X),
+                ]
+            },
+        ),
+        (Zero, Left) => {
+            Some(const { &[(c64(0.25, 0.), Z), (c64(0.5, 0.), Left), (c64(0., 0.25), X)] })
+        }
+        (Zero, Z) => Some(const { &[(c64(1., 0.), Zero)] }),
+        (Zero, Zero) => Some(const { &[(c64(1., 0.), Zero)] }),
+        (Zero, One) => None,
+        (One, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Y)] }),
+        (One, Plus) => {
+            Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., -0.25), Y)] })
+        }
+        (One, Minus) => {
+            Some(const { &[(c64(-0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., 0.25), Y)] })
+        }
+        (One, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), X)] }),
+        (One, Right) => {
+            Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., 0.25), X)] })
+        }
+        (One, Left) => Some(
+            const {
+                &[
+                    (c64(-0.25, 0.), Y),
+                    (c64(0.5, 0.), One),
+                    (c64(0., -0.25), X),
+                ]
+            },
+        ),
+        (One, Z) => Some(const { &[(c64(-1., 0.), One)] }),
+        (One, Zero) => None,
+        (One, One) => Some(const { &[(c64(1., 0.), One)] }),
+    }
+}