Skip to content

Instantly share code, notes, and snippets.

@joriki
Last active July 4, 2018 07:39
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save joriki/6b409b2636027e22c98c604e06f32547 to your computer and use it in GitHub Desktop.
Save joriki/6b409b2636027e22c98c604e06f32547 to your computer and use it in GitHub Desktop.
Calculate and check the probability for at least one of m buckets to contain at least k of n uniformly randomly placed balls; see https://math.stackexchange.com/questions/2839811.
import java.math.BigInteger;
public class Question2839811 {
final static int n = 180;
final static int m = 10;
final static int k = 24;
// XOR shift pseudorandom number generator because Java's built-in generator has regularities that mess up the result
static long x = 4;
static long cache;
static int ncached;
static int nbits;
static {
for (int bits = m - 1;bits != 0;bits >>= 1)
nbits++;
}
static long randomLong () {
x ^= (x << 21);
x ^= (x >>> 35);
x ^= (x << 4);
return x;
}
static int randomBucket () {
for (;;) {
if (ncached == 0) {
cache = randomLong ();
ncached = 64 / nbits;
}
long bucket = (cache & ((1L << nbits) - 1));
cache >>= nbits;
ncached--;
if (bucket < m)
return (int) bucket;
}
}
final static int ntrials = 100000000;
static BigInteger [] powers = new BigInteger [m + 1];
public static void main (String [] args) {
for (int i = 0;i <= m;i++)
powers [i] = BigInteger.ONE;
for (int j = 1;j <= n;j++)
for (int i = 0;i <= m;i++)
powers [i] = powers [i].multiply (BigInteger.valueOf (i));
BigInteger factor = BigInteger.ONE;
BigInteger sum = BigInteger.ZERO;
for (int i = 0;i <= 3;i++) {
factor = factor.negate ();
if (i != 0) {
BigInteger pi = recurse (i);
sum = sum.add (pi.multiply (factor));
System.out.println (i + " : " + pi.doubleValue () / powers [m].doubleValue ());
System.out.println (" " + pi);
System.out.println (" " + sum.doubleValue () / powers [m].doubleValue ());
System.out.println (" " + sum);
}
factor = factor.multiply (BigInteger.valueOf (m - i)).divide (BigInteger.valueOf (i + 1));
}
long count = 0;
long [] counts = new long [4];
outer:
for (int j = 0;j < ntrials;j++) {
int [] buckets = new int [m];
for (int i = 0;i < n;i++)
buckets [randomBucket ()]++;
for (int b : buckets)
if (b >= k) {
count++;
continue outer;
}
// for (int i = 0;i < 4;i++)
// if (buckets [i] >= k)
// counts [i]++;
// else
// continue outer;
}
for (int i = 0;i < 4;i++)
System.out.println (i + " : " + counts [i] / (double) ntrials);
System.out.println (count / (double) ntrials);
}
static BigInteger recurse (int d) {
return recurse (0,d,BigInteger.valueOf (m - d),powers [m - d]);
}
static BigInteger recurse (int indexSum,int d,BigInteger md,BigInteger product) {
d--;
BigInteger sum = BigInteger.ZERO;
for (int i = 0;i <= n;i++,indexSum++) {
if (i >= k)
sum = sum.add (d == 0 ? product : recurse (indexSum,d,md,product));
product = product.multiply (BigInteger.valueOf (n - indexSum));
if (product.signum () == 0)
break;
product = product.divide (BigInteger.valueOf (i + 1));
product = product.divide (md);
}
return sum;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment