From 3e583ec7b23ef18e3c52fd42125d127d5f1aa16c Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sat, 13 Apr 2024 21:54:31 +0300
Subject: [PATCH 01/12] Refactor occupations interface, remove redundant code

---
 src/manybody.jl | 74 ++++++++++++-------------------------------------
 1 file changed, 18 insertions(+), 56 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index dfb7f25a..afe6f795 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -19,7 +19,7 @@ Base.union(sv1::SortedVector{T}, svs::SortedVector{T}...) where {T} =
     SortedVector(union(sv1.sortedvector, (occ.sortedvector for occ in svs)...))
 
 # Special methods for fast operator construction
-allocate_buffer(sv::AbstractVector) = similar(first(sv))
+allocate_buffer(sv::AbstractVector) = ismutable(first(sv)) ? similar(first(sv)) : Ref(first(sv))
 function state_index(sv::SortedVector, state)
     ret = searchsortedfirst(sv.sortedvector, state, order = sv.ord)
     ret == length(sv) + 1 && return nothing
@@ -41,7 +41,7 @@ struct ManyBodyBasis{B,O,UT} <: Basis
     onebodybasis::B
     occupations::O
     occupations_hash::UT
-    function ManyBodyBasis{B,O}(onebodybasis::B, occupations::O) where {B,O<:AbstractVector{Vector{Int}}}
+    function ManyBodyBasis{B,O}(onebodybasis::B, occupations::O) where {B,O<:AbstractVector}
         h = hash(hash.(occupations))
         new{B,O,typeof(h)}(length(occupations), onebodybasis, occupations, h)
     end
@@ -94,26 +94,7 @@ end
 
 Creation operator for the i-th mode of the many-body basis `b`.
 """
-function create(::Type{T}, b::ManyBodyBasis, index) where {T}
-    Is = Int[]
-    Js = Int[]
-    Vs = T[]
-    # <{m}_i| at |{m}_j>
-    buffer = allocate_buffer(b.occupations)
-    for (i, occ_i) in enumerate(b.occupations)
-        if occ_i[index] == 0
-            continue
-        end
-        copyto!(buffer, occ_i)
-        buffer[index] -= 1
-        j = state_index(b.occupations, buffer)
-        j === nothing && continue
-        push!(Is, i)
-        push!(Js, j)
-        push!(Vs, sqrt(occ_i[index]))
-    end
-    return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
-end
+create(::Type{T}, b::ManyBodyBasis, index) where {T} = transition(T, b, index, ())
 create(b::ManyBodyBasis, index) = create(ComplexF64, b, index)
 
 """
@@ -121,26 +102,7 @@ create(b::ManyBodyBasis, index) = create(ComplexF64, b, index)
 
 Annihilation operator for the i-th mode of the many-body basis `b`.
 """
-function destroy(::Type{T}, b::ManyBodyBasis, index) where {T}
-    Is = Int[]
-    Js = Int[]
-    Vs = T[]
-    buffer = allocate_buffer(b.occupations)
-    # <{m}_j| a |{m}_i>
-    for (i, occ_i) in enumerate(b.occupations)
-        if occ_i[index] == 0
-            continue
-        end
-        copyto!(buffer, occ_i)
-        buffer[index] -= 1
-        j = state_index(b.occupations, buffer)
-        j === nothing && continue
-        push!(Is, j)
-        push!(Js, i)
-        push!(Vs, sqrt(occ_i[index]))
-    end
-    return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
-end
+destroy(::Type{T}, b::ManyBodyBasis, index) where {T} = transition(T, b, (), index)
 destroy(b::ManyBodyBasis, index) = destroy(ComplexF64, b, index)
 
 """
@@ -178,7 +140,7 @@ function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T}
     buffer = allocate_buffer(b.occupations)
     # <{m}_j| at_to a_from |{m}_i>
     for (i, occ_i) in enumerate(b.occupations)
-        C = state_transition!(buffer, occ_i, from, to)
+        C = state_transition!(buffer, occ_i, to, from)
         C === nothing && continue
         j = state_index(b.occupations, buffer)
         j === nothing && continue
@@ -236,7 +198,7 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::Operator)
         value = op.data[i, j]
         iszero(value) && continue
         for (m, occ) in enumerate(basis.occupations)
-            C = state_transition!(buffer, occ, i, j)
+            C = state_transition!(buffer, occ, j, i)
             C === nothing && continue
             n = state_index(basis.occupations, buffer)
             n === nothing && continue
@@ -255,7 +217,7 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType)
     buffer = allocate_buffer(basis.occupations)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(basis.occupations)
-            C = state_transition!(buffer, occ, row, column)
+            C = state_transition!(buffer, occ, column, row)
             C === nothing && continue
             n = state_index(basis.occupations, buffer)
             n === nothing && continue
@@ -277,7 +239,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
         value = op_data[i, j, k, l]
         iszero(value) && continue
         for (m, occ) in enumerate(basis.occupations)
-            C = state_transition!(buffer, occ, (i, j), (k, l))
+            C = state_transition!(buffer, occ, (k, l), (i, j))
             C === nothing && continue
             n = state_index(basis.occupations, buffer)
             n === nothing && continue
@@ -297,7 +259,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(basis.occupations)
             index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2+row])
-            C = state_transition!(buffer, occ, index[1:2], index[3:4])
+            C = state_transition!(buffer, occ, index[3:4], index[1:2])
             C === nothing && continue
             n = state_index(basis.occupations, buffer)
             n === nothing && continue
@@ -332,8 +294,8 @@ end
 
 onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state = states]
 
-get_value(state::Ket,  m, n) = conj(state.data[m]) * state.data[n]
-get_value(state::Operator,  m, n) = state.data[n, m]
+get_value(state::Ket, m, n) = conj(state.data[m]) * state.data[n]
+get_value(state::Operator, m, n) = state.data[n, m]
 function onebodyexpect_1(op::Operator, state)
     b = basis(state)
     occupations = b.occupations
@@ -344,7 +306,7 @@ function onebodyexpect_1(op::Operator, state)
         value = op.data[i, j]
         iszero(value) && continue
         for (m, occ) in enumerate(occupations)
-            C = state_transition!(buffer, occ, i, j)
+            C = state_transition!(buffer, occ, j, i)
             C === nothing && continue
             n = state_index(occupations, buffer)
             n === nothing && continue
@@ -361,7 +323,7 @@ function onebodyexpect_1(op::SparseOpPureType, state)
     result = complex(0.0)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(occupations)
-            C = state_transition!(buffer, occ, row, column)
+            C = state_transition!(buffer, occ, column, row)
             C === nothing && continue
             n = state_index(occupations, buffer)
             n === nothing && continue
@@ -371,16 +333,16 @@ function onebodyexpect_1(op::SparseOpPureType, state)
     result
 end
 
-Base.@propagate_inbounds function state_transition!(buffer, occ_m, at_indices, a_indices)
-    any(==(0), (occ_m[m] for m in at_indices)) && return nothing
+Base.@propagate_inbounds function state_transition!(buffer::Vector{Int}, occ::Vector{Int}, at_indices, a_indices)
+    any(==(0), (occ[m] for m in a_indices)) && return nothing
     result = 1
-    copyto!(buffer, occ_m)
-    for i in at_indices
+    copyto!(buffer, occ)
+    for i in a_indices
         result *= buffer[i]
         result == 0 && return nothing
         buffer[i] -= 1
     end
-    for i in a_indices
+    for i in at_indices
         buffer[i] += 1
         result *= buffer[i]
     end

From a5fd9683882c461acf70e11afdf4e1ef31e3a60c Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sat, 13 Apr 2024 22:16:09 +0300
Subject: [PATCH 02/12] Add fermion bitstrings

---
 src/manybody.jl | 57 ++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 54 insertions(+), 3 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index afe6f795..1c28f552 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -20,12 +20,13 @@ Base.union(sv1::SortedVector{T}, svs::SortedVector{T}...) where {T} =
 
 # Special methods for fast operator construction
 allocate_buffer(sv::AbstractVector) = ismutable(first(sv)) ? similar(first(sv)) : Ref(first(sv))
-function state_index(sv::SortedVector, state)
+function state_index(sv::SortedVector{T}, state::T) where {T}
     ret = searchsortedfirst(sv.sortedvector, state, order = sv.ord)
     ret == length(sv) + 1 && return nothing
     return sv.sortedvector[ret] == state ? ret : nothing
 end
-state_index(sv::AbstractVector, state) = findfirst(==(state), sv)
+state_index(sv::AbstractVector{T}, state::T) where {T} = findfirst(==(state), sv)
+state_index(occs, state::Base.RefValue) = state_index(occs, state[])
 
 """
     ManyBodyBasis(b, occupations)
@@ -48,6 +49,7 @@ struct ManyBodyBasis{B,O,UT} <: Basis
 end
 ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(onebodybasis, occupations)
 ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations))
+_vec2fb(mb::ManyBodyBasis) = ManyBodyBasis(mb.onebodybasis, _vec2fb.(mb.occupations))
 
 """
     fermionstates(Nmodes, Nparticles)
@@ -333,7 +335,8 @@ function onebodyexpect_1(op::SparseOpPureType, state)
     result
 end
 
-Base.@propagate_inbounds function state_transition!(buffer::Vector{Int}, occ::Vector{Int}, at_indices, a_indices)
+# Occupations as Vector{Int}
+Base.@propagate_inbounds function state_transition!(buffer, occ::Vector{Int}, at_indices, a_indices)
     any(==(0), (occ[m] for m in a_indices)) && return nothing
     result = 1
     copyto!(buffer, occ)
@@ -349,6 +352,54 @@ Base.@propagate_inbounds function state_transition!(buffer::Vector{Int}, occ::Ve
     return √result
 end
 
+# Occupations as bitstrings (fermions only)
+struct FermionBitstring{T<:Unsigned}
+    bits::T
+    n::Int
+    function FermionBitstring(bits::T, n::Int) where T<:Unsigned
+        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+        nrest = sizeof(bits) * 8 - n
+        new{T}(UInt((bits << nrest) >> nrest), n)
+    end
+end
+FermionBitstring(bits::Integer, n::Int) = FermionBitstring(unsigned(bits), n)
+
+Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
+    fb1.bits == fb2.bits && fb1.n == fb2.n
+Base.isless(fb1::FermionBitstring, fb2::FermionBitstring) =
+    fb1.bits < fb2.bits || fb1.bits == fb2.bits && fb1.n < fb2.n
+
+Base.getindex(fb::FermionBitstring, i::Int) = Bool((fb.bits >> (i - 1)) & 1)
+write_bit(fb::FermionBitstring, i::Int, value::Bool) =
+    value ? FermionBitstring{T}(fb.bits | (one(fb.bits) << (i - 1)), fb.n) :
+            FermionBitstring{T}(fb.bits & ~(one(fb.bits) << (i - 1)), fb.n)
+
+function _vec2fb(occ::Vector{Int})
+    n = length(occ)
+    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
+    bits = UInt(0)
+    for i in 1:n
+        if occ[i] != 0 && occ[i] != 1
+            throw(ArgumentError("Occupations must be 0 or 1"))
+        end
+        occ[i] == 1 && (bits |= UInt(1) << (i - 1))
+    end
+    FermionBitstring(bits, n)
+end
+
+Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstring, at_indices, a_indices)
+    for i in a_indices
+        occ[i] || return nothing
+        occ = write_bit(occ, i, false)
+    end
+    for i in at_indices
+        occ[i] && return nothing
+        occ = write_bit(occ, i, true)
+    end
+    buffer[] = occ
+    return 1
+end
+
 function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
     if index == Nmodes
         occupations[index] = Nparticles

From 6e7818ffaa12a36955e30fc35b05f01acef7dc94 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sat, 13 Apr 2024 23:45:37 +0300
Subject: [PATCH 03/12] Fix bitstrings to be 100% compatible with vectors

---
 src/manybody.jl | 60 +++++++++++++++++++++++++++++--------------------
 1 file changed, 36 insertions(+), 24 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index 1c28f552..d89b17da 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -49,7 +49,6 @@ struct ManyBodyBasis{B,O,UT} <: Basis
 end
 ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(onebodybasis, occupations)
 ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations))
-_vec2fb(mb::ManyBodyBasis) = ManyBodyBasis(mb.onebodybasis, _vec2fb.(mb.occupations))
 
 """
     fermionstates(Nmodes, Nparticles)
@@ -356,35 +355,30 @@ end
 struct FermionBitstring{T<:Unsigned}
     bits::T
     n::Int
-    function FermionBitstring(bits::T, n::Int) where T<:Unsigned
-        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
-        nrest = sizeof(bits) * 8 - n
-        new{T}(UInt((bits << nrest) >> nrest), n)
-    end
+end
+function FermionBitstring(bits::T, n::Int) where T<:Unsigned
+    n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+    nrest = sizeof(bits) * 8 - n
+    FermionBitstring{T}(UInt((bits << nrest) >> nrest), n)
 end
 FermionBitstring(bits::Integer, n::Int) = FermionBitstring(unsigned(bits), n)
 
-Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
+@inline Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
     fb1.bits == fb2.bits && fb1.n == fb2.n
-Base.isless(fb1::FermionBitstring, fb2::FermionBitstring) =
+@inline Base.isless(fb1::FermionBitstring, fb2::FermionBitstring) =
     fb1.bits < fb2.bits || fb1.bits == fb2.bits && fb1.n < fb2.n
+Base.show(io::IO, fb::FermionBitstring{T}) where {T} =
+    print(io, "FermionBitstring{$T}(0b", bitstring(fb.bits)[end-fb.n+1:end], ", ", fb.n, ")")
 
-Base.getindex(fb::FermionBitstring, i::Int) = Bool((fb.bits >> (i - 1)) & 1)
-write_bit(fb::FermionBitstring, i::Int, value::Bool) =
-    value ? FermionBitstring{T}(fb.bits | (one(fb.bits) << (i - 1)), fb.n) :
-            FermionBitstring{T}(fb.bits & ~(one(fb.bits) << (i - 1)), fb.n)
-
-function _vec2fb(occ::Vector{Int})
-    n = length(occ)
-    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
-    bits = UInt(0)
-    for i in 1:n
-        if occ[i] != 0 && occ[i] != 1
-            throw(ArgumentError("Occupations must be 0 or 1"))
-        end
-        occ[i] == 1 && (bits |= UInt(1) << (i - 1))
-    end
-    FermionBitstring(bits, n)
+Base.@propagate_inbounds function Base.getindex(fb::FermionBitstring, i::Int)
+    @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
+    Bool((fb.bits >> (fb.n - i)) & 1)
+end
+Base.@propagate_inbounds function write_bit(fb::FermionBitstring, i::Int, value::Bool)
+    @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
+    offset = fb.n - i
+    value ? FermionBitstring(fb.bits | (one(fb.bits) << offset), fb.n) :
+            FermionBitstring(fb.bits & ~(one(fb.bits) << offset), fb.n)
 end
 
 Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstring, at_indices, a_indices)
@@ -400,6 +394,24 @@ Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstrin
     return 1
 end
 
+# This function exists merely for testing purposes
+function _vec2fb(occ::Vector{Int})
+    n = length(occ)
+    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
+    bits = UInt(0)
+    for i in 1:n
+        if occ[i] != 0 && occ[i] != 1
+            throw(ArgumentError("Occupations must be 0 or 1"))
+        end
+        occ[i] == 1 && (bits |= UInt(1) << (i - 1))
+    end
+    FermionBitstring(bits, n)
+end
+function _vec2fb(mb::ManyBodyBasis)
+    new_sv = SortedVector(_vec2fb.(mb.occupations.sortedvector), mb.occupations.ord)
+    ManyBodyBasis(mb.onebodybasis, new_sv)
+end
+
 function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
     if index == Nmodes
         occupations[index] = Nparticles

From 50a98eed1bbddf220d9d15a724da2e83d3daca6a Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sun, 14 Apr 2024 02:12:25 +0300
Subject: [PATCH 04/12] Implement bitstrings into `fermionstates`, export, add
 docs

---
 src/QuantumOpticsBase.jl |  2 +-
 src/manybody.jl          | 76 +++++++++++++++++++++++++++++-----------
 2 files changed, 57 insertions(+), 21 deletions(-)

diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl
index 2008cf14..8e959490 100644
--- a/src/QuantumOpticsBase.jl
+++ b/src/QuantumOpticsBase.jl
@@ -54,7 +54,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
         #nlevel
                 NLevelBasis, transition, nlevelstate,
         #manybody
-                ManyBodyBasis, fermionstates, bosonstates,
+                ManyBodyBasis, FermionBitstring, fermionstates, bosonstates,
                 manybodyoperator, onebodyexpect,
         #metrics
                 tracenorm, tracenorm_h, tracenorm_nh,
diff --git a/src/manybody.jl b/src/manybody.jl
index d89b17da..dde440d3 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -51,16 +51,22 @@ ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(
 ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations))
 
 """
-    fermionstates(Nmodes, Nparticles)
-    fermionstates(b, Nparticles)
+    fermionstates([T, ]Nmodes, Nparticles)
+    fermionstates([T, ]b, Nparticles)
 
 Generate all fermionic occupation states for N-particles in M-modes.
 `Nparticles` can be a vector to define a Hilbert space with variable
-particle number.
+particle number. `T` is the type of the occupation states - default is `Vector{Int}`,
+but can also be `FermionBitstring`.
 """
-fermionstates(Nmodes::Int, Nparticles::Int) = SortedVector(_distribute_fermions(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]), Base.Reverse)
-fermionstates(Nmodes::Int, Nparticles::Vector{Int}) = union((fermionstates(Nmodes, N) for N in Nparticles)...)
-fermionstates(onebodybasis::Basis, Nparticles) = fermionstates(length(onebodybasis), Nparticles)
+function fermionstates(T::Type, Nmodes::Int, Nparticles::Int)
+    occ_buffer = zero(similar(T, Nmodes))
+    OT = typeof(occ_buffer)
+    SortedVector(_distribute_fermions(Nparticles, Nmodes, 1, occ_buffer, OT[]), Base.Reverse)
+end
+fermionstates(T::Type, Nmodes::Int, Nparticles::Vector{Int}) = union((fermionstates(T, Nmodes, N) for N in Nparticles)...)
+fermionstates(T::Type, onebodybasis::Basis, Nparticles) = fermionstates(T, length(onebodybasis), Nparticles)
+fermionstates(arg1, arg2) = fermionstates(Vector{Int}, arg1, arg2)
 
 """
     bosonstates(Nmodes, Nparticles)
