Create integer partitions in multiplicity representation. Either all partitions of n or into exactly k parts. The basis is algorithm Z from A. Zoghbi: Algorithms for generating integer partitions, Ottawa (1993) http://dx.doi.org/10.20381/ruor-11312

partitions.cxx

// Copyright (c) 2017 Frederik Beaujean ([email protected])

// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#include <cassert>
#include <algorithm>
#include <iostream>
#include <vector>

using namespace std;

using Int_t = int;
using vec = std::vector<Int_t>;

/*
 * The output has to be filtered because Algorithm Z keeps higher
 * coefficients so the naive \sum_i c_i*y_i may exceed `n`.
 */
void filter_partition(vec& y, vec& c, const Int_t n)
{
    assert(y.size() == c.size());
    assert(y.size() > 0);

    Int_t sum = 0;
    for(size_t i=1; i < y.size(); ++i) {
        // check if n already reached, then zero the rest
        sum += c[i]*y[i];
        if (sum > n) {
            fill(y.begin() + i, y.end(), 0.0);
            fill(c.begin() + i, c.end(), 0.0);
            break;
        }
    }
}

/*
 * Print a partition of `n` in multiplicity representation.
 *
 * The output has to be filtered because Algorithm Z keeps higher
 * coefficients so the naive sum may exceed `n`.
 */
void print_partition(vec& y, vec& c, bool welcome=false)
{
    assert(y.size() == c.size());
    assert(y.size() > 0);

    if (welcome)
        cout << "sum_i c_i*y_i" << endl;

    bool init = false;
    for(size_t i=1; i < y.size(); ++i) {
        // only print valid entries
        if (c[i] > 0 && y[i] > 0) {
            std::cout << (init? " + " : "") << c[i] << "*" << y[i];
            init = true;
        }
    }
    std::cout << std::endl;
}

/*!
 * Code up the update step inside the while loop of algorithm Z from
 *
 * A. Zoghbi: Algorithms for generating integer partitions, Ottawa (1993)
 *
 * to generate the next integer partition of `n` in multiplicity
 * representation `\sum_i c_i * y_i` such that `c_i` is the
 * multiplicity of part `y_i`.
 *
 */
void algorithmZ_update(const Int_t n, vec& y, vec& c, Int_t& h)
{
    assert(y.size() == c.size());
    assert(y.size() > 0);

    // running index
    Int_t i = h-1;
    // calculated part
    Int_t k = c[h];
    // calculated remainder
    Int_t r = c[h] * y[h];
    // calculate the remainder
    while (y[h] - y[i] < 2) {
        k += c[i];
        r += c[i] * y[i];
        --i;
    }

    // update current part when it equals 1
    if (c[i] == 1) {
        if (i != 0) {
            r += c[i] * y[i];
            ++y[i];
        } else {
            i = 1;
            y[i] = 1;
        }
    }
    // update current part != 1
    else {
        --c[i];
        r += y[i];
        ++i;
        y[i] = y[i-1] + 1;
    }

    // calculate next parts based on remainder left from previous update
    c[i] = k;
    r -= c[i] * y[i];
    h = i+1;

    // update last modified part if it's the remainder
    if (r == y[i]) {
        ++c[i];
        h = i;
    }
    // add new part with multiplicity 1
    else {
        y[h] = r;
        c[h] = 1;
    }
    filter_partition(y, c, n);
}

/*!
 * Print all partitions of `n` in multiplicity representation where
 * `c_i` is the multiplicity and `y_i` is the part.
 *
 * Based on the pseudocode of Algorithm Z in A. Zoghbi: Algorithms for
 * generating integer partitions, Ottawa (1993), which generates all
 * partitions of `n`.
 *
 * Example: partition(5)
 sum_i c_i*y_i
 1*5
 1*1 + 1*4
 1*2 + 1*3
 2*1 + 1*3
 1*1 + 2*2
 3*1 + 1*2
 5*1
 */
void partition(const Int_t n)
{
    assert(n > 0);

    vec y(n+1, 0), c(n+1, 0);

    // initialize
    y[0] = -1;
    c.at(0) = 1;
    y.at(1) = n;
    c[1] = 1;
    // start with only one part
    Int_t h = 1;

    print_partition(y, c, true);

    // quit if all parts are 1
    while (c[1] != n) {
        algorithmZ_update(n, y, c, h);
        print_partition(y, c);
    }
}

/*!
 * Print all partitions of `n` into exactly `s` parts in multiplicity
 * representation where `c_i` is the multiplicity and `y_i` is the
 * part.
 *
 * Based on the pseudocode of Algorithm Z in A. Zoghbi: Algorithms for
 * generating integer partitions, Ottawa (1993), which generates all
 * partitions of `n`. I made the necessary modifications to partition
 * into `s` parts.
 *
 * Example: partition(8, 3)
 sum_i c_i*y_i
 2*1 + 1*6
 1*1 + 1*2 + 1*5
 1*1 + 1*3 + 1*4
 2*2 + 1*4
 1*2 + 2*3
*/
void partition(const Int_t n, const Int_t s)
{
    assert(n > 0);
    assert(s > 0);
    assert(n >= s);

    // have at most s distinct parts + one buffer value at index 0
    vec y(s+1, 0), c(s+1, 0);

    y[0] = -1;
    c[0] = 1;

    // h is the number of distinct parts
    Int_t h = 0;

    // the largest part of any partition of n into s parts
    const auto maxPart = n-s+1;

    // define initial partition

    // If there is only one distinct part then maxPart divides n.  But
    // part(6,3) = 5+1 = 4+2 = 3+3 shows that it doesn't have to be
    // the first partition, rather it is always the last in our order
    // because it cannot be mutated further. So if first and last
    // agree, there is only one partition. From a quick glance at the
    // partition table or the recurrence relation defining it T(n, 1)
    // = T(n, n) = 1, T(n, k) = 0 (k>n), T(n, k) = T(n-1, k-1) +
    // T(n-k, k), this reduces the choices to s=1,n-1,n. But for
    // s=n-1, the partition is n = (n-2)*1 + 1*2 so has two distinct
    // parts if n > 2. For n=2, s=n-1=1 is no special case, and for
    // n=1 all options collapse.
    if (s == 1 || s == n) {
        y[1] = maxPart;
        c[1] = n / maxPart;
        h = 1;
    } else {
        // initialize to the partition where the largest part has the
        // maximum value and everything else is 1. There are two
        // distinct parts
        y[1] = 1;
        c[1] = s-1;
        y[2] = maxPart;
        c[2] = 1;
        h = 2;
    }
    print_partition(y, c, true);

    // quit if all parts differ by at most one
    while (y[h] - y[1] > 1) {
        algorithmZ_update(n, y, c, h);
        print_partition(y, c);
    }
}

int main()
{
    partition(6);
    cout << "-------------\n";
    for (Int_t i=1; i <= 6; ++i)
        partition(6, i);

    return 0;
}

// Local Variables:
// compile-command:"g++ -std=c++11 -g -O2 partitions.cxx -o partitions && ./partitions"
// End: