import java.util.Arrays;
import java.util.Comparator;
import java.util.Random;

public class SimulatedIcm {

	private Random rnd = new Random();

	public double[] getEquities(double[] payouts, double[] stacks,
			int maxIterations) {
		if (payouts.length > stacks.length)
			payouts = Arrays.copyOf(payouts, stacks.length);

		double[] equities = new double[stacks.length];
		double[] exp = getExponents(stacks);
		final double[] r = new double[stacks.length];

		Integer[] ids = new Integer[stacks.length];
		for (int i = 0; i < ids.length; i++)
			ids[i] = i;

		Comparator<Integer> comp = new Comparator<Integer>() {
			@Override
			public int compare(Integer i1, Integer i2) {
				if (r[i1] < r[i2])
					return 1;
				if (r[i1] > r[i2])
					return -1;
				return 0;
			}
		};

		for (int iteration = 0; iteration < maxIterations; iteration++) {
			for (int i = 0; i < r.length; i++)
				r[i] = Math.pow(rnd.nextDouble(), exp[i]);
			Arrays.sort(ids, comp);
			for (int i = 0; i < payouts.length; i++)
				equities[ids[i]] += payouts[i];
		}

		for (int i = 0; i < equities.length; i++)
			equities[i] /= (double) maxIterations;

		return equities;
	}

	private static double[] getExponents(double[] stacks) {
		double t = 0;
		for (double s : stacks)
			t += s;
		t /= (double) stacks.length;
		double[] exp = new double[stacks.length];
		for (int i = 0; i < stacks.length; i++)
			exp[i] = t / stacks[i];
		return exp;
	}

	public static void main(String[] args) {
		SimulatedIcm sicm = new SimulatedIcm();

		double[] payouts = new double[] { 50, 30, 20 };
		double[] stacks = new double[] { 1000, 2000, 3000, 4000, 5000, 6000,
				7000, 8000, 9000 };
		double[] eqs = sicm.getEquities(payouts, stacks, 1000000);
		System.out.println(Arrays.toString(eqs));
	}
}

from random import random

def trial(weights, payouts):
	scores = sorted((random() ** weight, i) for i, weight in enumerate(weights))
	results = [0] * len(payouts)
	for payout, score in zip(payouts, scores): results[score[1]] = payout
	return results

def sicm(stacks, payouts, trials):
	avg = sum(stacks) / float(len(stacks))
	payouts += [0] * (len(stacks) - len(payouts))
	payouts.sort()
	weights = [avg / s for s in stacks]
	return [sum(player) / float(trials) for player in zip(*(
		trial(weights, payouts) for i in xrange(trials)
	))]

# I guess I should add an equivalent test to be fair to the Java implementation...
if __name__ == '__main__':
	print sicm([1000 * n for n in range(1, 10)], [50, 30, 20], 1000000)

SICM <- function(stacks, prizes, sims=1000) {
  if(length(prizes) != length(stacks)) {stop("Number of prizes not equal to number of remaining stacks. You might need to add prizes of value 0.")}
  normstacks <- 1/(stacks / mean(stacks))
  res <- matrix(nrow=sims, ncol=length(stacks))
    for(i in 1:sims) {
      rnums <- runif(length(stacks))
      places <- rnums^normstacks
      res[i,] <- prizes[length(places)+1 - rank(places, ties.method="random")]
    }
  out <- data.frame("Start.Stack"=stacks, "ICM.EV"=colMeans(res))
  warning(paste("Errors are proportional to", 1/sqrt(sims)))
  return(out)
}

stacks <- c(rep(1500,10))
prizes <- c(50, 30, 20, rep(0,7))

# Ex 1
SICM(stacks, prizes)

# Ex 2
SICM(stacks, prizes, sims=10000) # changing the sims for more precise results

# Ex 3
SICM(stacks=c(5000, 3000, 2000), prizes=c(50, 30, 20))

Example output:

> SICM(stacks, prizes)
   Start.Stack ICM.EV
1         1500  10.76
2         1500   9.43
3         1500  10.38
4         1500  10.26
5         1500   9.28
6         1500   9.97
7         1500  10.88
8         1500   9.16
9         1500  10.28
10        1500   9.60
Warning message:
In SICM(stacks, prizes) : Errors are proportional to 0.0316227766016838

> SICM(stacks, prizes, sims=10000)
   Start.Stack ICM.EV
1         1500  9.838
2         1500  9.976
3         1500 10.060
4         1500 10.212
5         1500  9.542
6         1500  9.923
7         1500 10.652
8         1500  9.593
9         1500 10.113
10        1500 10.091
Warning message:
In SICM(stacks, prizes, sims = 10000) : Errors are proportional to 0.01


> SICM(c(5000, 3000, 2000), c(50, 30, 20), sims=1000)
  Start.Stack ICM.EV
1        5000  37.81
2        3000  32.85
3        2000  29.34
Warning message:
In SICM(c(5000, 3000, 2000), c(50, 30, 20), sims = 1000) :
  Errors are proportional to 0.0316227766016838

> x <- runif(1)
> q <- runif(1,0,2)
> time <- proc.time()
> for(i in 1:10000000){
+   x^(1/q)
+ }
> proc.time()-time
   user  system elapsed 
  23.05    0.03   23.14 
> 
> time <- proc.time()
> for(i in 1:10000000){
+   (1/q)*log(x)
+ }
> proc.time()-time
   user  system elapsed 
  15.45    0.00   15.52