A problem on constrained combinatorics

Question

Not sure if this is a proper place, but I really don't know where else to ask. I'm craving for an algorithm generating certain sequences of numbers (the problem comes from physics).

I'm looking for all such sequences of $a_n$ an $b_n$ that: $$ \begin{cases} \sum \limits_{n=1}^{K} (n a_n) + \sum \limits_{m=1}^{K} (m b_m) = K\\ \sum \limits_{n=1}^{K} a_n - \sum \limits_{m=1}^{K} b_m = Q\\ \end{cases} \quad,\quad \text{where}\qquad \begin{cases} K\in\mathbb{Z}^+\\ Q\in\mathbb{Z}\\ a_n,b_m= \{0,1\} \end{cases} $$

For my application, I am interested in finding all such sequences starting from $K=Q(\operatorname{mod}2)$ and up to some large value of $K$, so, I guess, it would make sense to solve the problem iteratively. I would really appreciate if someone could up with an efficient solution.

EXAMPLE

For $Q=0$, $K=2$ the only option is $a_1=b_2=1$. We shall write it as $\{\{1\},\{1\}\}$. The first sublist corresponds to $a$'s, the second $-$ to $b$'s. The $j$th position in the sublist denotes the value of $a_j$ or $b_j$.

For $Q=0$, $K=4$, the possible sets are: $\{a_3=1;b_1=1\}$, $\{a_1=1;b_3=1\}$, $\{a_2=1;b_2=1\}$, with all other coefficients vanishing.

Here's an example of the output: \begin{alignat}{9} &Q=0,K=2:\quad &&\{\{\{1\},\{1\}\}\}\\ &Q=0,K=4:\quad &&\{\{\{0,0,1\},\{1\}\},\{\{1\},\{0,0,1\}\},\{\{0,1\},\{0,1\}\}\} \\&...&& \end{alignat}

Where I omitted the vanishing elements following the last non-zero element.

P.S.

The particular form of the output is not important for me. I would be equally happy to get it, say, in the form of lists of numbers of non-vanishing elements: \begin{alignat}{9} &Q=0,K=2:\quad &&\{\{\{1\},\{1\}\}\}\\ &Q=0,K=4:\quad &&\{\{\{3\},\{1\}\},\{\{1\},\{3\}\},\{\{2\},\{2\}\}\} \\&...&& \end{alignat}

UPDATE

1. By efficient algorithm I mean the one which does not make any useless steps (like "going though all possible sequences and throwing away the ones which do not work"). Yes, the algorithm will work in exponential time, and there's nothing we can do about it.

2. I think, a friend of mine came up with a solution.

At the first step, we generate the table having $K$ rows which contains all the partitions of $k=1..K$ into $m$ different terms:

$$ \begin{pmatrix} \{1\} & & \\ \{2\} & & \\ \{3\} & \{1,2\} &\\ \{4\} & \{1,3\} &\\ \{5\} & \{1,4\},\{2,3\} &\\ \{6\} & \{1,5\},\{2,4\},\{3,3\} & \{1,2,3\} \end{pmatrix} $$

Then, we consider all the possible ways to split $K=K_1+K_2$. For each pair $\{K_1,K_2\}$, we consider partitions of $K_1$ and $K_2$ having $n_1$ and $n_2$ elements, correspondingly. Of those, we are interested only in such that $n_1-n_2 = Q$.

I'm wondering if anyone could help with writing the pseudocode for this solution? I'm too stupid to deal with all the For's... I feel like it should look smth like this:

for (all partitions of K into K1 and K2)
    for (all lengths of partitions of K1)
        for (all partitions of K1 of given length)
            for (all partitions of K2 having the length (Q-current length of K1 partition) )
                take the current partitions!
            end
        end
    end
end

MOTIVATION

For the curious ones: the problem is motivated by representing the fermionic part of the system considered in this paper, see eq. (3.5).