@@ -351,18 +357,45 @@ Base.@propagate_inbounds function state_transition!(buffer, occ::Vector{Int}, at
     return √result
 end
 
-# Occupations as bitstrings (fermions only)
+"""
+    FermionBitstring(bits, n)
+
+Fermionic occupation state represented as a bitstring. The bitstring `bits` is represented
+as an unsigned integer with the `n` least significant bits representing the occupation
+state. The bitstring is assumed to be in the canonical form where the bits are ordered
+from left to right with the most significant bit being the first mode.
+
+---
+
+    FermionBitstring(bits, n)
+
+Create a `FermionBitstring` from an unsigned integer `bits` with the `n` least significant
+bits representing the occupation state. All remaining bits are set to zero.
+"""
 struct FermionBitstring{T<:Unsigned}
     bits::T
     n::Int
-end
-function FermionBitstring(bits::T, n::Int) where T<:Unsigned
-    n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
-    nrest = sizeof(bits) * 8 - n
-    FermionBitstring{T}(UInt((bits << nrest) >> nrest), n)
+    function FermionBitstring{T}(bits, n) where T<:Unsigned
+        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+        new{T}(bits, n)
+    end
+    function FermionBitstring(bits::T, n::Int) where T<:Unsigned
+        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+        nrest = sizeof(bits) * 8 - n
+        FermionBitstring{T}(UInt((bits << nrest) >> nrest), n)
+    end
 end
 FermionBitstring(bits::Integer, n::Int) = FermionBitstring(unsigned(bits), n)
 
+Base.zero(fb::FermionBitstring) = FermionBitstring(zero(fb.bits), fb.n)
+Base.similar(::Type{FermionBitstring{T}}, n::Int) where {T} = FermionBitstring(zero(T), n)
+function Base.similar(::Type{FermionBitstring}, n::Int)
+    for type in (UInt8, UInt16, UInt32, UInt64, UInt128)
+        n ≤ sizeof(type) * 8 && return FermionBitstring(zero(type), n)
+    end
+    throw(ArgumentError("n must be less than 128"))
+end
+Base.copy(fb::FermionBitstring) = fb
 @inline Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
     fb1.bits == fb2.bits && fb1.n == fb2.n
 @inline Base.isless(fb1::FermionBitstring, fb2::FermionBitstring) =
@@ -380,7 +413,10 @@ Base.@propagate_inbounds function write_bit(fb::FermionBitstring, i::Int, value:
     value ? FermionBitstring(fb.bits | (one(fb.bits) << offset), fb.n) :
             FermionBitstring(fb.bits & ~(one(fb.bits) << offset), fb.n)
 end
-
+Base.@propagate_inbounds function write_bit(v::AbstractVector, i::Int, value::Bool)
+    setindex!(v, value, i)
+    return v
+end
 Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstring, at_indices, a_indices)
     for i in a_indices
         occ[i] || return nothing
@@ -429,14 +465,14 @@ function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
     if (Nmodes - index) + 1 < Nparticles
         return results
     end
-    if index == Nmodes
-        occupations[index] = Nparticles
-        push!(results, copy(occupations))
-    else
-        for n = min(1, Nparticles):-1:0
-            occupations[index] = n
-            _distribute_fermions(Nparticles - n, Nmodes, index + 1, occupations, results)
+    for new_index in index:Nmodes - Nparticles + 1
+        occupations = write_bit(occupations, new_index, true)
+        if Nparticles == 1
+            push!(results, copy(occupations))
+        else
+            _distribute_fermions(Nparticles - 1, Nmodes, new_index + 1, occupations, results)
         end
+        occupations = write_bit(occupations, new_index, false)
     end
     return results
 end

From 781d1f1164237ea5ff6423073ecf4e34e4d093a4 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sun, 14 Apr 2024 02:52:09 +0300
Subject: [PATCH 05/12] Add tests

---
 src/manybody.jl       | 18 ------------------
 test/test_manybody.jl | 31 +++++++++++++++++++++++++++++++
 2 files changed, 31 insertions(+), 18 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index dde440d3..304ba5af 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -430,24 +430,6 @@ Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstrin
     return 1
 end
 
-# This function exists merely for testing purposes
-function _vec2fb(occ::Vector{Int})
-    n = length(occ)
-    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
-    bits = UInt(0)
-    for i in 1:n
-        if occ[i] != 0 && occ[i] != 1
-            throw(ArgumentError("Occupations must be 0 or 1"))
-        end
-        occ[i] == 1 && (bits |= UInt(1) << (i - 1))
-    end
-    FermionBitstring(bits, n)
-end
-function _vec2fb(mb::ManyBodyBasis)
-    new_sv = SortedVector(_vec2fb.(mb.occupations.sortedvector), mb.occupations.ord)
-    ManyBodyBasis(mb.onebodybasis, new_sv)
-end
-
 function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
     if index == Nmodes
         occupations[index] = Nparticles
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index 0fb6d056..6ecc0d20 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -34,6 +34,20 @@ b = GenericBasis(Nmodes)
 @test ManyBodyBasis(b, bosonstates(b, 1)) == ManyBodyBasis(b, fermionstates(b, 1))
 @test ManyBodyBasis(b, bosonstates(b, 2)) != ManyBodyBasis(b, fermionstates(b, 2))
 
+function _vec2fb(occ::Vector{Int})
+    n = length(occ)
+    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
+    bits = UInt(0)
+    for i in 1:n
+        if occ[i] != 0 && occ[i] != 1
+            throw(ArgumentError("Occupations must be 0 or 1"))
+        end
+        occ[i] == 1 && (bits |= UInt(1) << (n - i))
+    end
+    FermionBitstring(bits, n)
+end
+@test collect(fermionstates(FermionBitstring, b, 3)) == _vec2fb.(fermionstates(b, 3))
+
 # Test basisstate
 b_mb = ManyBodyBasis(b, bosonstates(b, 2))
 psi_mb = basisstate(b_mb, [2, 0, 0, 0, 0])
@@ -140,6 +154,15 @@ y_ = sparse(y)
 @test 1e-12 > D(Y, manybodyoperator(b, y_))
 @test 1e-12 > D(X + Y, manybodyoperator(b, x_ + y_))
 
+# Same for fermions
+bf1 = ManyBodyBasis(b_single, fermionstates(b_single, [1, 2]))
+bf2 = ManyBodyBasis(b_single, fermionstates(FermionBitstring, b_single, [1, 2]))
+
+X1 = manybodyoperator(bf1, x)
+X2 = manybodyoperator(bf2, x)
+@test X1.data == X2.data
+@test 1e-12 > D(X2, manybodyoperator(bf2, x_))
+
 # Particle-particle interaction operator in second quantization
 x = randoperator(b_single ⊗ b_single)
 y = randoperator(b_single ⊗ b_single)
@@ -155,6 +178,14 @@ y_ = sparse(y)
 @test 1e-12 > D(Y, manybodyoperator(b, y_))
 @test 1e-12 > D(X + Y, manybodyoperator(b, x_ + y_))
 
+# Same for fermions
+X1 = manybodyoperator(bf1, x)
+X2 = manybodyoperator(bf2, x)
+Y2 = manybodyoperator(bf2, y)
+@test X1.data == X2.data
+@test 1e-12 > D(X2, manybodyoperator(bf2, x_))
+@test 1e-12 > D(X2 + Y2, manybodyoperator(bf2, x + y))
+
 # Test one-body expect
 x = randoperator(b_single)
 y = randoperator(b_single)

From 563c1a88eb1254ea06217fdd5299ee655b469cd4 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sun, 14 Apr 2024 18:52:38 +0300
Subject: [PATCH 06/12] Implement `OccupationNumbers` interface, correct
 fermion anticommutators

---
 src/manybody.jl       | 93 +++++++++++++++++++++++++++++++++----------
 test/test_manybody.jl |  2 +-
 2 files changed, 73 insertions(+), 22 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index 304ba5af..b237e97d 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -19,14 +19,14 @@ Base.union(sv1::SortedVector{T}, svs::SortedVector{T}...) where {T} =
     SortedVector(union(sv1.sortedvector, (occ.sortedvector for occ in svs)...))
 
 # Special methods for fast operator construction
-allocate_buffer(sv::AbstractVector) = ismutable(first(sv)) ? similar(first(sv)) : Ref(first(sv))
 function state_index(sv::SortedVector{T}, state::T) where {T}
     ret = searchsortedfirst(sv.sortedvector, state, order = sv.ord)
     ret == length(sv) + 1 && return nothing
     return sv.sortedvector[ret] == state ? ret : nothing
 end
 state_index(sv::AbstractVector{T}, state::T) where {T} = findfirst(==(state), sv)
-state_index(occs, state::Base.RefValue) = state_index(occs, state[])
+state_index(occs::AbstractVector{T}, state::Base.RefValue{T}) where {T} = state_index(occs, state[])
+state_index(sv::AbstractVector{T}, state::Any) where {T} = state_index(sv, convert(T, state))
 
 """
     ManyBodyBasis(b, occupations)
@@ -50,6 +50,9 @@ end
 ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(onebodybasis, occupations)
 ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations))
 
+allocate_buffer(occ) = similar(occ)
+allocate_buffer(mb::ManyBodyBasis) = allocate_buffer(first(mb.occupations))
+
 """
     fermionstates([T, ]Nmodes, Nparticles)
     fermionstates([T, ]b, Nparticles)
@@ -66,7 +69,7 @@ function fermionstates(T::Type, Nmodes::Int, Nparticles::Int)
 end
 fermionstates(T::Type, Nmodes::Int, Nparticles::Vector{Int}) = union((fermionstates(T, Nmodes, N) for N in Nparticles)...)
 fermionstates(T::Type, onebodybasis::Basis, Nparticles) = fermionstates(T, length(onebodybasis), Nparticles)
-fermionstates(arg1, arg2) = fermionstates(Vector{Int}, arg1, arg2)
+fermionstates(arg1, arg2) = fermionstates(OccupationNumbers{FermionStatistics,Int}, arg1, arg2)
 
 """
     bosonstates(Nmodes, Nparticles)
@@ -76,7 +79,11 @@ Generate all bosonic occupation states for N-particles in M-modes.
 `Nparticles` can be a vector to define a Hilbert space with variable
 particle number.
 """
-bosonstates(Nmodes::Int, Nparticles::Int) = SortedVector(_distribute_bosons(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]), Base.Reverse)
+bosonstates(Nmodes::Int, Nparticles::Int) = SortedVector(
+    _distribute_bosons(Nparticles, Nmodes, 1,
+    OccupationNumbers(BosonStatistics(), zeros(Int, Nmodes)),
+    OccupationNumbers{BosonStatistics, Int}[]),
+    Base.Reverse)
 bosonstates(Nmodes::Int, Nparticles::Vector{Int}) = union((bosonstates(Nmodes, N) for N in Nparticles)...)
 bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis), Nparticles)
 
