Random Golf of the Day #3: Integer Partitions 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 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.

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, &lt;s&gt;104&lt;/s&gt; &lt;s&gt;101&lt;/s&gt; 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=[1]+[0]*q*q;exec&quot;h[i]=h[i+~q]+h[i-i%q*q];i+=1;&quot;*r*q
while r:
 x=random()*sum(h[r*q:r*q-~m]);m=0
 while x&gt;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.