hilbert_space: Reduce hilbert_space::max_n_bits to 63

This way `hilbert_space::dim()` can return a valid value of type
`sv_index_type` even when all 63 bits are used up.

CHANGELOG.md

@@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
 
 ## [0.8.0] - Unreleased
 
+- Reduce the maximum allowed number of bits in the binary representation of
+ a basis state index (``hilbert_space::max_n_bits``) to 63. This way
+ ``hilbert_space::dim()`` can return a valid value of type ``sv_index_type``
+ even when all 63 bits are used up.
 - Fixed a negative index bug in ``n_fermion_sector_view``.
  Credits to Dr. Cezary Śliwa for providing the patch.
 - Whenever possible, use compiler intrinsics to speed up complex bit

doc/loperator/hilbert_space.rst

@@ -55,11 +55,15 @@ In the code, the basis vectors are represented by 64-bit unsigned integers
 (:type:`sv_index_type`).
 The binary form of :math:`|n\rangle_\mathcal{H}` is
 then a concatenation of binary forms of :math:`|n_i\rangle_{\mathcal{H}_i}`.
-For example, the following picture shows memory representation of basis
+For example, the following picture shows memory representation of the basis
 state :math:`|90\rangle_\mathcal{H} = |0\rangle_{\mathcal{H}_1} \otimes
 |5\rangle_{\mathcal{H}_2} \otimes |1\rangle_{\mathcal{H}_3} \otimes
 |1\rangle_{\mathcal{H}_4}`. Cells (bits) of the same color belong to the same
-elementary space and the higher blank cells 7-63 are unused -- set to zero.
+elementary space and the higher blank cells 7-62 are unused -- set to zero. The
+most significant bit 63 is reserved and can never be occupied by an elementary
+space. This restriction is necessary to make sure that the size of a Hilbert
+space with all the allowed 63 bits used up is still representable by a value
+of the type :type:`sv_index_type`.
 
 .. image:: ../images/basis_state.svg
  :width: 800
@@ -153,7 +157,7 @@ does not matter -- they will be reordered automatically.
  Construct from a list of elementary spaces. The elementary spaces need not
  be given in any particular order.
  Throws :struct:`hilbert_space_too_big` if all elementary spaces together
- would overflow the 64-bit integer type of the basis state index.
+ would exceed the 63-bit limit of the basis state index.
 
  .. function:: template<typename ScalarType, \
  typename ESConstructor = default_es_constructor> \
@@ -166,8 +170,7 @@ does not matter -- they will be reordered automatically.
  Construction of the elementary spaces is performed by the functor
  :expr:`es_constr`.
  Throws :struct:`hilbert_space_too_big` if all collected elementary spaces
- together would overflow the 64-bit integer type of the basis state index.
-
+ together would exceed the 63-bit limit of the basis state index.
 
  .. rubric:: Copy/move-constructors and assignments
 
@@ -191,7 +194,7 @@ does not matter -- they will be reordered automatically.
  :struct:`elementary_space_exists` if an elementary space equivalent to
  :expr:`es` is already part of the product.
  Throws :struct:`hilbert_space_too_big` if adding :expr:`es` into the product
- would overflow the 64-bit integer type of the basis state index.
+ would exceed the 63-bit limit of the basis state index.
 
  .. function:: bool has(elementary_space<IndexTypes...> const& es) const
 
@@ -233,8 +236,8 @@ does not matter -- they will be reordered automatically.
  The total number of used bits in the binary representation of a basis state
  index.
 
- .. function:: std::size_t dim() const
- friend std::size_t get_dim(hilbert_space const& hs)
+ .. function:: sv_index_type dim() const
+ friend sv_index_type get_dim(hilbert_space const& hs)
 
  The dimension of this Hilbert space computed as a product of dimensions
  of the elementary spaces.
@@ -266,7 +269,7 @@ does not matter -- they will be reordered automatically.
 
  .. struct:: hilbert_space_too_big : public std::runtime_error
 
- The total basis state index size exceeds 64 bits.
+ The total basis state index size exceeds 63 bits.
 
 .. function:: template<typename ScalarType, \
  typename... IndexTypes, \

include/libcommute/loperator/hilbert_space.hpp

@@ -46,7 +46,8 @@ template <typename... IndexTypes> class hilbert_space {
 
  // Size of a Hilbert space must be representable by an integer with
  // at most this number of bits.
- static constexpr int max_n_bits = std::numeric_limits<sv_index_type>::digits;
+ static constexpr int max_n_bits =
+ std::numeric_limits<sv_index_type>::digits - 1;
 
  // Compare two elementary_space_t objects wrapped in std::unique_ptr<T>
  struct less {
@@ -94,7 +95,7 @@ template <typename... IndexTypes> class hilbert_space {
  explicit hilbert_space_too_big(int n_bits)
  : std::runtime_error("Hilbert space size is not representable "
  "by a " +
- std::to_string(max_n_bits) +
+ std::to_string(max_n_bits + 1) +
  "-bit integer (n_bits = " + std::to_string(n_bits) +
  ")"),
  n_bits(n_bits) {}
@@ -227,8 +228,8 @@ template <typename... IndexTypes> class hilbert_space {
  int total_n_bits() const { return bit_range_end_ + 1; }
 
  // Dimension of this Hilbert space
- std::size_t dim() const { return std::size_t(1) << total_n_bits(); }
- friend std::size_t get_dim(hilbert_space const& hs) { return hs.dim(); }
+ sv_index_type dim() const { return sv_index_type(1) << total_n_bits(); }
+ friend sv_index_type get_dim(hilbert_space const& hs) { return hs.dim(); }
 
  // Apply functor `f` to all basis state indices
  template <typename Functor>

test/hilbert_space.cpp

@@ -348,17 +348,18 @@ TEST_CASE("Hilbert space", "[hilbert_space]") {
 
  SECTION("Very big Hilbert space") {
  hs_type hs1;
- for(int i = 0; i < 32; ++i)
+ for(int i = 0; i < 31; ++i)
  hs1.add(make_space_spin(3.0 / 2, "s", i));
 
  using ex_type = hs_type::hilbert_space_too_big;
- CHECK_THROWS_AS(hs1.add(make_space_spin(3.0 / 2, "s", 32)), ex_type);
+ CHECK_THROWS_AS(hs1.add(make_space_spin(3.0 / 2, "s", 31)), ex_type);
 
  auto expr = expression<double, std::string, int>(1);
- for(int i = 0; i < 32; ++i)
+ for(int i = 0; i < 31; ++i)
  expr *= S_p<4>("s", i);
+ expr *= S_p<2>("s", 31);
  hs_type hs2(expr);
- CHECK(hs2.total_n_bits() == 64);
+ CHECK(hs2.total_n_bits() == 63);
 
  expr *= S_p<4>("s", 32);
  CHECK_THROWS_AS(hs_type(expr), ex_type);