Extend solve_qubit to work beyond 63 qubits #240
Description
What should we add?
At the moment, solve_qubit() and project_operator_to_subspace() throw the following error:
if bitstring_matrix.shape[1] > 63:
raise ValueError("Bitstrings (rows) in bitstring_matrix must have length < 64.")
which is explained in the documentation:
.. note::
This function relies on ``jax`` to efficiently perform some calculations. ``jax``
converts the bit arrays to ``int64_t``, which means the bit arrays in
``bitstring_matrix`` may not have length greater than ``63``.
Whilst understanding where this limitation comes from, it would particularly high value to adapt the method to work for largest systems. For example, this would allow running SKQD of spin systems (e.g., https://quantum.cloud.ibm.com/learning/en/courses/quantum-diagonalization-algorithms/skqd) at utility scales of 100+ qubits.
I realise the scope of this change is potentially very large, if for example there is no int128_t type to use and a whole different backend is needed.