@@ -144,7 +151,7 @@ function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T}
     Is = Int[]
     Js = Int[]
     Vs = T[]
-    buffer = allocate_buffer(b.occupations)
+    buffer = allocate_buffer(b)
     # <{m}_j| at_to a_from |{m}_i>
     for (i, occ_i) in enumerate(b.occupations)
         C = state_transition!(buffer, occ_i, to, from)
@@ -200,7 +207,7 @@ end
 function manybodyoperator_1(basis::ManyBodyBasis, op::Operator)
     S = length(basis.onebodybasis)
     result = DenseOperator(basis)
-    buffer = allocate_buffer(basis.occupations)
+    buffer = allocate_buffer(basis)
     @inbounds for j = 1:S, i = 1:S
         value = op.data[i, j]
         iszero(value) && continue
@@ -221,7 +228,7 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType)
     Is = Int[]
     Js = Int[]
     Vs = ComplexF64[]
-    buffer = allocate_buffer(basis.occupations)
+    buffer = allocate_buffer(basis)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(basis.occupations)
             C = state_transition!(buffer, occ, column, row)
@@ -241,7 +248,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
     @assert S^2 == length(op.basis_l)
     result = DenseOperator(basis)
     op_data = reshape(op.data, S, S, S, S)
-    buffer = allocate_buffer(basis.occupations)
+    buffer = allocate_buffer(basis)
     @inbounds for l = 1:S, k = 1:S, j = 1:S, i = 1:S
         value = op_data[i, j, k, l]
         iszero(value) && continue
@@ -262,7 +269,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
     Is = Int[]
     Js = Int[]
     Vs = ComplexF64[]
-    buffer = allocate_buffer(basis.occupations)
+    buffer = allocate_buffer(basis)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(basis.occupations)
             index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2+row])
@@ -307,7 +314,7 @@ function onebodyexpect_1(op::Operator, state)
     b = basis(state)
     occupations = b.occupations
     S = length(b.onebodybasis)
-    buffer = allocate_buffer(occupations)
+    buffer = allocate_buffer(b)
     result = complex(0.0)
     for i = 1:S, j = 1:S
         value = op.data[i, j]
@@ -326,7 +333,7 @@ end
 function onebodyexpect_1(op::SparseOpPureType, state)
     b = basis(state)
     occupations = b.occupations
-    buffer = allocate_buffer(occupations)
+    buffer = allocate_buffer(b)
     result = complex(0.0)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(occupations)
@@ -341,20 +348,58 @@ function onebodyexpect_1(op::SparseOpPureType, state)
 end
 
 # Occupations as Vector{Int}
-Base.@propagate_inbounds function state_transition!(buffer, occ::Vector{Int}, at_indices, a_indices)
-    any(==(0), (occ[m] for m in a_indices)) && return nothing
-    result = 1
+struct OccupationNumbers{StatisticsT,T} <: AbstractVector{T}
+    statistics::StatisticsT
+    occupations::Vector{T}
+end
+Base.size(on::OccupationNumbers) = size(on.occupations)
+Base.@propagate_inbounds Base.getindex(on::OccupationNumbers, v...) = getindex(on.occupations, v...)
+Base.@propagate_inbounds Base.setindex!(on::OccupationNumbers, value, v...) =
+    setindex!(on.occupations, value, v...)
+Base.IndexStyle(::Type{<:OccupationNumbers}) = Base.IndexLinear()
+Base.similar(occ::OccupationNumbers, ::Type{T}, dims::Dims) where {T} =
+    OccupationNumbers(occ.statistics, similar(occ.occupations, T, dims))
+Base.similar(::Type{OccupationNumbers{StatisticsT,T}}, dims::Dims) where {StatisticsT,T} =
+    OccupationNumbers(StatisticsT(), similar(Vector{T}, dims))
+Base.convert(::Type{OccupationNumbers{StatisticsT,T}}, occ::AbstractVector) where {StatisticsT,T} =
+    OccupationNumbers(StatisticsT(), convert(Vector{T}, occ))
+
+struct FermionStatistics end
+struct BosonStatistics end
+
+Base.@propagate_inbounds function state_transition!(buffer, occ::OccupationNumbers{BosonStatistics},
+        at_indices, a_indices)
+    factor_sq = 1
     copyto!(buffer, occ)
     for i in a_indices
-        result *= buffer[i]
-        result == 0 && return nothing
+        factor_sq *= buffer[i]
+        factor_sq == 0 && return nothing
         buffer[i] -= 1
     end
     for i in at_indices
         buffer[i] += 1
-        result *= buffer[i]
+        factor_sq *= buffer[i]
+    end
+    return √factor_sq
+end
+
+Base.@propagate_inbounds function state_transition!(buffer, occ::OccupationNumbers{FermionStatistics},
+        at_indices, a_indices)
+    factor = 1
+    copyto!(buffer, occ)
+    for i in a_indices
+        buffer[i] == 0 && return nothing
+        buffer[i] = 0
+        ncomm = count(==(1), @view buffer[1:i-1])
+        isodd(ncomm) && (factor *= -1)
+    end
+    for i in at_indices
+        buffer[i] == 1 && return nothing
+        buffer[i] = 1
+        ncomm = count(==(1), @view buffer[1:i-1])
+        isodd(ncomm) && (factor *= -1)
     end
