# Random Golf of the Day #3: Integer Partitions

First off, you may treat this like any other code golf challenge, and answer it without worrying about the series at all. However, there is a leaderboard across all challenges. You can find the leaderboard along with some more information about the series in the first post.

Although I have a bunch of ideas lined up for the series, the future challenges are not set in stone yet. If you have any suggestions, please let me know on the relevant sandbox post.

# Hole 3: Integer Partitions

Time to increase the difficulty a bit.

A partition of a positive integer n is defined as a multiset of positive integers which sum to n. As an example if n = 5, the following partitions exist:

{1,1,1,1,1}
{2,1,1,1}
{2,2,1}
{3,1,1}
{3,2}
{4,1}
{5}


Note that these are multisets, so there is no order to them, {3,1,1}, {1,3,1} and {1,1,3} are all considered identical.

Your task is, given n, to generate a random partition of n. Here are the detailed rules:

• The distribution of partitions produced must be uniform. That is, in the above example, each partition should be returned with probability 1/7.

Of course, due to the technical limitations of PRNGs, perfect uniformity will be impossible. For the purpose of assessing uniformity of your submission, the following operations will be regarded as yielding perfectly uniform distributions:

• Obtaining a number from a PRNG (over any range), which is documented to be (approximately) uniform.
• Mapping a uniform distribution over a larger set of numbers onto a smaller set via modulo or multiplication (or some other operation which distributes values evenly). The larger set has to contain at least 1024 times as many possible values as the smaller set.
• Since the partitions are multisets you may return them in any order, and this order does not have to be consistent. However, for the purpose of the random distribution, the order is ignored. That is, in the above example, {3,1,1}, {1,3,1} and {1,1,3} together must have a probability of 1/7 of being returned.

• Your algorithm must have a deterministic runtime. In particular, you cannot generate random multisets and reject them if they don't sum to n.
• Your algorithm's time complexity must be polynomial in n. In particular, you cannot simply generate all partitions and select a random one (since the number of partitions grows exponentially with n). You may assume that the PRNG you're using can return uniformly distributed values in O(1) per value.
• You must not use any built-in function which solves this task.

You may write a full program or a function and take input via STDIN or closest alternative, command-line argument or function argument and produce output via return value or by printing to STDOUT (or closest alternative).

You may assume that n ≤ 65 (such that the number of partitions is less than 221). Output may be in any convenient, unambiguous list or string format.

If you submit a function, please consider also providing a little test program which calls the function a number of times and prints the results. It's okay if the parameters have to be tweaked in the code. This is just so people can check that the solution is at least approximately uniform.

This is code golf, so the shortest submission (in bytes) wins. And of course, the shortest submission per user will also enter into the overall leaderboard of the series.

The first post of the series generates a leaderboard.

# Language Name, N bytes


where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:

# Ruby, <s>104</s> <s>101</s> 96 bytes


(The language is not currently shown, but the snippet does require and parse it, and I may add a by-language leaderboard in the future.)

# Python 2, 179 bytes

from random import*
m=r=input();i=q=r+1;h=+*q*q;exec"h[i]=h[i+~q]+h[i-i%q*q];i+=1;"*r*q
while r:
x=random()*sum(h[r*q:r*q-~m]);m=0
while x>0:m+=1;x-=h[r*q+m]
print m;r-=m


I used formula (39) from this Knuth extract, which gives the number of partitions of n that have exactly m parts. This happens to equal the number of partitions of n that have m as the maximum element, which is the interpretation I'm using. The elements of the partition are generated from greatest to least. At each stage the formula is reused with the current remainder of n and maximum allowed element.

# Dyalog APL, 6759 51 bytes

p←{⍵,⊂1,⍨+/¨⌽⍵↑¨⍨⌽⍳⍴⍵}⍣⎕⊢⍬⋄f←{⍵=0:⍬⋄a,a∇⍵-a←{1++/(?+/⍵)>+\⍵}⍺↑⍵⊃p}⍨ (67 bytes)

p is a vector of vectors in which p[n][k] is the number of partitions of n into k summands, or equivalently: the number of partitions with greatest summand k. We build p by starting with the empty vector ⍬, reading n (the ⎕ reads input) and repeatedly applying the following:

