#include <float.h>

#include <stdint.h>

#include <stdlib.h>

// TODO: use large period prng

static uint64_t seed = -1;

static uint32_t rand32(void)

{

	seed = 6364136223846793005ULL*seed + 1;

	return seed >> 32;

}

static uint64_t rand64(void)

{

	uint64_t u = rand32();

	return u<<32 | rand32();

}

static double frand()

{

	return rand64() * 0x1p-64;

}

static float frandf()

{

	return rand32() * 0x1p-32f;

}

static long double frandl()

{

	return rand64() * 0x1p-64L

#if LDBL_MANT_DIG > 64

+ rand64() * 0x1p-128L

#endif

;

}

void t_randseed(uint64_t s)

{

	seed = s;

}

/* uniform random in [0,n), n > 0 must hold */

uint64_t t_randn(uint64_t n)

{

	uint64_t r, m;

	/* m is the largest multiple of n */

	m = -1;

	m -= m%n;

	while ((r = rand64()) >= m);

	return r%n;

}

/* uniform on [a,b], a <= b must hold */

uint64_t t_randint(uint64_t a, uint64_t b)

{

	uint64_t n = b - a + 1;

	if (n)

		return a + t_randn(n);

	return rand64();

}

/* shuffle the elements of p and q until the elements in p are well shuffled */

static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)

{

	size_t r;

	uint64_t t;

	while (np) {

		r = t_randn(nq+np--);

		t = p[np];

		if (r < nq) {

			p[np] = q[r];

			q[r] = t;

		} else {

			p[np] = p[r-nq];

			p[r-nq] = t;

		}

	}

}

/* shuffle the elements of p */

void t_shuffle(uint64_t *p, size_t n)

{

	shuffle2(p,0,n,0);

}

void t_randrange(uint64_t *p, size_t n)

{

	size_t i;

	for (i = 0; i < n; i++)

		p[i] = i;

	t_shuffle(p, n);

}

/* hash table insert, 0 means empty, v > 0 must hold, len is power-of-2 */

static int insert(uint64_t *tab, size_t len, uint64_t v)

{

	size_t i = v & (len-1);

	size_t j = 1;

	while (tab[i]) {

		if (tab[i] == v)

			return -1;

		i += j++;

		i &= len-1;

	}

	tab[i] = v;

	return 0;

}

/* choose k unique numbers from [0,n), k <= n */

int t_choose(uint64_t n, size_t k, uint64_t *p)

{

	uint64_t *tab;

	size_t i, j, len;

	if (n < k)

		return -1;

	if (n < 16) {

		/* no alloc */

		while (k)

			if (t_randn(n--) < k)

				p[--k] = n;

		return 0;

	}

	if (k < 8) {

		/* no alloc, n > 15 > 2*k */

		for (i = 0; i < k;) {

			p[i] = t_randn(n);

			for (j = 0; p[j] != p[i]; j++);

			if (j == i)

				i++;

		}

		return 0;

	}

	// TODO: if k < n/k use k*log(k) solution without alloc

	if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {

		/* allocation is n-k < 4*k */

		tab = malloc((n-k) * sizeof *tab);

		if (!tab)

			return -1;

		for (i = 0; i < k; i++)

			p[i] = i;

		for (; i < n; i++)

			tab[i-k] = i;

		if (k < n-k)

			shuffle2(p, tab, k, n-k);

		else

			shuffle2(tab, p, n-k, k);

		free(tab);

		return 0;

	}

	/* allocation is 2*k <= len < 4*k */

	for (len = 16; len < 2*k; len *= 2);

	tab = calloc(len, sizeof *tab);

	if (!tab)

		return -1;

	for (i = 0; i < k; i++)

		while (insert(tab, len, t_randn(n)+1));

	for (i = 0; i < len; i++)

		if (tab[i])

			*p++ = tab[i]-1;

	free(tab);

	return 0;

}