-    return √result
+    return factor
 end
 
 """
@@ -382,7 +427,7 @@ struct FermionBitstring{T<:Unsigned}
     function FermionBitstring(bits::T, n::Int) where T<:Unsigned
         n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
         nrest = sizeof(bits) * 8 - n
-        FermionBitstring{T}(UInt((bits << nrest) >> nrest), n)
+        FermionBitstring{T}(UInt((bits << nrest) >>> nrest), n)
     end
 end
 FermionBitstring(bits::Integer, n::Int) = FermionBitstring(unsigned(bits), n)
@@ -396,6 +441,7 @@ function Base.similar(::Type{FermionBitstring}, n::Int)
     throw(ArgumentError("n must be less than 128"))
 end
 Base.copy(fb::FermionBitstring) = fb
+allocate_buffer(fb::FermionBitstring) = Ref(fb)
 @inline Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
     fb1.bits == fb2.bits && fb1.n == fb2.n
 @inline Base.isless(fb1::FermionBitstring, fb2::FermionBitstring) =
@@ -405,7 +451,7 @@ Base.show(io::IO, fb::FermionBitstring{T}) where {T} =
 
 Base.@propagate_inbounds function Base.getindex(fb::FermionBitstring, i::Int)
     @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
-    Bool((fb.bits >> (fb.n - i)) & 1)
+    Bool((fb.bits >>> (fb.n - i)) & 1)
 end
 Base.@propagate_inbounds function write_bit(fb::FermionBitstring, i::Int, value::Bool)
     @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
@@ -418,16 +464,21 @@ Base.@propagate_inbounds function write_bit(v::AbstractVector, i::Int, value::Bo
     return v
 end
 Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstring, at_indices, a_indices)
+    factor = 1
     for i in a_indices
         occ[i] || return nothing
         occ = write_bit(occ, i, false)
+        ncomm = count_ones(occ.bits >>> (occ.n - i + 1))
+        isodd(ncomm) && (factor *= -1)
     end
     for i in at_indices
         occ[i] && return nothing
         occ = write_bit(occ, i, true)
+        ncomm = count_ones(occ.bits >>> (occ.n - i + 1))
+        isodd(ncomm) && (factor *= -1)
     end
     buffer[] = occ
-    return 1
+    return factor
 end
 
 function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index 6ecc0d20..77cc86aa 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -34,7 +34,7 @@ b = GenericBasis(Nmodes)
 @test ManyBodyBasis(b, bosonstates(b, 1)) == ManyBodyBasis(b, fermionstates(b, 1))
 @test ManyBodyBasis(b, bosonstates(b, 2)) != ManyBodyBasis(b, fermionstates(b, 2))
 
-function _vec2fb(occ::Vector{Int})
+function _vec2fb(occ::AbstractVector{Int})
     n = length(occ)
     n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
     bits = UInt(0)

From e7d30fcc20e0ebe68757a7966428bc8793038aca Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sun, 14 Apr 2024 23:05:27 +0300
Subject: [PATCH 07/12] Add tests

---
 src/manybody.jl       | 20 +++++++++++++++-----
 test/test_manybody.jl | 24 ++++++++++++++++++++++++
 2 files changed, 39 insertions(+), 5 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index b237e97d..cb712a74 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -440,6 +440,16 @@ function Base.similar(::Type{FermionBitstring}, n::Int)
     end
     throw(ArgumentError("n must be less than 128"))
 end
+function Base.convert(::Type{FermionBitstring{T}}, v::AbstractVector) where {T}
+    n = length(v)
+    n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+    bits = zero(T)
+    for i in 1:n
+        v[i] in (0, 1) || throw(ArgumentError("Occupations must be 0 or 1"))
+        v[i] == 1 && (bits |= one(T) << (n - i))
+    end
+    FermionBitstring{T}(bits, n)
+end
 Base.copy(fb::FermionBitstring) = fb
 allocate_buffer(fb::FermionBitstring) = Ref(fb)
 @inline Base.:(==)(fb1::FermionBitstring, fb2::FermionBitstring) =
@@ -498,13 +508,13 @@ function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
     if (Nmodes - index) + 1 < Nparticles
         return results
     end
+    if Nparticles == 0
+        push!(results, copy(occupations))
+        return results
+    end
     for new_index in index:Nmodes - Nparticles + 1
         occupations = write_bit(occupations, new_index, true)
-        if Nparticles == 1
-            push!(results, copy(occupations))
-        else
-            _distribute_fermions(Nparticles - 1, Nmodes, new_index + 1, occupations, results)
-        end
+        _distribute_fermions(Nparticles - 1, Nmodes, new_index + 1, occupations, results)
         occupations = write_bit(occupations, new_index, false)
     end
     return results
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index 77cc86aa..14650a81 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -137,6 +137,30 @@ d3 = destroy(b_mb, 3)
 c2 = create(b_mb, 2)
 @test t22_33 ≈ c2 * c2 * d3 * d3
 
+# Test commutator and anticommutator
+state = basisstate(b_mb, [1, 0, 0, 0])
+a1 = destroy(b_mb, 1)
+at2 = create(b_mb, 2)
+at3 = create(b_mb, 3)
+@test 1e-12 > D(a1 * at2 * state, at2 * a1 * state)
+@test 1e-12 > D(at2 * at3 * state, at3 * at2 * state)
+
+f_mb = ManyBodyBasis(b, fermionstates(b, [0, 1, 2]))
+state = basisstate(f_mb, [1, 0, 0, 0])
+a1 = destroy(f_mb, 1)
+at2 = create(f_mb, 2)
+at3 = create(f_mb, 3)
+@test 1e-12 > D(a1 * at2 * state, -at2 * a1 * state)
+@test 1e-12 > D(at2 * at3 * state, -at3 * at2 * state)
+
+f_mb2 = ManyBodyBasis(b, fermionstates(FermionBitstring, b, [0, 1, 2]))
+state = basisstate(f_mb2, [1, 0, 0, 0])
+a1 = destroy(f_mb2, 1)
+at2 = create(f_mb2, 2)
+at3 = create(f_mb2, 3)
+@test 1e-12 > D(a1 * at2 * state, -at2 * a1 * state)
+@test 1e-12 > D(at2 * at3 * state, -at3 * at2 * state)
+
 # Single particle operator in second quantization
 b_single = GenericBasis(Nmodes)
 b = ManyBodyBasis(b_single, bosonstates(b_single, [1, 2]))

From fd3f1cd508fb2179a18bef664ff0377942082729 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Sun, 14 Apr 2024 23:18:26 +0300
Subject: [PATCH 08/12] Refactor parameter names

---
 src/manybody.jl | 163 ++++++++++++++++++++++++------------------------
 1 file changed, 81 insertions(+), 82 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index cb712a74..036b12c7 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -1,11 +1,11 @@
 struct SortedVector{T, OT} <: AbstractVector{T}
     sortedvector::Vector{T}
-    ord::OT
-    function SortedVector(occ::AbstractVector{T}, ord::OT=Base.Order.Forward) where {T, OT}
-        if issorted(occ, order=ord)
-            new{T, OT}(occ, ord)
+    order::OT
+    function SortedVector(occ::AbstractVector{T}, order::OT=Base.Order.Forward) where {T, OT}
+        if issorted(occ, order=order)
+            new{T, OT}(occ, order)
         else
-            new{T, OT}(sort(occ, order=ord), ord)
+            new{T, OT}(sort(occ, order=order), order)
         end
     end
 end
@@ -16,17 +16,17 @@ Base.@propagate_inbounds function Base.getindex(sv::SortedVector, i::Int)
     return sv.sortedvector[i]
 end
 Base.union(sv1::SortedVector{T}, svs::SortedVector{T}...) where {T} =
-    SortedVector(union(sv1.sortedvector, (occ.sortedvector for occ in svs)...))
+    SortedVector(union(sv1.sortedvector, (occ.sortedvector for occ in svs)...), sv1.order)
 
 # Special methods for fast operator construction
-function state_index(sv::SortedVector{T}, state::T) where {T}
-    ret = searchsortedfirst(sv.sortedvector, state, order = sv.ord)
+function state_index(sv::SortedVector{T}, occ::T) where {T}
+    ret = searchsortedfirst(sv.sortedvector, occ, order = sv.order)
     ret == length(sv) + 1 && return nothing
-    return sv.sortedvector[ret] == state ? ret : nothing
+    return sv.sortedvector[ret] == occ ? ret : nothing
 end
-state_index(sv::AbstractVector{T}, state::T) where {T} = findfirst(==(state), sv)
-state_index(occs::AbstractVector{T}, state::Base.RefValue{T}) where {T} = state_index(occs, state[])
-state_index(sv::AbstractVector{T}, state::Any) where {T} = state_index(sv, convert(T, state))
+state_index(occupations::AbstractVector{T}, state::T) where {T} = findfirst(==(state), occupations)
+state_index(occupations::AbstractVector{T}, state::Base.RefValue{T}) where {T} = state_index(occupations, state[])
+state_index(occupations::AbstractVector{T}, state::Any) where {T} = state_index(occupations, convert(T, state))
 
 """
     ManyBodyBasis(b, occupations)
@@ -90,85 +90,85 @@ bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis),
 ==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash == b2.occupations_hash && b1.onebodybasis == b2.onebodybasis
 
 """
-    basisstate([T=ComplexF64,] b::ManyBodyBasis, occupation::Vector)
+    basisstate([T=ComplexF64,] mb::ManyBodyBasis, occupation::Vector)
 
 Return a ket state where the system is in the state specified by the given
 occupation numbers.
 """
-function basisstate(::Type{T}, basis::ManyBodyBasis, occupation::Vector) where {T}
-    index = state_index(basis.occupations, occupation)
+function basisstate(::Type{T}, mb::ManyBodyBasis, occupation::Vector) where {T}
+    index = state_index(mb.occupations, occupation)
     if isa(index, Nothing)
         throw(ArgumentError("Occupation not included in many-body basis."))
     end
-    basisstate(T, basis, index)
+    basisstate(T, mb, index)
 end
 
 """
-    create([T=ComplexF64,] b::ManyBodyBasis, index)
+    create([T=ComplexF64,] mb::ManyBodyBasis, index)
 
 Creation operator for the i-th mode of the many-body basis `b`.
 """
-create(::Type{T}, b::ManyBodyBasis, index) where {T} = transition(T, b, index, ())
-create(b::ManyBodyBasis, index) = create(ComplexF64, b, index)
+create(::Type{T}, mb::ManyBodyBasis, index) where {T} = transition(T, mb, index, ())
+create(mb::ManyBodyBasis, index) = create(ComplexF64, mb, index)
 
 """
-    destroy([T=ComplexF64,] b::ManyBodyBasis, index)
+    destroy([T=ComplexF64,] mb::ManyBodyBasis, index)
 
 Annihilation operator for the i-th mode of the many-body basis `b`.
 """
