Arbitrary Randomness (Speed edition) Given integer n, calculate a set of n random unique integers in range 1..n^2 (inclusive) such that the sum of the set is equal to n^2

Random, in this case, means uniformly random between valid outputs. Each valid output for a given n must have a uniform chance of being generated.

For example, n=3 should have a 1/3 chance each of outputting 6, 1, 2, 3, 5, 1, or 4, 3, 2. As this is a set, order is irrelevant, 4, 3, 2 is identical to 3, 2, 4

Scoring

The winner is the program that can compute the highest n in under 60 seconds.
Note: To prevent possible partial hardcoding, all entries must be under 4000 bytes

Testing

All code will be run on my local Windows 10 machine (Razer Blade 15, 16GB RAM, Intel i7-8750H 6 cores, 4.1GHz, GTX 1060 in case you want to abuse the GPU) so please provide detailed instructions to run your code on my machine.
At request, entries can be run either via Debian on WSL, or on a Xubuntu Virtual Machine (both on the same machine as above)

Submissions will be run 50 times in succession, final score will be an average of all 50 results.
Rust, n ≈ 1400

How to run

Build with cargo build --release, and run with target/release/arbitrary-randomness n.

This program runs fastest with lots of memory (as long as it’s not swapping, of course).  You can adjust its memory usage by editing the MAX_BYTES constant, currently set at 8 GiB.

How it works

The set is constructed by a sequence of binary decisions (each number is either inside or outside the set), each of whose probabilities are computed combinatorially by counting the number of possible sets constructible after each choice using dynamic programming.

Memory usage for large n is reduced by using a version of this binomial partitioning strategy.

Cargo.toml

[package]
name = "arbitrary-randomness"
version = "0.1.0"
authors = ["Anders Kaseorg &lt;andersk@mit.edu&gt;"]

[dependencies]
rand = "0.6"


src/main.rs

extern crate rand;

use rand::prelude::*;
use std::env;
use std::f64;
use std::mem;

const MAX_BYTES: usize = 8 &lt;&lt; 30; // 8 gibibytes

fn ln_add_exp(a: f64, b: f64) -&gt; f64 {
    if a &gt; b {
        (b - a).exp().ln_1p() + a
    } else {
        (a - b).exp().ln_1p() + b
    }
}

fn split(steps: usize, memory: usize) -&gt; usize {
    if steps == 1 {
        return 0;
    }
    let mut u0 = 0;
    let mut n0 = f64::INFINITY;
    let mut u1 = steps;
    let mut n1 = -f64::INFINITY;
    while u1 - u0 &gt; 1 {
        let u = (u0 + u1) / 2;
        let k = (memory * steps) as f64 / u as f64;
        let n = (0..memory)
            .map(|i| (k - i as f64) / (i as f64 + 1.))
            .product();
        if n &gt; steps as f64 {
            u0 = u;
            n0 = n;
        } else {
            u1 = u;
            n1 = n;
        }
    }
    if n0 - (steps as f64) &lt;= steps as f64 - n1 {
        u0
    } else {
        u1
    }
}

fn gen(n: usize, rng: &amp;mut impl Rng) -&gt; Vec&lt;usize&gt; {
    let s = n * n.wrapping_sub(1) / 2;
    let width = n.min(MAX_BYTES / ((s + 1) * mem::size_of::&lt;f64&gt;()));
    let ix = |m: usize, k: usize| m + k * (s + 1);
    let mut ln_count = vec![-f64::INFINITY; ix(0, width)];
    let mut checkpoints = Vec::with_capacity(width);
    let mut a = Vec::with_capacity(n);
    let mut m = s;
    let mut x = 1;

    for k in (1..=n).rev() {
        let i = loop {
            let i = checkpoints.len();
            let k0 = *checkpoints.last().unwrap_or(&amp;0);
            if k0 == k {
                checkpoints.pop();
                break i - 1;
            }
            if i == 0 {
                ln_count[ix(0, i)] = 0.;
                for m in 1..=s {
                    ln_count[ix(m, i)] = -f64::INFINITY;
                }
            } else {
                for m in 0..=s {
                    ln_count[ix(m, i)] = ln_count[ix(m, i - 1)];
                }
            }
            let k1 = k - split(k - k0, width - 1 - i);
            for step in k0 + 1..=k1 {
                for m in step..=s {
                    ln_count[ix(m, i)] = ln_add_exp(ln_count[ix(m - step, i)], ln_count[ix(m, i)]);
                }
            }
            if k1 == k {
                break i;
            }
            checkpoints.push(k1);
        };

        while m &gt;= k &amp;&amp; rng.gen_bool((ln_count[ix(m - k, i)] - ln_count[ix(m, i)]).exp()) {
            m -= k;
            x += 1;
        }
        a.push(x);
        x += 1;
    }
    a
}

fn main() {
    if let [_, n] = &amp;env::args().collect::&lt;Vec&lt;_&gt;&gt;()[..] {
        let n = n.parse().unwrap();
        let mut rng = StdRng::from_entropy();
        println!("{:?}", gen(n, &amp;mut rng));
    } else {
        panic!("expected one argument");
    }
}


Try it online!

(Note: the TIO version has a few modifications.  First, the memory limit is reduced to 1 GiB. Second, since TIO doesn’t let you write a Cargo.toml and depend on external crates like rand, I instead pulled in drand48 from the C library using the FFI.  I didn’t bother to seed it, so the TIO version will produce the same result on every run.  Don’t use the TIO version for official benchmarking.)