1 #if !defined HAVE_PARTITION_DIST_ASC_SUBSET_LEX_H__ 2 #define HAVE_PARTITION_DIST_ASC_SUBSET_LEX_H__ 3 // This file is part of the FXT library. 4 // Copyright (C) 2012, 2013, 2014, 2015, 2019 Joerg Arndt 5 // License: GNU General Public License version 3 or later, 6 // see the file COPYING.txt in the main directory. 7 8 9 #include "comb/is-partition-asc.h" 10 11 #include "comb/comb-print.h" 12 #include "comb/is-partition-asc.h" 13 14 #include "fxttypes.h" 15 16 //#include <cmath> 17 18 19 class partition_dist_asc_subset_lex 20 // Integer partitions into distinct parts. 21 // Representation as list of parts in increasing order. 22 // Subset-lex order. 23 // The length of consecutive partitions changes by at most one. 24 // Only the last two positions in a partition at the end change. 25 // Loopless algorithm. 26 // See Joerg Arndt, Subset-lex: did we miss an order?, (2014) 27 // http://arxiv.org/abs/1405.6503 28 { 29 public: 30 ulong *a_; // partition: a[1] + a[2] + ... + a[m] = n 31 ulong n_; // integer partitions of n 32 ulong m_; // current partition has m parts 33 34 partition_dist_asc_subset_lex(const partition_dist_asc_subset_lex&) = delete; 35 partition_dist_asc_subset_lex & operator = (const partition_dist_asc_subset_lex&) = delete; 36 37 public: partition_dist_asc_subset_lex(ulong n)38 explicit partition_dist_asc_subset_lex(ulong n) 39 { 40 n_ = n; 41 42 // We need floor((sqrt(1+8*n)-1)/2) elements: 43 ulong k = 1, s = n_; 44 while ( s >= ( k + (k+1) ) ) { s -= k; k += 1; } 45 a_ = new ulong[k+1]; 46 a_[0] = 0; // sentinel: read (once) by test for right-extension 47 48 first(); 49 } 50 ~partition_dist_asc_subset_lex()51 ~partition_dist_asc_subset_lex() 52 { delete [] a_; } 53 data()54 const ulong * data() const { return a_ + 1; } 55 num_parts()56 ulong num_parts() const { return m_; } 57 first()58 void first() 59 { 60 if ( n_ == 0 ) { m_=0; return; } 61 a_[1] = n_; 62 m_ = 1; 63 } 64 65 next()66 ulong next() 67 // Return number of parts of generated partition. 68 // Return zero if the current is the last partition. 69 // Loopless algorithm. 70 { 71 ulong y = a_[m_-1]; // may read sentinel a[0] 72 ulong z = a_[m_]; 73 74 if ( z >= 2*y + 3 ) // can extend to the right 75 { // [*, Y, Z] --> [*, Y, Y+1, Z-1] 76 y += 1; 77 a_[m_] = y; 78 a_[m_+1] = z - y; // >= y 79 ++m_; 80 return m_; 81 } 82 else // add to the left 83 { 84 z -= 1; y += 1; 85 86 if ( z > y ) // add one unit to the left 87 { // [*, Y, Z] --> [*, Y+1, Z-1] 88 89 if ( m_<=1 ) return 0; // current is last 90 91 a_[m_-1] = y; 92 a_[m_] = z; 93 return m_; 94 } 95 else // add to one unit second left and combine last with second last 96 { // [*, X, Y, Z] --> [*, X+1, Y+Z] 97 98 if ( m_<=2 ) return 0; // current is last 99 100 a_[m_-2] += 1; 101 a_[m_-1] += z; 102 --m_; 103 return m_; 104 } 105 } 106 } 107 108 OK()109 bool OK() const 110 { 111 if ( ! is_partition_asc(data(), num_parts(), n_) ) return false; 112 113 for (ulong j=1; j<m_; ++j) // parts distinct? 114 if ( a_[j] == a_[j+1] ) 115 return false; 116 117 return true; 118 } 119 120 void print(const char *bla, bool dfz=false) const 121 { print_vec(bla, data(), m_, dfz); } 122 }; 123 // ------------------------- 124 125 126 #endif // !defined HAVE_PARTITION_DIST_ASC_SUBSET_LEX_H__ 127