-destroy(::Type{T}, b::ManyBodyBasis, index) where {T} = transition(T, b, (), index)
-destroy(b::ManyBodyBasis, index) = destroy(ComplexF64, b, index)
+destroy(::Type{T}, mb::ManyBodyBasis, index) where {T} = transition(T, mb, (), index)
+destroy(mb::ManyBodyBasis, index) = destroy(ComplexF64, mb, index)
 
 """
-    number([T=ComplexF64,] b::ManyBodyBasis, index)
+    number([T=ComplexF64,] mb::ManyBodyBasis, index)
 
 Particle number operator for the i-th mode of the many-body basis `b`.
 """
-function number(::Type{T}, b::ManyBodyBasis, index) where {T}
-    diagonaloperator(b, T[occ[index] for occ in b.occupations])
+function number(::Type{T}, mb::ManyBodyBasis, index) where {T}
+    diagonaloperator(mb, T[occ[index] for occ in mb.occupations])
 end
-number(b::ManyBodyBasis, index) = number(ComplexF64, b, index)
+number(mb::ManyBodyBasis, index) = number(ComplexF64, mb, index)
 
 """
-    number([T=ComplexF64,] b::ManyBodyBasis)
+    number([T=ComplexF64,] mb::ManyBodyBasis)
 
 Total particle number operator.
 """
-function number(::Type{T}, b::ManyBodyBasis) where {T}
-    diagonaloperator(b, T[sum(occ) for occ in b.occupations])
+function number(::Type{T}, mb::ManyBodyBasis) where {T}
+    diagonaloperator(mb, T[sum(occ) for occ in mb.occupations])
 end
-number(b::ManyBodyBasis) = number(ComplexF64, b)
+number(mb::ManyBodyBasis) = number(ComplexF64, mb)
 
 """
-    transition([T=ComplexF64,] b::ManyBodyBasis, to, from)
+    transition([T=ComplexF64,] mb::ManyBodyBasis, to, from)
 
 Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
 
 Note that `to` and `from` can be collections of indices. The resulting operator in this case
 will be equal to ``a^\\dagger_{to_1} a^\\dagger_{to_2} \\ldots a_{from_2} a_{from_1}``.
 """
-function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T}
+function transition(::Type{T}, mb::ManyBodyBasis, to, from) where {T}
     Is = Int[]
     Js = Int[]
     Vs = T[]
-    buffer = allocate_buffer(b)
+    buffer = allocate_buffer(mb)
     # <{m}_j| at_to a_from |{m}_i>
-    for (i, occ_i) in enumerate(b.occupations)
+    for (i, occ_i) in enumerate(mb.occupations)
         C = state_transition!(buffer, occ_i, to, from)
         C === nothing && continue
-        j = state_index(b.occupations, buffer)
+        j = state_index(mb.occupations, buffer)
         j === nothing && continue
         push!(Is, j)
         push!(Js, i)
         push!(Vs, C)
     end
-    return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
+    return SparseOperator(mb, sparse(Is, Js, Vs, length(mb), length(mb)))
 end
-transition(b::ManyBodyBasis, to, from) = transition(ComplexF64, b, to, from)
+transition(mb::ManyBodyBasis, to, from) = transition(ComplexF64, mb, to, from)
 
 # Calculate many-Body operator from one-body operator
 """
-    manybodyoperator(b::ManyBodyBasis, op)
+    manybodyoperator(mb::ManyBodyBasis, op)
 
 Create the many-body operator from the given one-body operator `op`.
 
@@ -192,70 +192,69 @@ where ``X`` is the N-particle operator, ``x`` is the one-body operator and
 ``|u⟩`` are the one-body states associated to the
 different modes of the N-particle basis.
 """
-function manybodyoperator(basis::ManyBodyBasis, op)
+function manybodyoperator(mb::ManyBodyBasis, op)
     @assert op.basis_l == op.basis_r
-    if op.basis_l == basis.onebodybasis
-        result = manybodyoperator_1(basis, op)
-    elseif op.basis_l == basis.onebodybasis ⊗ basis.onebodybasis
-        result = manybodyoperator_2(basis, op)
+    if op.basis_l == mb.onebodybasis
+        result = manybodyoperator_1(mb, op)
+    elseif op.basis_l == mb.onebodybasis ⊗ mb.onebodybasis
+        result = manybodyoperator_2(mb, op)
     else
         throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis."))
     end
     result
 end
 
-function manybodyoperator_1(basis::ManyBodyBasis, op::Operator)
-    S = length(basis.onebodybasis)
-    result = DenseOperator(basis)
-    buffer = allocate_buffer(basis)
+function manybodyoperator_1(mb::ManyBodyBasis, op::Operator)
+    S = length(mb.onebodybasis)
+    result = DenseOperator(mb)
+    buffer = allocate_buffer(mb)
     @inbounds for j = 1:S, i = 1:S
         value = op.data[i, j]
         iszero(value) && continue
-        for (m, occ) in enumerate(basis.occupations)
+        for (m, occ) in enumerate(mb.occupations)
             C = state_transition!(buffer, occ, j, i)
             C === nothing && continue
-            n = state_index(basis.occupations, buffer)
+            n = state_index(mb.occupations, buffer)
             n === nothing && continue
             result.data[m, n] += C * value
         end
     end
     return result
 end
-manybodyoperator_1(basis::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyoperator_1(basis, dagger(op)))
+manybodyoperator_1(mb::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyoperator_1(mb, dagger(op)))
 
-function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType)
-    N = length(basis)
+function manybodyoperator_1(mb::ManyBodyBasis, op::SparseOpPureType)
+    N = length(mb)
     Is = Int[]
     Js = Int[]
     Vs = ComplexF64[]
-    buffer = allocate_buffer(basis)
+    buffer = allocate_buffer(mb)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
-        for (m, occ) in enumerate(basis.occupations)
+        for (m, occ) in enumerate(mb.occupations)
             C = state_transition!(buffer, occ, column, row)
             C === nothing && continue
-            n = state_index(basis.occupations, buffer)
+            n = state_index(mb.occupations, buffer)
             n === nothing && continue
             push!(Is, m)
             push!(Js, n)
             push!(Vs, C * value)
         end
     end
-    return SparseOperator(basis, sparse(Is, Js, Vs, N, N))
+    return SparseOperator(mb, sparse(Is, Js, Vs, N, N))
 end
 
-function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
-    S = length(basis.onebodybasis)
-    @assert S^2 == length(op.basis_l)
-    result = DenseOperator(basis)
+function manybodyoperator_2(mb::ManyBodyBasis, op::Operator)
+    S = length(mb.onebodybasis)
+    result = DenseOperator(mb)
     op_data = reshape(op.data, S, S, S, S)
-    buffer = allocate_buffer(basis)
+    buffer = allocate_buffer(mb)
     @inbounds for l = 1:S, k = 1:S, j = 1:S, i = 1:S
         value = op_data[i, j, k, l]
         iszero(value) && continue
-        for (m, occ) in enumerate(basis.occupations)
+        for (m, occ) in enumerate(mb.occupations)
             C = state_transition!(buffer, occ, (k, l), (i, j))
             C === nothing && continue
-            n = state_index(basis.occupations, buffer)
+            n = state_index(mb.occupations, buffer)
             n === nothing && continue
             result.data[m, n] += C * value
         end
@@ -263,26 +262,26 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
     return result
 end
 
-function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
-    N = length(basis)
-    S = length(basis.onebodybasis)
+function manybodyoperator_2(mb::ManyBodyBasis, op::SparseOpType)
+    N = length(mb)
+    S = length(mb.onebodybasis)
     Is = Int[]
     Js = Int[]
     Vs = ComplexF64[]
-    buffer = allocate_buffer(basis)
+    buffer = allocate_buffer(mb)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
-        for (m, occ) in enumerate(basis.occupations)
+        for (m, occ) in enumerate(mb.occupations)
             index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2+row])
             C = state_transition!(buffer, occ, index[3:4], index[1:2])
             C === nothing && continue
-            n = state_index(basis.occupations, buffer)
+            n = state_index(mb.occupations, buffer)
             n === nothing && continue
             push!(Is, m)
             push!(Js, n)
             push!(Vs, C * value)
         end
     end
-    return SparseOperator(basis, sparse(Is, Js, Vs, N, N))
+    return SparseOperator(mb, sparse(Is, Js, Vs, N, N))
 end
 
 
@@ -308,13 +307,13 @@ end
 
 onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state = states]
 
-get_value(state::Ket, m, n) = conj(state.data[m]) * state.data[n]
-get_value(state::Operator, m, n) = state.data[n, m]
+matrix_element(state::Ket, m, n) = conj(state.data[m]) * state.data[n]
+matrix_element(state::Operator, m, n) = state.data[n, m]
 function onebodyexpect_1(op::Operator, state)
-    b = basis(state)
-    occupations = b.occupations
-    S = length(b.onebodybasis)
-    buffer = allocate_buffer(b)
+    mb = basis(state)
+    occupations = mb.occupations
+    S = length(mb.onebodybasis)
+    buffer = allocate_buffer(mb)
     result = complex(0.0)
     for i = 1:S, j = 1:S
         value = op.data[i, j]
@@ -324,16 +323,16 @@ function onebodyexpect_1(op::Operator, state)
             C === nothing && continue
             n = state_index(occupations, buffer)
             n === nothing && continue
-            result += C * value * get_value(state, m, n)
+            result += C * value * matrix_element(state, m, n)
         end
     end
     result
 end
 
 function onebodyexpect_1(op::SparseOpPureType, state)
-    b = basis(state)
-    occupations = b.occupations
-    buffer = allocate_buffer(b)
+    mb = basis(state)
+    occupations = mb.occupations
+    buffer = allocate_buffer(mb)
     result = complex(0.0)
     @inbounds for (row, column, value) in zip(findnz(op.data)...)
         for (m, occ) in enumerate(occupations)
@@ -341,7 +340,7 @@ function onebodyexpect_1(op::SparseOpPureType, state)
             C === nothing && continue
             n = state_index(occupations, buffer)
             n === nothing && continue
-            result += C * value * get_value(state, m, n)
+            result += C * value * matrix_element(state, m, n)
         end
     end
     result
@@ -438,7 +437,7 @@ function Base.similar(::Type{FermionBitstring}, n::Int)
     for type in (UInt8, UInt16, UInt32, UInt64, UInt128)
         n ≤ sizeof(type) * 8 && return FermionBitstring(zero(type), n)
     end
-    throw(ArgumentError("n must be less than 128"))
+    throw(ArgumentError("n must be less than 128; got $n"))
 end
 function Base.convert(::Type{FermionBitstring{T}}, v::AbstractVector) where {T}
     n = length(v)

From 7a28e0e59b78369260b532bded59b2b37d86a3d2 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Mon, 15 Apr 2024 19:32:47 +0300
Subject: [PATCH 09/12] Refactor names once again, add `BigInt` support

---
 src/manybody.jl       | 75 +++++++++++++++++++------------------------
 test/test_manybody.jl | 18 +++--------
 2 files changed, 38 insertions(+), 55 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index 036b12c7..56a7a098 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -24,9 +24,9 @@ function state_index(sv::SortedVector{T}, occ::T) where {T}
     ret == length(sv) + 1 && return nothing
     return sv.sortedvector[ret] == occ ? ret : nothing
 end
-state_index(occupations::AbstractVector{T}, state::T) where {T} = findfirst(==(state), occupations)
-state_index(occupations::AbstractVector{T}, state::Base.RefValue{T}) where {T} = state_index(occupations, state[])
-state_index(occupations::AbstractVector{T}, state::Any) where {T} = state_index(occupations, convert(T, state))
+state_index(occupations::AbstractVector{T}, occ::T) where {T} = findfirst(==(occ), occupations)
+state_index(occupations::AbstractVector{T}, occ::Base.RefValue{T}) where {T} = state_index(occupations, occ[])
+state_index(occupations::AbstractVector{T}, occ::Any) where {T} = state_index(occupations, convert(T, occ))
 
 """
     ManyBodyBasis(b, occupations)
@@ -145,7 +145,7 @@ number(mb::ManyBodyBasis) = number(ComplexF64, mb)
 Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
 
 Note that `to` and `from` can be collections of indices. The resulting operator in this case
-will be equal to ``a^\\dagger_{to_1} a^\\dagger_{to_2} \\ldots a_{from_2} a_{from_1}``.
+will be equal to ``\\ldots a^\\dagger_{to_2} a^\\dagger_{to_1} \\ldots a_{from_2} a_{from_1}``.
 """
 function transition(::Type{T}, mb::ManyBodyBasis, to, from) where {T}
     Is = Int[]
@@ -153,8 +153,8 @@ function transition(::Type{T}, mb::ManyBodyBasis, to, from) where {T}
     Vs = T[]
     buffer = allocate_buffer(mb)
     # <{m}_j| at_to a_from |{m}_i>
-    for (i, occ_i) in enumerate(mb.occupations)
-        C = state_transition!(buffer, occ_i, to, from)
+    for (i, occ) in enumerate(mb.occupations)
+        C = state_transition!(buffer, occ, to, from)
         C === nothing && continue
         j = state_index(mb.occupations, buffer)
         j === nothing && continue
@@ -402,36 +402,30 @@ Base.@propagate_inbounds function state_transition!(buffer, occ::OccupationNumbe
 end
 
 """
-    FermionBitstring(bits, n)
+    FermionBitstring{T}
 
-Fermionic occupation state represented as a bitstring. The bitstring `bits` is represented
-as an unsigned integer with the `n` least significant bits representing the occupation
-state. The bitstring is assumed to be in the canonical form where the bits are ordered
-from left to right with the most significant bit being the first mode.
+Bitstring representation of a fermionic occupation state.
 
 ---
 
     FermionBitstring(bits, n)
-
-Create a `FermionBitstring` from an unsigned integer `bits` with the `n` least significant
-bits representing the occupation state. All remaining bits are set to zero.
 """
-struct FermionBitstring{T<:Unsigned}
+struct FermionBitstring{T}
     bits::T
     n::Int
-    function FermionBitstring{T}(bits, n) where T<:Unsigned
-        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
-        new{T}(bits, n)
-    end
-    function FermionBitstring(bits::T, n::Int) where T<:Unsigned
-        n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
-        nrest = sizeof(bits) * 8 - n
-        FermionBitstring{T}(UInt((bits << nrest) >>> nrest), n)
+    function FermionBitstring{T}(bits, n::Integer) where {T}
+        T<:Unsigned && n > sizeof(T) * 8 &&
+            throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
+        mask = T(1) << n - 1
+        new{T}(bits & mask, Int(n))
     end
 end
-FermionBitstring(bits::Integer, n::Int) = FermionBitstring(unsigned(bits), n)
+function FermionBitstring(bits::T, n::Integer) where T
+    FermionBitstring{T}(bits, n)
+end
 
 Base.zero(fb::FermionBitstring) = FermionBitstring(zero(fb.bits), fb.n)
+Base.length(fb::FermionBitstring) = fb.n
 Base.similar(::Type{FermionBitstring{T}}, n::Int) where {T} = FermionBitstring(zero(T), n)
 function Base.similar(::Type{FermionBitstring}, n::Int)
     for type in (UInt8, UInt16, UInt32, UInt64, UInt128)
@@ -439,15 +433,12 @@ function Base.similar(::Type{FermionBitstring}, n::Int)
     end
     throw(ArgumentError("n must be less than 128; got $n"))
 end
-function Base.convert(::Type{FermionBitstring{T}}, v::AbstractVector) where {T}
-    n = length(v)
-    n > sizeof(T) * 8 && throw(ArgumentError("n must be less than $(sizeof(T) * 8)"))
-    bits = zero(T)
-    for i in 1:n
-        v[i] in (0, 1) || throw(ArgumentError("Occupations must be 0 or 1"))
-        v[i] == 1 && (bits |= one(T) << (n - i))
+function Base.convert(::Type{T}, v::AbstractVector) where {T<:FermionBitstring}
+    new_v = similar(T, length(v))
+    for i in 1:length(new_v)
+        new_v = Base.setindex(new_v, Bool(v[i]), i)
     end
-    FermionBitstring{T}(bits, n)
+    return new_v
 end
 Base.copy(fb::FermionBitstring) = fb
 allocate_buffer(fb::FermionBitstring) = Ref(fb)
@@ -462,27 +453,23 @@ Base.@propagate_inbounds function Base.getindex(fb::FermionBitstring, i::Int)
     @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
     Bool((fb.bits >>> (fb.n - i)) & 1)
 end
-Base.@propagate_inbounds function write_bit(fb::FermionBitstring, i::Int, value::Bool)
+Base.@propagate_inbounds function Base.setindex(fb::FermionBitstring, value::Bool, i::Int)
     @boundscheck i in 1:fb.n || throw(BoundsError(fb, i))
     offset = fb.n - i
     value ? FermionBitstring(fb.bits | (one(fb.bits) << offset), fb.n) :
             FermionBitstring(fb.bits & ~(one(fb.bits) << offset), fb.n)
 end
-Base.@propagate_inbounds function write_bit(v::AbstractVector, i::Int, value::Bool)
-    setindex!(v, value, i)
-    return v
-end
 Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstring, at_indices, a_indices)
     factor = 1
     for i in a_indices
         occ[i] || return nothing
-        occ = write_bit(occ, i, false)
+        occ = setocc(occ, i, false)
         ncomm = count_ones(occ.bits >>> (occ.n - i + 1))
         isodd(ncomm) && (factor *= -1)
     end
     for i in at_indices
         occ[i] && return nothing
-        occ = write_bit(occ, i, true)
+        occ = setocc(occ, i, true)
         ncomm = count_ones(occ.bits >>> (occ.n - i + 1))
         isodd(ncomm) && (factor *= -1)
     end
@@ -502,7 +489,11 @@ function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
     end
     return results
 end
-
+Base.@propagate_inbounds setocc(fb::FermionBitstring, i::Int, value::Bool) = Base.setindex(fb, value, i)
+Base.@propagate_inbounds function setocc(v::AbstractVector, i::Int, value::Bool)
+    setindex!(v, value, i)
+    return v
+end
 function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
     if (Nmodes - index) + 1 < Nparticles
         return results
@@ -512,9 +503,9 @@ function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
         return results
     end
     for new_index in index:Nmodes - Nparticles + 1