In two dimensions, the state of the system containing fermions and neutral bosons is described by a list of numbers saying how many particles of each sort with given momentum we have: $$ |n_1,n_2,\ldots,n_N; \overline{n}_1,\overline{n}_2,\ldots,\overline{n}_{\overline{N}}; \widetilde{n}_1^{\widetilde{m}_1},\widetilde{n}_2^{\widetilde{m}_2},\ldots,\widetilde{n}_{\widetilde{N}}^{\widetilde{m}_{\widetilde{N}}}\rangle $$ Here $n_k$ is the number of fermions with momentum $k$, $\overline{n}_k$ is the number of anti-fermions with momentum $k$. Both can be either $0$ or $1$, due to the Pauli exclusion principle (or, better say, due to the spin-statistics theorem). As for bosons of each momentum, we can have an arbitrary number of those. The momenta are discretized due to the fact that the system is considered on a finite interval of space (much like in the case of the Fourier series). The momenta are positive because we are in the light-cone coordinates: if you are an observer is flying to the left with the speed of light, all the massive particles seem for you to be right-movers. We are interested in finding all the states with a given charge (the difference between the numbers of fermions and anti-fermions) and total momentum (the sum of momenta of all the particles in the state).

The first step is trivial, the total momentum can be separated into momenta of the bosonic and fermionic parts. Here I only asked about the latter, since the former is trivially given by all possible partitions of an integer.

Answer 1

In the following code I solve your problem. In the following code ba is the concatenation of binary vector b and $\alpha$ each one of them is of size K. So ba has length 2K. So vector a occupies the K least significant bits of ba and b occupies the most significant bits of b.

You can imagine the flow of the program as considering every possibility (0 or 1) for every bit starting from the least significant bit. If you are familiar with trees in cs I search the entire possible space in an preorder tree walk and I stop searching further the space under some certain conditions. These are:

1. In every level (n) down the tree I calculate maximum sum of remaining bits, if all of them are to be included, this maximum sum I call it rest. If sum+rest < K then there are not enough a or b occurrences left so the sum will ever be K so there is no point to continue further searching.
2. Since n_a - n_b = Q there are must be at least Q occurrences of $\alpha$ terms in the partition. In other words a vector must have at least Q bits set. So n_a is currently the bits set in $\alpha$ and K-n are the remaining maximum possible set bits. So if maximum remaining bits + n_a aren't enough to reach Q (that is if K-n + n_a
3. In a similar manner if the remaining occurrences of 'b' terms which are at most 2*K-n (when n > K) aren't enough to boost b_a such that n_a - n_b be exactly Q there is no point to traverse the tree downwards further.
4. And of course if at some point I hit the exactly sum ==K with n_a-n_b ==Q there is no point to continue further, as by doing so the remaining partitions might exceed the already achieved sum and Q.

>

#include <stdio.h>
#include <stdlib.h>


//Bad!! global variables
const int K = 50;
const int Q = 4;

void proc(int current_sum, int rest, int current_n_a, int current_n_b, int* ba, int n, int value);
void print_partition(int* ba);

int main(int argc, char const *argv[])
{
    int *ba;
    int rest;
    ba = calloc(2*K, sizeof(int));

    rest = 2 * (K/2 * (K+1)); //initialize to maximum possible sum if ba={1,1,1,1,1,1,1,1,...}
    proc(0, rest, 0, 0, ba, 1, 0);
    proc(0, rest, 0, 0, ba, 1, 1);
    return 0;
}


void proc(int current_sum, int rest, int current_n_a, int current_n_b, int* ba, int n, int value){
    int sum, n_a, n_b;
    ba[n-1] = value;
    n_a = current_n_a;
    n_b = current_n_b;
    if (n <= K){
        sum = current_sum + n*value;
        rest -= n;
        if (value){
            n_a = current_n_a + 1;
        }
    }else{
        sum = current_sum + (n-K)*value;
        rest -= (n-K);
        if (value){
            n_b = current_n_b + 1;
        }
    }

    if ((sum == K) && (n_a - n_b == Q)){
         print_partition(ba);   //we are done here
    }else{ 
        if ( (n < 2*K) && (sum < K) && (sum + rest >= K) ){ //backtracking conditions
            if ( (n <= K) && ( (K-n + n_a) >= Q ) /**/){ //ensure that there are enough occurences of a's left  so there are at least Q a's
                //recurse
                proc(sum, rest, n_a, n_b, ba, n+1, 0);
                proc(sum, rest, n_a, n_b, ba, n+1, 1);
            }
            if ( (n > K) && ((n_a - (2*K-n + n_b)) <= Q) /**/) { //ensure that there are enough possible occurences of b's left so n_b can be high enough in order n_a - n_b = Q 
                //recurse
                proc(sum, rest, n_a, n_b, ba, n+1, 0);
                proc(sum, rest, n_a, n_b, ba, n+1, 1);
            }
        }
    }
    ba[n-1] = 0; //reset value for the calling function     
}

