Partitions in co-lexicographic order (PARI/GP algorithm without recursion)

Question

It is likely that many of you are familiar with partitions and the ways they are represented (I mean the lexicographic order and its variations).

Here are some partitions in co-lexicographic order:

[1]

[1, 1]
[2]

[1, 1, 1]
[2, 1]
[3]

[1, 1, 1, 1]
[2, 1, 1]
[3, 1]
[2, 2]
[4]

[1, 1, 1, 1, 1]
[2, 1, 1, 1]
[3, 1, 1]
[2, 2, 1]
[4, 1]
[3, 2]
[5]

Do you see the pattern? It looks like we only need to know

[1]
[2]
[3]
[2, 2]
[4]
[3, 2]
[5]

For a specific partition of n we just concatenate it with the vector of ones. Here the number of ones is equal to n minus sum over the all parts of a given partition.

This idea is not new (see A210941).

The problem is to generate a given row of A210941 without recursion.

I tried to write an efficient algorithm that works without recursion and this is what I got:

\\ Main function is list(vec).
\\ Input is a vector of the indices k.
\\ Output is a vector of the corresponding partitions.
\\ Here b4(k) is computed without recursion.
\\ So this program is good when you need to compute a few number of partitions.
\\ Especially when k is large.

b1(n) = {my(A = 1, B, v1); v1 = [1];
until(v1[A]>n, v1 = concat(v1, 0); B = 1;
while((binomial(B+1, 2) - binomial(B\2+1, 2)) <= A,
v1[A+1] += (-1)^((B-1)\2)*v1[A - binomial(B+1, 2) + binomial(B\2+1, 2) + 1]; B++); A++);
A-1}
list(vec) = {my(b2(n) = my(v1);
v1 = vector(n, i, vector(i\2+1, j, i==1));
for(i = 2, n, v1[i][i\2+1] = 1; v1[i][1] = vecsum(v1[i-1]);
for(j = 2, i\2, v1[i][j] = v1[i-1][j-1] - if((i-j) >= 2*(j-1), v1[i-j][j-1])));
for(i = 2, n, forstep(j = i\2, 1, -1, v1[i][j] += v1[i][j+1]));
v1,
vv1 = b2(b1(vecmax(vec))),
b3(n) = my(A = 0); until(vv1[A][1] > n, A++); A-1,
b4(n) = my(A = b3(n),
n = if(n == vv1[A][1], n, vv1[A][1] + vv1[A+1][1] - n),
B = n - vv1[A][1], C, v1);
v1 = []; while(!(B == 0), C = 1;
until(vv1[A+1][C]<=B, C++); v1 = concat(v1, C-1);
n -= vv1[A][1] - vv1[A - C + 1][1] + vv1[A+1][C];
A = b3(n); B = n - vv1[A][1]); v1 = Vecrev(v1);
v1 = concat(A + (#v1 > 0), v1));
vec = vector(#vec, i, b4(vec[i]))}

I need to evaluate the speed of this algorithm. To do this, it is desirable to compare it with something, so I am interested in alternative algorithms that work without recursion.

If you don't have any ideas for an alternative algorithm, please write in the comments how you rate the speed of my algorithm.

Answer 1

I don't have an alternative algorithm to benchmark against (I'll try to conjure up something after work tonight), but I can offer a way to analyze the efficiency of your algorithm. Its a bit too long for a comment, sorry.

What we can do is empirically estimate the time complexity of your solution given a decent spread of sample inputs. Start by running your algorithm for various input values n, and measure its runtime t, starting small and doubling the input until it begins to take too long to run (this can be seconds or minutes, depending on your patience).

Take your table of these n / t pairs, and plot them on a log-log plot. Then, find a best-fit line through your points (you may want to discard the smaller n values, since they'll be noisy, and skewed by lower order terms).

If your algorithm is approximately O(n^a) time complexity, then t ~= c * n^a + b for large n. Taking the log of both sides, log(t) = log(c * n^a + b), but the n^a term will wash out the constant term for large enough n, so we drop it: log(t) = log(c * n^a) = log(c) + a * log(n). You can see this is a linear equation in terms of log(t) and log(n), which is precisely what we've graphed. The slope of your best-fit line is (approximately) the observed degree of your time complexity, and the intercept is the log of the constant factor. This of course assumes you have large enough sample points that the observed times reflect the true time complexity, and that the time complexity is polynomial (excluding logarithmic factors).

If you believe your algorithm is exponential in nature, you can still use this trick but instead you have log(t) = log(c * a^n) = log(c) + log(a) * n, so only take the logarithm on the time domain, and then your base a is the exp of your slope.