-        occupations = write_bit(occupations, new_index, true)
+        occupations = setocc(occupations, new_index, true)
         _distribute_fermions(Nparticles - 1, Nmodes, new_index + 1, occupations, results)
-        occupations = write_bit(occupations, new_index, false)
+        occupations = setocc(occupations, new_index, false)
     end
     return results
 end
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index 14650a81..d42fcca6 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -34,19 +34,11 @@ b = GenericBasis(Nmodes)
 @test ManyBodyBasis(b, bosonstates(b, 1)) == ManyBodyBasis(b, fermionstates(b, 1))
 @test ManyBodyBasis(b, bosonstates(b, 2)) != ManyBodyBasis(b, fermionstates(b, 2))
 
-function _vec2fb(occ::AbstractVector{Int})
-    n = length(occ)
-    n > sizeof(UInt) * 8 && throw(ArgumentError("n must be less than $(sizeof(UInt) * 8)"))
-    bits = UInt(0)
-    for i in 1:n
-        if occ[i] != 0 && occ[i] != 1
-            throw(ArgumentError("Occupations must be 0 or 1"))
-        end
-        occ[i] == 1 && (bits |= UInt(1) << (n - i))
-    end
-    FermionBitstring(bits, n)
-end
-@test collect(fermionstates(FermionBitstring, b, 3)) == _vec2fb.(fermionstates(b, 3))
+@test collect(fermionstates(FermionBitstring, b, 3)) ==
+    convert.(FermionBitstring, fermionstates(b, 3))
+b2 = GenericBasis(130)
+@test collect(fermionstates(FermionBitstring{BigInt}, b2, 2)) ==
+    convert.(FermionBitstring{BigInt}, fermionstates(b2, 2))
 
 # Test basisstate
 b_mb = ManyBodyBasis(b, bosonstates(b, 2))

From 2a4d7fca181646224d8de9d25849a5ef464381a3 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Mon, 15 Apr 2024 19:49:33 +0300
Subject: [PATCH 10/12] Performance

---
 src/manybody.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/manybody.jl b/src/manybody.jl
index 56a7a098..33540702 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -360,6 +360,8 @@ Base.similar(occ::OccupationNumbers, ::Type{T}, dims::Dims) where {T} =
     OccupationNumbers(occ.statistics, similar(occ.occupations, T, dims))
 Base.similar(::Type{OccupationNumbers{StatisticsT,T}}, dims::Dims) where {StatisticsT,T} =
     OccupationNumbers(StatisticsT(), similar(Vector{T}, dims))
+Base.isless(occ1::OccupationNumbers, occ2::OccupationNumbers) = occ1.occupations < occ2.occupations
+Base.copyto!(dest::OccupationNumbers, src::OccupationNumbers) = copyto!(dest.occupations, src.occupations)
 Base.convert(::Type{OccupationNumbers{StatisticsT,T}}, occ::AbstractVector) where {StatisticsT,T} =
     OccupationNumbers(StatisticsT(), convert(Vector{T}, occ))
 
@@ -368,6 +370,7 @@ struct BosonStatistics end
 
 Base.@propagate_inbounds function state_transition!(buffer, occ::OccupationNumbers{BosonStatistics},
         at_indices, a_indices)
+    any(==(0), (occ[m] for m in a_indices)) && return nothing
     factor_sq = 1
     copyto!(buffer, occ)
     for i in a_indices
@@ -384,6 +387,7 @@ end
 
 Base.@propagate_inbounds function state_transition!(buffer, occ::OccupationNumbers{FermionStatistics},
         at_indices, a_indices)
+    any(==(0), (occ[m] for m in a_indices)) && return nothing
     factor = 1
     copyto!(buffer, occ)
     for i in a_indices

From 94b58eba7dbdc99b44822fe1fabeeacc1dc63b74 Mon Sep 17 00:00:00 2001
From: "A. Yavorsky" <iavorskii.ar@phystech.edu>
Date: Tue, 21 May 2024 23:10:07 +0300
Subject: [PATCH 11/12] Add occupations type for `bosonstates` constructor

---
 src/manybody.jl       | 40 +++++++++++++++++++++-------------------
 test/test_manybody.jl |  2 ++
 2 files changed, 23 insertions(+), 19 deletions(-)

diff --git a/src/manybody.jl b/src/manybody.jl
index 33540702..f36550a3 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -59,8 +59,8 @@ allocate_buffer(mb::ManyBodyBasis) = allocate_buffer(first(mb.occupations))
 
 Generate all fermionic occupation states for N-particles in M-modes.
 `Nparticles` can be a vector to define a Hilbert space with variable
-particle number. `T` is the type of the occupation states - default is `Vector{Int}`,
-but can also be `FermionBitstring`.
+particle number. `T` is the type of the occupation states - default is
+`OccupationNumbers{FermionStatistics,Int}`, but can be any occupations type.
 """
 function fermionstates(T::Type, Nmodes::Int, Nparticles::Int)
     occ_buffer = zero(similar(T, Nmodes))
@@ -72,20 +72,22 @@ fermionstates(T::Type, onebodybasis::Basis, Nparticles) = fermionstates(T, lengt
 fermionstates(arg1, arg2) = fermionstates(OccupationNumbers{FermionStatistics,Int}, arg1, arg2)
 
 """
-    bosonstates(Nmodes, Nparticles)
-    bosonstates(b, Nparticles)
+    bosonstates([T, ]Nmodes, Nparticles)
+    bosonstates([T, ]b, Nparticles)
 
 Generate all bosonic occupation states for N-particles in M-modes.
 `Nparticles` can be a vector to define a Hilbert space with variable
-particle number.
+particle number. `T` is the type of the occupation states - default is
+`OccupationNumbers{BosonStatistics,Int}`, but can be any occupations type.
 """
-bosonstates(Nmodes::Int, Nparticles::Int) = SortedVector(
-    _distribute_bosons(Nparticles, Nmodes, 1,
-    OccupationNumbers(BosonStatistics(), zeros(Int, Nmodes)),
-    OccupationNumbers{BosonStatistics, Int}[]),
-    Base.Reverse)
-bosonstates(Nmodes::Int, Nparticles::Vector{Int}) = union((bosonstates(Nmodes, N) for N in Nparticles)...)
-bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis), Nparticles)
+function bosonstates(T::Type, Nmodes::Int, Nparticles::Int)
+    occ_buffer = zero(similar(T, Nmodes))
+    OT = typeof(occ_buffer)
+    SortedVector(_distribute_bosons(Nparticles, Nmodes, 1, occ_buffer, OT[]), Base.Reverse)
+end
+bosonstates(T::Type, Nmodes::Int, Nparticles::Vector{Int}) = union((bosonstates(T, Nmodes, N) for N in Nparticles)...)
+bosonstates(T::Type, onebodybasis::Basis, Nparticles) = bosonstates(T, length(onebodybasis), Nparticles)
+bosonstates(arg1, arg2) = bosonstates(OccupationNumbers{BosonStatistics,Int}, arg1, arg2)
 
 ==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash == b2.occupations_hash && b1.onebodybasis == b2.onebodybasis
 
@@ -481,23 +483,23 @@ Base.@propagate_inbounds function state_transition!(buffer, occ::FermionBitstrin
     return factor
 end
 
+Base.@propagate_inbounds setocc(fb::FermionBitstring, i::Int, value) = Base.setindex(fb, value, i)
+Base.@propagate_inbounds function setocc(v::AbstractVector, i::Int, value)
+    setindex!(v, value, i)
+    return v
+end
 function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
     if index == Nmodes
-        occupations[index] = Nparticles
+        occupations = setocc(occupations, index, Nparticles)
         push!(results, copy(occupations))
     else
         for n = Nparticles:-1:0
-            occupations[index] = n
+            occupations = setocc(occupations, index, n)
             _distribute_bosons(Nparticles - n, Nmodes, index + 1, occupations, results)
         end
     end
     return results
 end
-Base.@propagate_inbounds setocc(fb::FermionBitstring, i::Int, value::Bool) = Base.setindex(fb, value, i)
-Base.@propagate_inbounds function setocc(v::AbstractVector, i::Int, value::Bool)
-    setindex!(v, value, i)
-    return v
-end
 function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
     if (Nmodes - index) + 1 < Nparticles
         return results
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index d42fcca6..137243a5 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -36,6 +36,8 @@ b = GenericBasis(Nmodes)
 
 @test collect(fermionstates(FermionBitstring, b, 3)) ==
     convert.(FermionBitstring, fermionstates(b, 3))
+@test collect(bosonstates(Vector{Int}, b, 3)) ==
+    convert.(Vector{Int}, bosonstates(b, 3))
 b2 = GenericBasis(130)
 @test collect(fermionstates(FermionBitstring{BigInt}, b2, 2)) ==
     convert.(FermionBitstring{BigInt}, fermionstates(b2, 2))

From 4bb56d76ac0fe2a2f6616ef04c9e22e52b546b13 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <github.acc@krastanov.org>
Date: Mon, 10 Jun 2024 14:36:50 -0400
Subject: [PATCH 12/12] bump version number

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a3216b37..96c20987 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.4.21"
+version = "0.4.22"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"