void print_partition(int *ba){
    static int count = 1;
    int i;
    printf("#:%d    ", count);
    printf("{");
    printf("a:{");
    for (i = 0; i < K; i++){
        if (ba[i]){
            printf("%d,", i+1);
        }
    }
    printf("}b:{");
    for (i = 0; i < K; i++){
        if (ba[K+i]){
            printf("%d,", i+1);
        }
    }
    printf("}}\n");
    count++;

}

This code is what I would say a quick 'n' dirty solution. Sure there are ways to improve this code removing global variables, and improve the clarity of this code, anyone feel free to edit my answer for better code-clarity performance etc.

Answer 2

Rephrasing, $$A, B \subseteq \{1, 2, \ldots, K\} \\ \sum_{a \in A}a + \sum_{b \in B}b = K \\ |A| - |B| = Q$$

So as your friend suggests, this can be tackled by looking at partitions of an integer into a given number of distinct parts.

for K1 = 0 to K:
    for N1 = 0 to K1:
        for A in distinct_partitions(K1, N1):
            for B in distinct_partitions(K - K1, N1 - Q):
                process(A, B)

The easiest way to generate all partitions of an integer is a recursion with a second parameter: the largest part. Similarly, the easiest way to generate all partitions of an integer into a given number of distinct parts is a recursion with a third parameter: the largest part.

distinct_partitions(n, k):
    # Special case
    if n == 0 and k == 0:
        yield []

    for m = 1 to n:
        yield from distinct_partitions_inner(n, k, m)


distinct_partitions_inner(n, k, m):
    # Base case
    if k == 1:
        if n == m:
            yield [n]
        return

    # Parts are distinct, positive, and no greater than m
    if k > m:
        return

    # The smallest possible partition is 1+2+...+(k-1)+m
    if n < m + (k-1) * k / 2:
        return

    # The largest possible partition is m + (m-1) + ... + (m-k+1)
    if n > m * k - (k-1) * k / 2:
        return

    # m is the largest part, so we have n-m left to assign to smaller parts
    for i = 1 to m-1:
        for subpartition in distinct_partitions_inner(n-m, k-1, i):
            yield [m] + subpartition

The early aborts can be adjusted according to aesthetics to try to abort one level up the call stack.

---

Note that the table of the number of Fock states given in the Pauli and Brodsky paper you link contains a few errors. I derived their generating function and generated some data using the following Python code:

import math


def convolveP(a):
    ''' Convolves a generating function with the g.f. of the partition numbers '''
    for n in range(0, len(a)):
        for j in range(int((math.sqrt(1 + 24*n) + 1) / 6)):
            t = 1 + j * (3*j + 5) // 2
            a[n] += (-1)**j * a[n - t]
            if t + j < n:
                a[n] += (-1)**j * a[n - t - j - 1]
    return a


def fock(Q, n):
    triangles = [0] * n
    for j in range(int((math.sqrt(1 + 8*n) - 2*Q + 1) / 2)):
        triangles[(Q+j) * (Q+j+1) // 2] = (-1) ** j
    fermionic = convolveP(triangles)
    return convolveP(fermionic)


for Q in range(0, 5):
    print(fock(Q, 20))

This gave enough terms to find the sequences for the first four values of $Q$ in the Online Encyclopedia of Integer Sequences:

• $Q=0$: A006330
• $Q=1$: A001523 (correcting $a(0)$ to $0$: see comment by Michael Somos)
• $Q=2$: A288578, offset so that the first non-zero term is at $K=3$
• $Q=3$: A288579, offset so that the first non-zero term is at $K=6$

All four of those reference this paper, section 2.1, which gives some similar g.f.s (offset by a power of $q$) and a combinatorial explanation. I find that the explanation doesn't match up, but if I figure out the discrepancy and a suitable bijection to extract the partitions for fermions, antifermions, and bosons then I'll edit, because that might give a more economical generation procedure.