#include #include #include #include #include #include #include #include #include /* A269134 Number of combinatory separations of normal multisets of weight n. A multiset is "normal" if it spans an initial interval of positive integers. For example {1,1,2,2,2,3} is a normal multiset. Shorthand: 112223. The "type" of a multiset of integers is the unique normal multiset that has the same sequence of multiplicities when its entries are taken in increasing order. For example, the type of 335556 is 112223. Let h be a multiset. Then a pair "h<=g" is a "combinatory separation" iff g is the "types" of a multiset partition of h in at least one way. There are 14 distinct combinatory separations of normal multisets of weight 3: 111<={111}, 111<={1,11}, 111<={1,1,1}; 112<={112}, 112<={1,11}, 112<={1,12}, 112<={1,1,1}; 122<={122}, 122<={1,11}, 122<={1,12}, 122<={1,1,1}; 123<={123}, 123<={1,12}, 123<={1,1,1}. Note that different multiset partitions can lead to the same combinatory separation, e.g. {1,23}, {2,13} and {3,12} all lead to 123<={1,12}. */ /* This program iterates all "g" of total weight n, and then finds and counts all "h" that form a combinatory separation. There are 2^(n-1) normal multisets of weight n. These are stored as n-bit integers. For example, 112223 is stored as 0b101001 where the "1" bits indicate the positions where the element increases. These are A034691 normal multiset partitions, "g". These are stored as a vector of (value,count) pairs, in value order, each pair packed into 4 bytes. Given "g", the set of "h" is stored as a bitset where "h" is in the mathematical set iff (h-2^(n-1)) is in the bitset. This is cached up to A269134Calculator::cacheN. The set is calculated recursively (see A269134Calculator::calculate) as follows: If "g" is empty, the result is {{}}. If "g" has one element m, the result set is {m}. Otherwise iterate each non-empty subset of g: Extract all "1"s from the chosen subset to form a prefix. Define g2 from g without the extracted "1"s. Normalize each element. Recursively find the result set for g2. Sub-result = {prefix union {e+1 : e in h2} : h2 in result(g2)} Result for g is the union of the sub-results calculated above. Worked example g={11,12}: Subset {11}. Prefix=11. g2={12}. Result(g2)={12}. Sub-result(g)={1123}. Subset {12}. Prefix=1. g2={11,1}. Result(g2)={111,112,122}. Sub-result(g)={1222,1223,1233}. Subset {11,12}. Prefix=111. g2={1}. Result(g2)={1}. Sub-result(g)={1112}. Result(g)={1112,1123,1222,1223,1233}. These are the h such that h<=g is a combinatory separation. */ struct NormalMultiset { uint32_t value; operator uint32_t() const { return value; } int bitLength() const { return 32 - std::countl_zero(value); } NormalMultiset withoutTopBit() const { return NormalMultiset{value - (uint32_t{1} << (bitLength()-1))}; } friend std::ostream& operator<<(std::ostream&, NormalMultiset); }; std::ostream& operator<<(std::ostream& ostream, NormalMultiset nm) { char x = 0; for (int k = nm.bitLength()-1; k >= 0; k--) { x += (nm.value >> k) & 1; char c = ((x <= 9) ? '0' : char{'A' - 10}) + x; ostream << c; } return ostream; } struct NormalMultisetPartitionElement { NormalMultisetPartitionElement(uint32_t normalMultiset, uint8_t count) : data((normalMultiset << 5) + count) {} NormalMultiset normalMultiset() const { return NormalMultiset{data >> 5}; } uint8_t count() const { return data & 31; } void add(uint8_t delta) { if (count() + delta > 31) throw std::out_of_range("Out of range"); data += delta; } operator uint32_t() const { return data; } auto operator==(const NormalMultisetPartitionElement& other) const { return data == other.data; } auto operator<=>(const NormalMultisetPartitionElement& other) const { return data <=> other.data; } friend std::ostream& operator<<(std::ostream&, NormalMultisetPartitionElement); private: uint32_t data; }; std::ostream& operator<<(std::ostream& ostream, NormalMultisetPartitionElement elem) { ostream << elem.normalMultiset(); for (uint8_t i = 1; i < elem.count(); i++) ostream << ',' << elem.normalMultiset(); return ostream; } class NormalMultisetPartition { public: static void generate(uint8_t n, std::invocable auto fn); void generateSubsets(std::invocable auto fn) const; const std::vector& data() const { return data_; } uint8_t n() const { uint8_t s = 0; for (auto elem : data_) s += elem.normalMultiset().bitLength() * elem.count(); return s; } void add(uint32_t normalMultiset, uint8_t count); friend std::ostream& operator<<(std::ostream&, NormalMultisetPartition const&); auto operator==(const NormalMultisetPartition& other) const { return data_ == other.data_; } auto operator<=>(const NormalMultisetPartition& other) const { return data_ <=> other.data_; } private: std::vector data_; void gen(uint8_t n, uint8_t kMin, uint32_t nmMin, std::regular_invocable auto fn); }; void NormalMultisetPartition::generate(uint8_t n, std::invocable auto fn) { NormalMultisetPartition nmp; if (n) { nmp.gen(n, 1, 1, fn); } else { fn(nmp); } } void NormalMultisetPartition::generateSubsets(std::invocable auto fn) const { std::vector used(data_.size(), 0); std::vector nm2(data_.size()); for (size_t i = 0; i < data_.size(); i++) nm2[i] = data_[i].normalMultiset().withoutTopBit(); while (true) { for (size_t i = 0; true; i++) { if (i == used.size()) return; if (used[i] < data_[i].count()) { used[i]++; break; } used[i] = 0; } NormalMultisetPartition nmp2; for (size_t i = 0; i < data_.size(); i++) { if (used[i] && nm2[i]) { nmp2.add(nm2[i], used[i]); } if (used[i] < data_[i].count()) { nmp2.add(data_[i].normalMultiset(), data_[i].count() - used[i]); } } fn(nmp2); } } void NormalMultisetPartition::add(uint32_t normalMultiset, uint8_t count) { for (auto it = data_.begin(); true; ++it) { if (it == data_.end() || normalMultiset < (*it).normalMultiset()) { data_.emplace(it, normalMultiset, count); break; } else if (normalMultiset == (*it).normalMultiset()) { (*it).add(count); break; } } } void NormalMultisetPartition::gen(uint8_t n, uint8_t kMin, uint32_t nmMin, std::regular_invocable auto fn) { for (uint8_t k = kMin; k <= n; k++) { for (uint32_t nm = (k==kMin) ? nmMin : uint32_t{1} << (k-1); nm != (uint32_t{1} << k); nm++) { for (uint8_t c = (n/k) - !!(n%k); c >= 1; c--) { uint8_t n2 = n - k*c; data_.emplace_back(nm, c); if (n2) { gen(n2, k, nm+1, fn); } else { fn(*this); } data_.pop_back(); } } } } std::ostream& operator<<(std::ostream& ostream, NormalMultisetPartition const& nmp) { if (nmp.data_.empty()) return ostream << "()"; char c = '('; for (auto x : nmp.data_) { ostream << c << x; c = ','; } return ostream << ')'; } template <> struct std::hash { size_t operator()(const NormalMultisetPartition& nmp) const { size_t hash = 0; for (auto elem : nmp.data()) hash = hash * 31 + (uint32_t)elem; return hash; } }; //Reference to an object which may be temporary, stored within this object, or may exist elsewhere. template class Ref { public: Ref() : optional(std::in_place), ref(*optional) {} Ref(T& t) : ref(t) {} Ref(T&& t) : optional(t), ref(*optional) {} Ref(Ref& other) { if (other.optional.has_value()) { Ref(*(other.optional)); } else { Ref(*(other.ref)); } } Ref(Ref&& other) { if (other.optional.has_value()) { Ref(std::move(*(other.optional))); } else { Ref(*(other.ref)); } } constexpr const T* operator->() const noexcept { return &ref.get(); } constexpr T* operator->() noexcept { return &ref.get(); } constexpr const T& operator*() const& noexcept { return ref.get(); } constexpr T& operator*() & noexcept { return ref.get(); } private: std::optional optional; std::reference_wrapper ref; }; class Bitset { public: typedef uint64_t data_type; typedef size_t size_type; static constexpr unsigned int bits = std::numeric_limits::digits; Bitset() {} Bitset(data_type data) { if (data) this->data.emplace_back(data); } size_type capacityBits() const { return data.size() * bits; } bool empty() const { return data.empty(); } size_type count() const { size_type result = 0; for (data_type x : data) { result += std::popcount(x); } return result; } bool get(size_type pos) const; Bitset& set(size_type pos, bool value = true); void unionWithOffset(const Bitset& source, size_type bitOffset); public: //private std::vector data; void ensureSize(size_type reqdSize) { if (reqdSize > data.size()) data.resize(reqdSize); } }; bool Bitset::get(size_type pos) const { size_type idx = pos / bits; auto offset = pos % bits; return (idx < data.size()) && (data[idx] & (data_type{1} << offset)); } Bitset& Bitset::set(size_type pos, bool value) { size_type idx = pos / bits; auto shlOffset = pos % bits; ensureSize(idx + 1); if (value) { data[idx] |= data_type{1} << shlOffset; } else { data[idx] &= ~(data_type{1} << shlOffset); } return *this; }; void Bitset::unionWithOffset(const Bitset& source, size_type bitOffset) { if (source.empty()) return; size_type idxOffset = bitOffset / bits; auto shlOffset = bitOffset % bits; auto shrOffset = bits - shlOffset; size_type reqdSize = source.data.size() + idxOffset; if (shlOffset == 0) { ensureSize(reqdSize); for (size_type size = source.data.size(), i = 0, j = idxOffset; i < size; i++, j++) { data[j] |= source.data[i]; } } else { bool overflow = !!(source.data.back() >> shrOffset); reqdSize += overflow; ensureSize(reqdSize); data[idxOffset] |= source.data[0] << shlOffset; size_type size = source.data.size(), i = 1, j = idxOffset+1; for (; i < size; i++, j++) { data[j] |= (source.data[i-1] << shlOffset) | (source.data[i] >> shrOffset); } if (overflow) { data[j] |= source.data[i] >> shrOffset; } } } class A269134Calculator { public: const Ref bitset(const NormalMultisetPartition& nmp); private: static const uint8_t cacheN = 16; //2 GB at 16. Factor of 5 each additional +1. std::unordered_map cache; void calculate(const NormalMultisetPartition& nmp, Bitset& target, Bitset::size_type offset); }; const Ref A269134Calculator::bitset(const NormalMultisetPartition& nmp) { if (nmp.n() <= cacheN) { Bitset& result = cache[nmp]; if (result.empty()) { calculate(nmp, result, 0); } return Ref(result); } else { Bitset result; calculate(nmp, result, 0); return Ref(std::move(result)); } } void A269134Calculator::calculate(const NormalMultisetPartition& nmp, Bitset& target, Bitset::size_type offset) { if (nmp.data().empty()) { target.set(offset); } else if (nmp.data().size() == 1 && nmp.data()[0].count() == 1) { target.set(offset + nmp.data()[0].normalMultiset().withoutTopBit()); } else { nmp.generateSubsets([this,offset,&target](const NormalMultisetPartition& nmp2) { auto n2 = nmp2.n(); Bitset::size_type offset2 = offset + (n2 ? Bitset::size_type{1}<<(n2-1) : 0); if (n2 <= cacheN) { auto b2 = bitset(nmp2); target.unionWithOffset(*b2, offset2); } else { calculate(nmp2, target, offset2); } }); } } int main() { A269134Calculator a269134; for (uint8_t n = 0; n <= 27; n++) { uint_fast64_t a = 0; NormalMultisetPartition::generate(n, [&a,&a269134,n](const NormalMultisetPartition& nmp) { Ref bitset = a269134.bitset(nmp); auto p = bitset->count(); a += p; }); std::cout << +n << ' ' << a << std::format(" {:%H:%M:%OS}", std::chrono::system_clock::now()) << '\n' << std::flush; } }