29 #pragma implementation
32 #ifndef _mpqc_src_lib_chemistry_qc_nbody_string_h
33 #define _mpqc_src_lib_chemistry_qc_nbody_string_h
43 #include <unordered_set>
45 #define BOOST_DYNAMIC_BITSET_DONT_USE_FRIENDS // need this to expose dynamic_bitset's naughty bits
46 #include <boost/dynamic_bitset.hpp>
47 #include <boost/functional/hash.hpp>
55 template <
size_t Ns=64ul>
58 typedef std::bitset<Ns> parent_type;
59 typedef size_t state_index_type;
73 for(
auto v : occupied_states) {
91 std::bitset<Ns>::reset();
107 return std::bitset<Ns>::count();
116 (*this)[from] =
false;
136 size_t count(
size_t pos1,
size_t pos2)
const {
140 tmp <<= (Ns - (pos2 - pos1 - 1));
144 size_t hash_value()
const {
145 std::hash<std::bitset<Ns> > h;
149 template <
size_t M>
bool operator==(
const FermionOccupationNBitString<M>&
other)
const {
150 return nstates_ ==
other.nstates_ &&
151 static_cast<parent_type>(*
this) == static_cast<parent_type>(
other);
155 MPQC_ASSERT(nstates_+
other.nstates_ <= Ns);
158 for(
size_t pos1=0; pos1<nstates_; ++pos1, ++pos)
159 result[pos] = (*
this)[pos1];
160 for(
size_t pos2=0; pos2<
other.nstates_; ++pos2, ++pos)
161 result[pos] =
other[pos2];
169 template <
size_t N1,
size_t N2>
170 FermionOccupationNBitString<N1+N2>
operator+(
const FermionOccupationNBitString<N1>& s1,
171 const FermionOccupationNBitString<N2>& s2) {
172 FermionOccupationNBitString<N1+N2> result;
174 for(
size_t pos1=0; pos1<N1; ++pos1, ++pos)
175 result[pos] = s1[pos1];
176 for(
size_t pos2=0; pos2<N2; ++pos2, ++pos)
177 result[pos] = s2[pos2];
181 template <
class CharT,
class Traits,
size_t Ns>
182 std::basic_ostream<CharT, Traits>&
183 operator<<(std::basic_ostream<CharT, Traits>& os,
184 const FermionOccupationNBitString<Ns>& x)
186 os << std::bitset<Ns>(x);
196 typedef boost::dynamic_bitset<> parent_type;
197 typedef size_t state_index_type;
213 boost::dynamic_bitset<>(N, 0ul)
215 for(
auto v : occupied_states) {
233 boost::dynamic_bitset<>::reset();
241 return boost::dynamic_bitset<>::count();
250 (*this)[from] =
false;
270 size_t count(
size_t pos1,
size_t pos2)
const {
274 tmp <<= (size() - (pos2 - pos1 - 1));
278 size_t hash_value()
const {
279 typedef std::vector<parent_type::block_type, parent_type::allocator_type> m_bits_type;
280 boost::hash<m_bits_type> h;
285 return size() ==
other.size() &&
286 static_cast<parent_type>(*
this) == static_cast<parent_type>(
other);
292 FermionOccupationDBitString
operator+(
const FermionOccupationDBitString& s1,
293 const FermionOccupationDBitString& s2) {
294 const size_t s1_size = s1.size();
295 const size_t s2_size = s2.size();
296 FermionOccupationDBitString result1(s1);
297 FermionOccupationDBitString result2(s2);
298 result1.resize(s1_size + s2_size);
299 result2.resize(s1_size + s2_size);
301 return result1 | result2;
308 template <
class CharT,
class Traits>
309 std::basic_ostream<CharT, Traits>&
310 operator<<(std::basic_ostream<CharT, Traits>& os,
311 const FermionOccupationDBitString& x)
313 os << boost::dynamic_bitset<>(x);
322 typedef size_t state_index_type;
326 Block(
size_t o,
bool n =
true) : offset(o), length(1ul), occ(n) {}
327 Block(
size_t o,
size_t l,
bool n =
true) : offset(o), length(l), occ(n) {}
331 bool operator<(
const Block&
other)
const {
return offset <
other.offset; }
333 operator long int()
const {
334 long int result = offset;
337 return occ ? -result : result;
340 typedef std::set<Block> Blocks;
347 blocks_.insert(
Block(0,Ns,
false));
355 blocks_.insert(
Block(0,Ns,
false));
357 for(
auto v : occupied_states) {
367 std::swap(blocks_,blocks);
377 return this->
count() == 0;
385 blocks_.insert(
Block(0,nstates_,
false));
402 for(
auto v : blocks_) {
403 if (v.occ) c += v.length;
414 auto ub_iter = blocks_.upper_bound(
Block(i));
415 auto result_iter = --ub_iter;
416 return result_iter->occ;
430 auto ub_iter = blocks_.upper_bound(
Block(from));
431 auto result_iter = --ub_iter;
432 MPQC_ASSERT(result_iter->occ ==
true);
434 if (result_iter->length == 1) {
435 auto uocc_block = this->block_erase(result_iter);
437 else if (from == result_iter->offset) {
438 auto occ_block = this->block_pop_front(result_iter);
439 if (occ_block != blocks_.begin()) {
440 auto uocc_block = this->block_push_back(--occ_block);
443 auto uocc_block = blocks_.insert(occ_block,
Block(0, 1,
false));
445 else if (from == result_iter->offset + result_iter->length - 1) {
446 auto occ_block = this->block_pop_back(result_iter);
447 if (++occ_block != blocks_.end()) {
448 auto uocc_block = this->block_push_front(occ_block);
451 auto uocc_block = blocks_.insert(occ_block,
Block(0, 1,
false));
454 this->block_split(result_iter, from);
471 auto ub_iter = blocks_.upper_bound(
Block(to));
472 auto result_iter = --ub_iter;
473 MPQC_ASSERT(result_iter->occ ==
false);
475 if (result_iter->length == 1) {
476 auto occ_block = this->block_replace(result_iter,
Block(to,1,
true));
478 else if (to == result_iter->offset) {
479 auto uocc_block = this->block_pop_front(result_iter);
480 if (uocc_block != blocks_.begin()) {
481 auto occ_block = this->block_push_back(--uocc_block);
484 auto occ_block = blocks_.insert(uocc_block,
Block(0, 1,
true));
486 else if (to == result_iter->offset + result_iter->length - 1) {
487 auto uocc_block = this->block_pop_back(result_iter);
488 if (++uocc_block != blocks_.end()) {
489 auto occ_block = this->block_push_front(uocc_block);
492 auto occ_block = blocks_.insert(uocc_block,
Block(nstates_-1, 1,
true));
495 this->block_split(result_iter, to);
507 size_t count(
size_t pos1,
size_t pos2)
const {
509 MPQC_ASSERT(pos1 <= pos2);
510 auto ub1_iter = blocks_.upper_bound(
Block(pos1));
511 auto blk1_iter = --ub1_iter;
512 auto ub2_iter = blocks_.upper_bound(
Block(pos2));
513 auto blk2_iter = --ub2_iter;
514 size_t c = (blk1_iter->occ ? (blk1_iter->offset + blk1_iter->length - pos1 - 1) : 0);
515 c += (blk2_iter->occ ? (pos2 - blk2_iter->offset) : 0);
516 ++blk1_iter; --blk2_iter;
517 for(
auto i=blk1_iter; i!=blk2_iter; ++i)
524 boost::dynamic_bitset<> to_bitset()
const {
525 boost::dynamic_bitset<> result(nstates_);
526 for(
auto v : blocks_) {
528 for(
size_t l=0, pos=v.offset;
529 l<v.length; ++l, ++pos) {
539 Blocks result_blocks;
540 auto w = result_blocks.begin();
541 auto u =
other.blocks().begin();
542 auto v = blocks().begin();
545 size_t red_end = red->offset+red->length;
546 size_t black_end = black->offset+black->length;
547 const size_t size2 =
size() * 2;
548 while (red_end + black_end != size2) {
549 if (red_end <= black_end) {
550 w = result_blocks.insert(w,
Block(red->offset, red->length, red->occ ^ black->occ));
552 red_end = red->offset+red->length;
553 if (red->offset == black_end) {
555 black_end = black->offset+black->length;
559 w = result_blocks.insert(w,
Block(red->offset, black_end - red->offset, red->occ ^ black->occ));
561 black_end = black->offset+black->length;
562 std::swap(red,black);
563 std::swap(red_end,black_end);
566 w = result_blocks.insert(w,
Block(red->offset, red->length, red->occ ^ black->occ));
571 size_t hash_value()
const {
572 boost::hash<Blocks> h;
577 MPQC_ASSERT(nstates_ ==
other.nstates_);
578 return blocks_ ==
other.blocks_;
587 const size_t other_offset = nstates_;
588 nstates_ +=
other.nstates_;
589 MPQC_ASSERT(not
other.blocks_.empty());
590 MPQC_ASSERT(not blocks_.empty());
591 auto back_iter = (--blocks_.end());
592 if (back_iter->occ ==
other.blocks_.begin()->occ) {
593 Block merged_block(back_iter->offset, back_iter->length+
other.blocks_.begin()->length, back_iter->occ);
594 blocks_.erase(back_iter);
595 auto current_iter = blocks_.insert(merged_block).first;
596 for(
auto i = ++
other.blocks_.begin();
597 i!=
other.blocks_.end();
599 Block offset_block(i->offset + other_offset, i->length, i->occ);
600 current_iter = blocks_.insert(current_iter, offset_block);
604 auto current_iter = back_iter;
605 for(
auto i =
other.blocks_.begin();
606 i!=
other.blocks_.end();
608 Block offset_block(i->offset + other_offset, i->length, i->occ);
609 current_iter = blocks_.insert(current_iter, offset_block);
618 Blocks& blocks() {
return blocks_; }
619 const Blocks& blocks()
const {
return blocks_; }
621 Blocks::const_iterator block_pop_front(Blocks::const_iterator block) {
622 Block popped_block(*block);
623 --popped_block.length;
624 ++popped_block.offset;
625 blocks_.erase(block);
626 const auto result_iter = blocks_.insert(popped_block).first;
629 Blocks::const_iterator block_pop_back(Blocks::const_iterator block) {
633 Blocks::const_iterator block_push_front(Blocks::const_iterator block) {
634 Block popped_block(*block);
635 ++popped_block.length;
636 --popped_block.offset;
637 blocks_.erase(block);
638 const auto result_iter = blocks_.insert(popped_block).first;
641 Blocks::const_iterator block_push_back(Blocks::const_iterator block) {
646 Blocks::const_iterator block_replace(Blocks::const_iterator block,
647 const Block& new_block) {
648 blocks_.erase(block);
649 auto result = blocks_.insert(new_block);
650 MPQC_ASSERT(result.second ==
true);
654 Blocks::const_iterator block_erase(Blocks::const_iterator block) {
655 const bool have_prev = (block != blocks_.begin());
656 Blocks::const_iterator next_block(block); ++next_block;
657 const bool have_next = (next_block != blocks_.end());
660 Blocks::const_iterator prev_block(block); --prev_block;
661 prev_block->length += block->length;
662 blocks_.erase(block);
664 prev_block->length += next_block->length;
665 blocks_.erase(next_block);
669 else if (have_next) {
670 Block new_block(*next_block);
671 const size_t length = block->length;
672 new_block.offset -= length;
673 new_block.length += length;
674 blocks_.erase(block);
675 blocks_.erase(next_block);
676 auto result = blocks_.insert(new_block).first;
681 blocks_.insert(Block(0,nstates_,
false));
682 return blocks_.begin();
685 return blocks_.begin();
688 Blocks::const_iterator block_split(Blocks::const_iterator block,
690 const size_t orig_length = block->length;
691 const size_t offset = block->offset;
692 block->length = (pos - block->offset);
694 Block split_block(pos, 1, !block->occ);
695 auto iter = blocks_.insert(block, split_block);
697 Block remainder_block(pos+1, offset+orig_length-pos-1, block->occ);
698 blocks_.insert(iter, remainder_block);
705 template <
class CharT,
class Traits>
706 std::basic_ostream<CharT, Traits>&
707 operator<<(std::basic_ostream<CharT, Traits>& os,
708 const FermionOccupationBlockString& x)
714 FermionOccupationBlockString
operator+(
const FermionOccupationBlockString& s1,
715 const FermionOccupationBlockString& s2) {
716 FermionOccupationBlockString result(s1);
725 template <
typename FString>
728 typedef std::vector<FString> container_type;
729 typedef typename container_type::iterator iterator;
730 typedef typename container_type::const_iterator const_iterator;
731 typedef FString value_type;
736 const_iterator insert(
const FString& s) {
737 strings_.push_back(s);
738 return --strings_.end();
740 const_iterator insert(FString&& s) {
741 strings_.push_back(s);
742 return --strings_.end();
745 const_iterator begin()
const {
746 return strings_.cbegin();
748 const_iterator end()
const {
749 return strings_.cend();
752 size_t size()
const {
753 return strings_.size();
757 container_type strings_;
760 template <
typename FString>
763 typedef std::unordered_set<FString> container_type;
764 typedef typename container_type::iterator iterator;
765 typedef typename container_type::const_iterator const_iterator;
766 typedef FString value_type;
771 const_iterator insert(
const FString& s) {
772 auto result = strings_.insert(s);
773 MPQC_ASSERT(result.second ==
true);
776 const_iterator insert(FString&& s) {
777 auto result = strings_.insert(s);
778 MPQC_ASSERT(result.second ==
true);
782 const_iterator begin()
const {
783 return strings_.cbegin();
785 const_iterator end()
const {
786 return strings_.cend();
789 size_t size()
const {
790 return strings_.size();
794 container_type strings_;
805 template <
size_t Nb,
typename FString>
807 template <
typename FString>
810 static const size_t Rank = 1;
811 typedef typename FString::state_index_type state_index_type;
820 size_t to()
const {
return to_[0]; }
826 const std::array<state_index_type,1>&
cre()
const {
return to_; }
832 size_t from()
const {
return from_[0]; }
838 const std::array<state_index_type,1>&
ann()
const {
return from_; }
847 if (os[from_[0]] ==
true && (os[to_[0]] ==
false || from_ == to_)) {
848 os.remove(from_[0]).add(to_[0]);
861 if (os[from_[0]] ==
true && (os[to_[0]] ==
false || from_ == to_)) {
862 os.remove(from_[0]).add(to_[0]);
864 if (to_[0] > from_[0]) {
865 const bool sign_changes = os.count(from_[0],to_[0]) & size_t(1);
868 const bool sign_changes = os.count(to_[0],from_[0]) & size_t(1);
883 std::pair<bool,FString >
operator()(
const FString& os)
const {
885 const bool sign_changes = this->
apply_sign(result);
886 return std::make_pair(sign_changes,result);
890 std::array<state_index_type, Rank> from_;
891 std::array<state_index_type, Rank> to_;
894 template <
size_t Nb,
typename FString>
895 class FermionBasicNCOper {
897 static const size_t Rank = Nb;
898 typedef typename FString::state_index_type state_index_type;
900 template <
typename Int> FermionBasicNCOper(
const std::array<Int,Rank>& t,
const std::array<Int,Rank>& f) : to_(t) , from_(f) {
907 const std::array<state_index_type,Rank>&
cre()
const {
return to_; }
913 const std::array<state_index_type,Rank>&
ann()
const {
return from_; }
922 for(
size_t r=0; r<
Rank; ++r)
923 os.remove(from_[r]).add(to_[r]);
933 bool sign_changes =
false;
934 for (
size_t r=0; r<
Rank; ++r) {
935 os.remove(from_[r]).add(to_[r]);
937 if (to_[r] > from_[r]) {
938 sign_changes ^= os.count(from_[r],to_[r]) & size_t(1);
940 sign_changes ^= os.count(to_[r],from_[r]) & size_t(1);
951 std::pair<bool,FString >
operator()(
const FString& os)
const {
953 const bool sign_changes = this->
apply_sign(result);
954 return std::make_pair(sign_changes,result);
958 std::array<state_index_type, Rank> from_;
959 std::array<state_index_type, Rank> to_;
962 template <
class CharT,
class Traits,
size_t Nb,
typename FString>
963 std::basic_ostream<CharT, Traits>&
964 operator<<(std::basic_ostream<CharT, Traits>& os,
965 const FermionBasicNCOper<Nb, FString>& o)
967 os << o.to() <<
"<-" << o.from();
977 template <
typename FString,
unsigned int R,
bool GenerateStrings>
980 typedef typename FString::state_index_type state_index_type;
981 static const unsigned int Rank = R;
984 rstr_(ref_string), str_(ref_string)
994 const FString& operator*()
const {
1007 operator bool()
const {
1012 if (rstr_ ==
other.rstr_) {
1013 if (more_ !=
other.more_)
1016 return str_ ==
other.str_;
1023 static StringReplacementListIterator begin(
const FString& ref_string) {
1024 return StringReplacementListIterator(ref_string);
1027 static StringReplacementListIterator end(
const FString& ref_string) {
1028 return StringReplacementListIterator(ref_string,
false);
1035 FermionBasicNCOper<Rank,FString> oper_;
1039 nparticles_ = rstr_.count();
1040 if (nparticles_ == rstr_.size()
1050 StringReplacementListIterator(
const FString& ref_string,
1051 bool more) : rstr_(ref_string), str_(ref_string), more_(more) {
1060 template <
typename FermionStringSet>
1063 typedef typename FermionStringSet::value_type string_type;
1071 MPQC_ASSERT(nparticles_ <= nstates_);
1074 void operator()(FermionStringSet& sset) {
1078 for(
size_t k=nstates_-nparticles_+1; k>=1; --k) {
1084 do_iter(s0, sset, nstates_-k, nparticles_-1);
1092 static void do_iter(
const string_type& s0,
1093 FermionStringSet& sset,
1094 size_t nstates,
size_t nparticles) {
1096 if (nparticles == 1ul) {
1097 string_type s1(nstates);
1099 string_type s2(s1 + s0);
1102 for (
size_t pos = 1; pos < nstates; ++pos) {
1103 s2.remove(pos-1).add(pos);
1109 for (
size_t k = nstates - nparticles + 1; k >= 1; --k) {
1112 string_type s2(s1 + s0);
1113 do_iter(s2, sset, nstates - k, nparticles - 1);
1127 template <
size_t Ns>
1128 class hash<
sc::FermionOccupationNBitString<Ns> > {
1132 return s.hash_value();
1140 class hash<
sc::FermionOccupationDBitString > {
1144 return s.hash_value();
1152 class hash<
sc::FermionOccupationBlockString> {
1156 return s.hash_value();
1162 #endif // end of header guard