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);