{⍵,⊂1,⍨+/¨⌽⍵↑¨⍨⌽⍳⍴⍵}
⍴⍵   ⍝ the current length, initially 0
⍳⍴⍵   ⍝ 1 2 ... length
⌽⍳⍴⍵   ⍝ length ... 2 1
⍵↑¨⍨       ⍝ take length elements from p, length-1 from p, etc
⍝ padded with 0-s, e.g. if p was (,1)(1 1)(1 1 1)(1 2 1 1)(1 2 2 1 1):
⍝ we get:     (1 0 0 0 0)(1 1 0 0)(1 1 1)(1 2)(,1)
⌽           ⍝ reverse it: (,1)(1 2)(1 1 1)(1 1 0 0)(1 0 0 0 0)
+/¨            ⍝ sum each:   1 3 3 2 1
1,⍨               ⍝ append 1:   1 3 3 2 1 1
⍵,⊂                  ⍝ append the above to the vector of vectors


After n applications (⍣⎕), we have built p.

f picks a random partition. n f k is a random partition of at most k summands. f n is n f n.

{⍵=0:⍬⋄a,a∇⍵-a←{1++/(?+/⍵)>+\⍵}⍺↑⍵⊃p}⍨
⍨ ⍝ "selfie" -- use n as k if no k is provided
⍵=0:⍬                                 ⍝ if n=0 return empty
⍵⊃p   ⍝ pick the n-th element of p
⍺↑      ⍝ take k elements from that
{1++/(?+/⍵)>+\⍵}        ⍝ use them as weights to pick a random number 1...k
{           +\⍵}        ⍝   partial sums of weights
{    (?+/⍵)    }        ⍝   a random number 1...sum of weights
{    (?+/⍵)>+\⍵}        ⍝   which partial sums is it greater than?
{  +/          }        ⍝   count how many "greater than"-s
{1+            }        ⍝   we're off by one
a←                        ⍝ this will be the greatest number in our partition
a∇⍵-a                         ⍝ recur with n1=n-a and k1=a
a,                              ⍝ prepend a


Some improvements:

• inline p at the cost of slightly worse (but still good enough) performance

• in the computation of p rearrange ⌽ and 1, to save a character

• turn {1++/(?+/⍵)>+\⍵} into a train with 1+ in front: 1+(+/(?+/)>+\)

