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