Random Golf of the Day #3: Integer Partitions

Question

About the Series

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 2²¹). 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.

Leaderboard

The first post of the series generates a leaderboard.

To make sure that your answers show up, please start every answer with a headline, using the following Markdown template:

# 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.)

Answer 1

Python 2, 179 bytes

from random import*
m=r=input();i=q=r+1;h=[1]+[0]*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.

Answer 2

Dyalog APL, 67 59 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[1], length-1 from p[2], 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)

test n=5; test n=65; freq stats

Answer 3

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 [0] (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.

Answer 4

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;
    }

}

Answer 5

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

Answer 6

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 [0] if negative, or [1] if zero. ";
}2j
);p                  " Discard the total and print. ";

Answer 7

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.

Answer 8

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.

Answer 9

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;
        }
        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;
}

Adapted from https://stackoverflow.com/a/2163753/4230423 and golfed

How this works: We can calculate how many partitions of an integer n there are in O(n²) time. As a side effect, this produces a table of size O(n²) 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 :)

Answer 10

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],[1]+P][random()*(p+q)<p]
print d[N*R+M][1]

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.

Answer 11

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[0] = 999;
}

int main()
{
    int result[100], 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[65][64]), 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 {
        pushad;

        // 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:
        add ebp, [esi + ebx];
        sub bh, dl;
        jae sum_loop;

        add esi, 4;
        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;
        add edi, 4;
        sub ch, al;
        ja n_loop;

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

        dec edx;
        mov dh, (MAX_M + 1) * MAX_M * 4 / 256;
        add esp, edx;
        popad;
        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[0] = 999; // have to "use" the array to prevent optimizing it out
      }