Speed up LinearQubitOperator matvec

[blazingphoenix7]

Second of the operator hot-path findings from #1407; the first is #1415.

LinearQubitOperator._matvec applies a QubitOperator to a state vector, and it
is the inner loop behind generate_linear_qubit_operator and the eigensolvers
built on top of it. For each term it took the input vector, split it into 2**k
pieces to reach the first non-identity qubit, split each of those in half, applied
the 2x2 Pauli to the pairs, and concatenated everything back, once per Pauli in
the term. That is a lot of Python-level list building and array copying for what
is, per term, a single reindexing of the amplitudes.

A Pauli string acts on the state vector as a signed permutation, so it can be
written directly on the flat index. Number the amplitudes 0 .. 2**n - 1 and put
qubit q at bit n - 1 - q (the ordering the split version already used, and the
one qubit_operator_sparse uses). Then for one term:

• the X and Y qubits flip their bit, so amplitude c moves to c ^ x_mask;
• the Y and Z qubits give a factor of -1 whenever the matching bit of c is set,
  that is (-1) raised to the parity of popcount(c & z_mask);
• each Y gives an extra factor of 1j.

So the term becomes out[idx ^ x_mask] += coeff * 1j**(y_count % 4) * sign * x,
where sign is a vector of +-1 from the bit parity and idx ^ x_mask is a
permutation of the indices, which is what makes the fancy-indexed += safe. The
masks are cheap integer folds over the term and the per-amplitude work is a few
vectorized numpy operations with no intermediate lists.

This is an implementation change only. Public behavior is unchanged: complex
output, inputs of shape (dim,) or (dim, 1), the identity term, operators
declared on more qubits than they act on, and the empty operator all behave as
before. Only the serial _matvec is touched; ParallelLinearQubitOperator is
left alone.

The parity of popcount uses a small xor-fold helper on uint64 rather than
numpy.bitwise_count, because the project still supports numpy 1.26 (the
max_compat CI job pins numpy==1.26.4) and bitwise_count needs numpy 2.0.

Correctness

linear_qubit_operator_test.py passes unchanged. On top of that I compared the
new matvec against qubit_operator_sparse(op) @ x over random operators for
n = 2..12 (mixed X/Y/Z, term counts including the empty and identity-only
operators), real and complex x, shapes (dim,) and (dim, 1), and n_qubits
oversized by 2, plus a scipy.sparse.linalg.eigsh run checked against dense
eigenvalues. Worst absolute differences:

random real x, shape (dim,)          5.4e-15
random complex x, shape (dim,)       9.1e-15
random complex x, shape (dim, 1)     9.1e-15
oversized n_qubits by 2              2.4e-15
identity-only and empty operator     0
eigsh vs dense (3 smallest)          8.9e-15

black and mypy are clean on the file.

Motivation

This is what #499 asked for. That issue called for a better _matvec, measured
that one was worth 80-100x at system sizes of 16 qubits and up, and noted that
ParallelLinearQubitOperator could then be deleted as obsolete once the serial
path was fast. The sketch there routed each Pauli through Cirq's
apply_unitary; no port ever landed and the issue was closed for inactivity in
2025. This PR gets the same win with plain numpy on the flat index, so it adds
no dependency, and the 120x at 16 qubits below bears out the issue's estimate.

The parallel path is also where the deadlock possibility in #1405 lives, and
with the serial matvec this fast there is rarely a reason to reach for it.
Removing it, as #499 suggested, would be a separate change; this PR fixes
neither issue.

Benchmarks

Random QubitOperators with mixed X/Y/Z terms applied to a complex state vector:

n_qubits   terms    before     after    speedup
   16       100     20.5 s     0.17 s      120x
   18       100     87.0 s     1.5 s        60x
   20        50      157 s     3.0 s        52x

Measured on my laptop (Intel Core Ultra 5 225, Windows, Python 3.12, numpy 2.5).
I ran the base/patch comparison twice in alternation with fixed seeds, taking
medians over three runs at n = 16 and single runs at n = 18 and 20 since the old
code takes minutes there; the per-round speedups agreed within a few percent, so
the table pools both rounds.

[gemini-code-assist]

Code Review

This pull request refactors the _matvec method in LinearQubitOperator to use bitwise operations and masks, introducing a new _bit_parity helper function to apply operators as signed permutations instead of recursively splitting vectors. The review feedback points out that using 1j ** (y_count % 4) introduces numerical inaccuracies due to floating-point approximations and suggests using a lookup table [1, 1j, -1, -1j][y_count % 4] to keep the real and imaginary parts exact.

[gemini-code-assist]

Using 1j ** (y_count % 4) introduces numerical inaccuracies (spurious imaginary parts of order $10^{-16}$) due to floating-point approximations in Python's complex exponentiation. For example, 1j ** 2 evaluates to (-1+1.2246467991473532e-16j) instead of exactly -1. In quantum simulation, keeping real/imaginary parts exact is important to avoid accumulating numerical noise.

Using a pre-defined lookup table like [1, 1j, -1, -1j][y_count % 4] is 100% exact and also faster.

Suggested change

	amplitudes = coefficient * (1j ** (y_count % 4)) * vec
	amplitudes = coefficient * [1, 1j, -1, -1j][y_count % 4] * vec

[blazingphoenix7]

Thanks, took the table. One note for the thread: 1j ** (y_count % 4) is actually already exact here since the exponent is a Python int, so CPython evaluates it by repeated multiplication rather than the polar form the example uses. The table makes that explicit though, and keeps it exact if the exponent ever comes through as a numpy scalar.