• make f an anonymous function and supply ⎕ (eval'ed input) as an argument to obtain a complete program

{⍵=0:⍬⋄a,a∇⍵-a←1+(+/(?+/)>+\)⍺↑⍵⊃{⍵,⊂⌽1,+/¨⍵↑¨⍨⌽⍳⍴⍵}⍣⍵⊢⍬}⍨⎕ (59 bytes)

Test with n=5

Test with n=65

And the following link runs n=5 thousands of times and gathers statistics about the frequency of each partition: ⎕rl←0 ⋄ {⍺,⍴⍵}⌸ {⍵=0:⍬⋄a,a∇⍵-a←1+(+/(?+/)>+\)⍺↑⍵⊃{⍵,⊂⌽1,+/¨⍵↑¨⍨⌽⍳⍴⍵}⍣⍵⊢⍬}⍨ ¨10000⍴5

More improvements, with help from Roger Hui:

• replace {⍵=0:A⋄B} with {×⍵:B⋄A}. Signum (×⍵) returns true for ⍵>0 and false for ⍵=0.

• replace (+/(?+/)>+\) with +/b<?⊃⌽b←+\, it saves a character

• use a matrix instead of vector of vectors to compute p: replace ⍵⊃{⍵,⊂⌽1,+/¨⍵↑¨⍨⌽⍳⍴⍵}⍣⍵⊢⍬ with ⊃↓(0,⍨⊢⍪⍨1 1⍉+\)⍣⍵⍪1.

{×⍵:a,a∇⍵-a←1++/b<?⊃⌽b←+\⍺↑⊃↓(0,⍨⊢⍪⍨1 1⍉+\)⍣⍵⍪1⋄⍬}⍨ (51 bytes)

• How does one get help from Roger Hui? Mar 1, 2015 at 13:52
• Write a toy APL interpreter to get yourself hired in the same company as him. Pose the above expression as a challenge, promise a pint of beer for every character he takes out. Then profit: fewer characters and more booze as he doesn't drink beer.
– ngn
Mar 1, 2015 at 20:57
• I see. That's a neat strategy, let's see if I can reproduce that… Can you ask him if Dyalog APL is going to get something like J's u/\. y any time soon? Mar 1, 2015 at 21:04
– ngn
Mar 2, 2015 at 19:04
• Thank you for asking him. Now I wonder if it's possible in linear time, too. Mar 2, 2015 at 19:26

## GolfScript, 90 bytes

~[[[1.]]]\({..[[{{(\{)}%+}%1$,1$,-=}%[1,]@0=+{1+}%]zip{{(\.,/*~}%.,.rand@=+}:^%]\+}*0=^(;


Online demo

This is an adaptation of my (simpler) partition counting code which instead of merely tracking a count tracks both a count and a uniformly selected one of the counted elements.

Side-by-side comparison of the two:

~[[[1.]]]\({..[[{{(\{)}%+}%1$,1$,-=}%[1,]@0=+{1+}%]zip{{(\.,/*~}%.,.rand@=+}:^%]\+}*0=^(;
[[ 1  ]]\({..[[{          1$,1$,-=}%  0 @0=+     ]zip{{+}*                }:^%]\+}*0=^


Differences:

• The initial ~ is because this is a program rather than a snippet.
• The [1.] replacing 1 corresponds to the change in what's tracked.
• The additional {(\{)}%+}% increments each element in that partition, and the {1+}% adds 1 to the partition.
• 0 becomes  (golfed to 1,) as part of the change in what's tracked, but because it needs to remain an array when prepended to another one it needs the extra [ ].
• The simple sum {+}* becomes a weighted selection from the partitions, combined with a summing of their count.
• The (; removes the count from the output and puts the partition into a nice format.

### Test framework

;7000,{;
'5'

~[[[1.]]]\({..[[{{(\{)}%+}%1$,1$,-=}%[1,]@0=+{1+}%]zip{{(\.,/*~}%.,.rand@=+}:^%]\+}*0=^(;

}%
:RESULTS
.&${ RESULTS.[2$]--,' '\n
}/


Tweak the initial 7000 if you want to run a different number of trials. Note that this is far too slow for an online demo.

# Pyth, 64 bytes

Uses https://stackoverflow.com/a/2163753/4230423 except that a)No cache since Pyth automatically memoizes, b)Prints each instead of appending to list, and c) is translated to Pyth.

M?smg-Gddr1hhS,GHG1Akd,QOgQQWQFNr1hhS,QkKg-QNNI<dKB-=dK)N=kN-=QN


I'll post an explanation of this when I have the time, but here is the corresponding python code:

g=lambda G,H: sum(map(lambda d:g(G-d, d), range(1, (H if H<G else G) + 1))) if G else 1
Q=input()
k,d = Q,random.randrange(g(Q, Q))
while Q:
for N in range(1, min(k, Q) + 1):
K = g(Q-N, N)
if d < K:
break
d -= K
print N
k=N
Q -= N


Edit: I finally got around to doing the explanation:

M                Lambda g(G,H)
?         G     If G truthy
s              Sum
m             Map
g            Recursive call
-Gdd        G-d,d
r            Range
1           1 to
h           +1
hS         First element of sorted (does min)
,GH       From G and H
1             Else 1
A                Double assign
kd              Vars k and d
,               To vals
Q              Q (evaled input)
O              Randrange 0 till val
gQQ           Call g(Q, Q)
WQ               While Q is truthy
FN              For N in
r              Range
1             From one
h             Till +1
hS,QK        Min(Q,K)
Kg             K=g(
-QN           Q-N
N             N
I<dK           If d<k
B             Break (implicit close paren)
-=dk           Subtracts d-=k
)               Close out for loop
N               Prints N
=kN             Set k=N
-=QN            Subtracts Q-=N


# Java, 285 267 bytes

int[][]p;void p(int n){p=new int[n+1][n+1];int a=n,b=k(n,a),c,d;for(b*=Math.random();n>0;System.out.print(c+" "),n-=a=c)for(c=0;c++<(a<n?a:n)&b>=(d=k(n-c,c));b-=d);}int k(int n,int k){if(p[n][k]<1)for(int a=0,b=0;b<k&b++<n;p[n][k]=a)a+=k(n-b,b);return n>0?p[n][k]:1;}


This is the same method as TheBestOne's answer, but it uses a simple array instead of a map. Also, instead of returning the random partition as a List, it prints them to the console.

Below is a test program that runs it 100000 times. For the example n=5, all sets were within 0.64% of a perfect 1/7 on my last run.

public class Partition {
public static void main(String[] args) {
Partition p = new Partition();
for(int i=0;i<100000;i++){
p.p(5);
System.out.println();
}
}

int[][]p;

void p(int n){
p=new int[n+1][n+1];
int a=n,b=k(n,a),c,d;
for(b*=Math.random();n>0;System.out.print(c+" "),n-=a=c)
for(c=0;c++<(a<n?a:n)&b>=(d=k(n-c,c));b-=d);
}

int k(int n,int k){
if(p[n][k]<1)
for(int a=0,b=0;b<k&b++<n;p[n][k]=a)
a+=k(n-b,b);
return n>0?p[n][k]:1;
}

}

• Although you've golfed the Math.min call down to (k<n?k:n), you can go further by ditching it entirely and just doing two checks: b<k&b++<n. You can also easily ditch the n>0 part of the loop conditional (since n>0&b<n reduces to b<n when b is guaranteed non-negative). Feb 24, 2015 at 12:10
• @PeterTaylor Thanks. Taking another look let me get rid of the extra return statement and the separate int declaration also. Feb 24, 2015 at 13:50

# CJam, 64 56 bytes

ri_L{_0>{\,f{)_@1$-j+}{)@)2$+:Umr@<@@?U+}*}{!a\;}?}2j);p


You can test it with this script:

ria100*{_L{_0>{\,f{)_@1$-j+}{)@)2$+:Umr@<@@?U+}*}{!a\;}?}2j);}%__|\f{_,\2$a-,-}2/p  ### Explanation ri_ " Read an integer and duplicate. "; L{ " Create a memoized function of the maximum and the sum, which returns a random partition, and the total number of partitions as the last item. "; _0> " If sum > 0: "; { \,f{ " For I in 0..max-1: "; )_@1$-   " Stack: I+1 I+1 sum-I-1 ";
j+       " Recursively call with the two parameters, and prepend I+1. ";
}
{            " Reduce on the results: ";
)@)2\$+   " Stack: partition1 total1 partition2 total1+total2 ";
:Umr     " U = total1+total2, then generate a random number smaller than that. ";
@<@@?    " If it is <total1, choose partition1, else choose partition2. ";
U+       " Append the total back to the array. ";
}*
}
{!a\;}?          " Else return  if negative, or  if zero. ";
}2j
);p                  " Discard the total and print. ";

• You should remove the incorrect "not golfed very well" part of your answer ;) Feb 26, 2015 at 20:09
• @anatolyg Removed. But I believe it is still possible to remove some bytes. I'm just too lazy to do that. Feb 26, 2015 at 22:52

# Octave, 200

function r=c(m)r=[];a=eye(m);a(:,1)=1;for(i=3:m)for(j=2:i-1)a(i,j)=a(i-1,j-1)+a(i-j,j);end;end;p=randi(sum(a(m,:)));while(m>0)b=a(m,:);c=cumsum(b);x=min(find(c>=p));r=[r x];p=p-c(x)+b(x);m=m-x;end;end


Ungolfed:

function r=c(m)
r=[];
a=eye(m);
a(:,1)=1;
for(i=3:m)
for(j=2:i-1)
a(i,j)=a(i-1,j-1)+a(i-j,j);
end;
end;
p=randi(sum(a(m,:)));
while(m>0)
b=a(m,:);
c=cumsum(b);
x=min(find(cumsum(b)>=p));
r=[r x];
p=p-c(x)+b(x);
m=m-x;
end
end


Construct a square matrix where each cell (m,n) reflects the number of partitions of m whose largest number is n, according to the Knuth extract @feersum so kindly cited. For example, 5,2 gives us 2 because there are two valid partitions 2,2,1 and 2,1,1,1. 6,3 gives us 3 for 3,1,1,1,3,2,1 and 3,3.

We can now deterministically find the p'th partition. Here, we are generating p as a random number but you can alter the script slightly so p is a parameter:

function r=c(m,p)
r=[];
a=eye(m);
a(:,1)=1;
for(i=3:m)
for(j=2:i-1)
a(i,j)=a(i-1,j-1)+a(i-j,j);
end;
end;
while(m>0)
b=a(m,1:m);
c=cumsum(b);
x=min(find(c>=p));
r=[r x];
p=p-c(x)+b(x);
m=m-x;
end
end


We can now deterministically show that each outcome is solely dependent on p:

octave:99> for(i=1:7)
> c(5,i)
> end
ans =

1   1   1   1   1

ans =

2   1   1   1

ans =

2   2   1

ans =

3   1   1

ans =

3   2

ans =

4   1

ans =  5


Thus, going back to the original where p is randomly generated, we can be assured each outcome is equally likely.

• I'm not sure about your 5,2 example. Shouldn't the two partitions be (2,2,1) and (2,1,1,1,1) (since the two you've listed have numbers greater than 2). Feb 24, 2015 at 19:38
• You're right, I got things twisted up. There are two partitions with two components, and two partitions starting with 2. I meant the latter. Feb 24, 2015 at 20:05

# R, 198 bytes

function(m){r=c();a=diag(m);a[,1]=1;for(i in 3:m)for(j in 2:(i-1))a[i,j]=a[i-1,j-1]+a[i-j,j];p=sample(sum(a[m,]),1);while(m>0){b=a[m,];c=cumsum(b);x=min(which(c>=p));r=c(r,x);p=p-c[x]+b[x];m=m-x};r}


Ungolfed:

f <- function(m) {
r <- c()
a <- diag(m)
a[, 1] <- 1
for (i in 3:m)
for (j in 2:(i-1))
a[i, j] <- a[i-1, j-1] + a[i-j, j]
p <- sample(sum(a[m, ]), 1)
while (m > 0) {
b <- a[m, ]
c <- cumsum(b)
x <- min(which(c >= p))
r <- c(r, x)
p <- p - c[x] + b[x]
m <- m - x
}
return(r)
}


It follows the same structure as @dcsohl's great solution in Octave, and is thus also based on the Knuth extract posted by @feersum.

I'll edit this later if I can come up with a more creative solution in R. In the meantime, any input is of course welcome.

# Java, 392 bytes

import java.util.*;Map a=new HashMap();List a(int b){List c=new ArrayList();int d=b,e=b(b,d),f=(int)(Math.random()*e),g,i;while(b>0){for(g=0;g++<Math.min(d, b);f-=i){i=b(b-g,g);if(f<i)break;}c.add(g);d=g;b-=g;}return c;}int b(int b,int c){if(b<1)return 1;List d=Arrays.asList(b,c);if(a.containsKey(d))return(int)a.get(d);int e,f;for(e=f=0;f++<Math.min(c, b);)e+=b(b-f,f);a.put(d,e);return e;}


Call with a(n). Returns a List of Integers

Indented:

import java.util.*;

Map a=new HashMap();

List a(int b){
List c=new ArrayList();
int d=b,e=b(b,d),f=(int)(Math.random()*e),g,i;
while(b>0){
for(g=0;g++<Math.min(d, b);f-=i){
i=b(b-g,g);
if(f<i)
break;
}
d=g;
b-=g;
}
return c;
}

int b(int b,int c){
if(b<1)
return 1;
List d=Arrays.asList(b,c);
if(a.containsKey(d))
return(int)a.get(d);
int e,f;
for(e=f=0;f++<Math.min(c, b);)
e+=b(b-f,f);
a.put(d,e);
return e;
}


How this works: We can calculate how many partitions of an integer n there are in O(n2) time. As a side effect, this produces a table of size O(n2) which we can then use to generate the kth partition of n, for any integer k, in O(n) time.

So let total = the number of partitions. Pick a random number k from 0 to total - 1. Generate the kth partition.

As usual, suggestions are welcome :)

# Python 2, 173 bytes

from random import*
N,M=input__
R=67;d=[(0,[])]*R*R
for k in range(R*R):p,P=d[k+~R];q,Q=d[k-k%R*R];d[k]=p+q+0**k,[[x+1 for x in Q],+P][random()*(p+q)<p]
print d[N*R+M]


Recursively makes a dictionary d, with keys k representing a pair (n,m) by k=67*n+m (using the guaranteed n<=65). The entry is the tuple of the number of partition of n into m parts, and a random such partition. The counts are computed by the recursive formula (thanks to feersum for pointing it out)

f(n,m) = f(n-1,m-1) + f(n,n-m),

and the random partition is updated by picking one of the two of its branches with probability proportional to its count. The update is done by adding is done by appending a 1 for the first branch and incrementing every element for the second.

I had a lot of trouble getting out-of-bounds values of m and n to give counts of zero. At first, I used a dictionary that defaults to a count of 0 and an empty list. Here, I'm using a list and padding it with this default entry instead. Negative indices cause the list to be read from its end, which gives a default entry is nothing near the end as ever reached, and wraparounds only touch a region where m>n.

# 80386 machine code, 105 bytes

Hexdump of the code:

60 8b fa 81 ec 00 41 00 00 33 c0 8b f4 33 d2 42
89 14 06 42 33 ed 8b d8 03 2c 1e 2a fa 73 f9 83
c6 04 89 2c 06 42 3b d1 76 ea fe c4 3a e1 76 db
33 d2 0f c7 f0 f7 f5 86 e9 85 d2 74 1b 33 c0 8d
34 0c 39 14 86 77 03 40 eb f8 2b 54 86 fc 40 89
07 83 c7 04 2a e8 77 e1 42 89 17 83 c7 04 fe cd
7f f7 4a b6 41 03 e2 61 c3


As a C function: void random_partition(int n, int result[]);. It returns the result as a list of numbers in the supplied buffer; it doesn't mark the end of the list in any way, but the user can discover the end by accumulating the numbers - the list ends when the sum is equal to n.

How to use (in Visual Studio):

#include <stdio.h>

__declspec(naked) void __fastcall random_partiton(int n, int result[])
{
#define a(byte) __asm _emit 0x ## byte
a(60) a(8b) a(fa) a(81) a(ec) a(00) a(41) a(00) a(00) a(33) a(c0) a(8b) a(f4) a(33) a(d2) a(42)
a(89) a(14) a(06) a(42) a(33) a(ed) a(8b) a(d8) a(03) a(2c) a(1e) a(2a) a(fa) a(73) a(f9) a(83)
a(c6) a(04) a(89) a(2c) a(06) a(42) a(3b) a(d1) a(76) a(ea) a(fe) a(c4) a(3a) a(e1) a(76) a(db)
a(33) a(d2) a(0f) a(c7) a(f0) a(f7) a(f5) a(86) a(e9) a(85) a(d2) a(74) a(1b) a(33) a(c0) a(8d)
a(34) a(0c) a(39) a(14) a(86) a(77) a(03) a(40) a(eb) a(f8) a(2b) a(54) a(86) a(fc) a(40) a(89)
a(07) a(83) a(c7) a(04) a(2a) a(e8) a(77) a(e1) a(42) a(89) a(17) a(83) a(c7) a(04) a(fe) a(cd)
a(7f) a(f7) a(4a) a(b6) a(41) a(03) a(e2) a(61) a(c3)
}

void make_stack() // see explanations about stack below
{
volatile int temp[65 * 64];
temp = 999;
}

int main()
{
int result, j = 0, n = 64, counter = n;
make_stack(); // see explanations about stack below

random_partiton(n, result);

while (counter > 0)
{
printf("%d ", result[j]);
counter -= result[j];
++j;
}
putchar('\n');
}


Example output (with n = 64):

21 7 4 4 3 3 3 3 2 2 2 2 2 1 1 1 1 1 1

This requires much explanations...

Of course I used the algorithm that everyone else used as well; there wasn't any choice with the requirement about complexity. So I don't have to explain the algorithm too much. Anyway:

I denote by f(n, m) the number of partitionings of n elements using parts no greater than m. I store them in a 2-D array (declared in C as f), where the first index is n, and the second m-1. I decided that supporting n=65 was too much trouble, so abandoned it...

Here is C code that calculates this table:

#define MAX_M 64
int f[(MAX_M + 1) * MAX_M];
int* f2;
int c; // accumulates the numbers needed to calculate f(n, m)
int m;
int k; // f(k, m), for various values of k, are accumulated
int n1;

for (n1 = 0; n1 <= n; ++n1)
{
f2 = f;
f2[n1 * MAX_M] = 1;
for (m = 2; m <= n; ++m)
{
c = 0;
k = n1;
while (k >= 0)
{
c += f2[k * MAX_M];
k -= m;
}
++f2;
f2[n1 * MAX_M] = c;
}
}


This code has some obfuscated style, so it can be converted to assembly language easily. It calculates the elements up to f(n, n), which is the number of partitionings of n elements. When this code is finished, the temporary variable c contains the needed number, which can be used to select a random partitioning:

int index = rand() % c;


Later, this index is converted to the required format (list of numbers) using the generated table.

do {
if (index == 0)
break;

m = 0;
f2 = &f[n * MAX_M];
while (f2[m] <= index)
{
++m;
}

index -= f2[m-1];
++m;
*result++ = m;
n -= m;
} while (n > 0);

do {
*result++ = 1;
--n;
} while (n > 0);


This code is also optimized for conversion to assembly language. There is a small "bug": if the partitioning doesn't contain any 1 number at the end, the last loop encounters n = 0, and outputs an unneeded 1 element. It doesn't hurt, however, because the printing code tracks the sum of the number, and doesn't print this extraneous number.

When converted to inline assembly, this code looks like this:

__declspec(naked) void _fastcall random_partition_asm(int n, int result[])
{
_asm {

// ecx = n
// edx = m
// bh = k; ebx = k * MAX_M * sizeof(int)
// ah = n1; eax = n1 * MAX_M * sizeof(int)
// esp = f
// ebp = c
// esi = f2
// edi = result

mov edi, edx;
sub esp, (MAX_M + 1) * MAX_M * 4; // allocate space for table
xor eax, eax;
row_loop:
mov esi, esp;
xor edx, edx;
inc edx;
mov dword ptr [esi + eax], edx;
inc edx;

col_loop:
xor ebp, ebp;
mov ebx, eax;

sum_loop:
sub bh, dl;
jae sum_loop;

mov [esi + eax], ebp;
inc edx;
cmp edx, ecx;
jbe col_loop;

inc ah;
cmp ah, cl;
jbe row_loop;

// Done calculating the table

// ch = n; ecx = n * MAX_M * sizeof(int)
// eax = m
// ebx =
// edx = index
// esp = f
// esi = f2
// ebp = c
// edi = result

xor edx, edx;
rdrand eax; // generate a random number
div ebp; // generate a random index in the needed range
xchg ch, cl; // multiply by 256

n_loop:
test edx, edx;
jz out_trailing;
xor eax, eax;
lea esi, [esp + ecx];

m_loop:
cmp [esi + eax * 4], edx;
ja m_loop_done;
inc eax;
jmp m_loop;
m_loop_done:

sub edx, [esi + eax * 4 - 4];
inc eax;
mov [edi], eax;
sub ch, al;
ja n_loop;

out_trailing:
inc edx;
out_trailing_loop:
mov dword ptr [edi], edx;
dec ch;
jg out_trailing_loop;

dec edx;
mov dh, (MAX_M + 1) * MAX_M * 4 / 256;
ret;
}
}


Some fun things to note:

• Generating a random number takes just 3 bytes of machine code (rdrand instruction)
• By a coincidence, the size of the table is 64, so the size of one row is 256 bytes. I use this to hold row indices in "high-byte" registers like ah, which gives me automatic multiplication by 256. To take advantage of this, I sacrificed the support for n = 65. I hope I can be excused for this sin...
• Allocating space on stack is performed by subtracting 0x4100 from the stack pointer register esp. This is a 6-byte instruction! When adding this number back, I managed to do it in 5 bytes:

    dec edx; // here edx = 1 from earlier calculations
mov dh, (MAX_M + 1) * MAX_M * 4 / 256; // now edx = 0x4100
add esp, edx; // this deallocates space on stack

• When debugging this function in MS Visual Studio, I found out that it crashes when it writes data to the space that it allocated on stack! After some digging around, I discovered some sort of stack overrun protection: OS seems to allocate only a very limited range of virtual addresses for stack; if a function accesses an address too far away, OS assumes it's an overrun and kills the program. However, if a function has many local variables, OS does some extra "magic" to make it work. So I have to call an empty function that has a large array allocated on stack. After this function returns, extra stack VM pages are allocated and can be used.

    void make_stack()
{
volatile int temp[65 * 64];
temp = 999; // have to "use" the array to prevent optimizing it out
}