function win_prob(nVacancies, p, listSize, nVoters, nSimulations)
   v = 0.5*(flipud(eye(length(p),1)) - sqrt(p))
   Q = eye(length(v)) - (2/(v'*v)[1,1]) .* (v*v')
   M = diagm(sqrt(p./nVoters)) * Q[:,1:(end-1)]

   # Generate random columns of vote shares;
   # each column represents a multinomial experiment.
   f = M * randn((length(p)-1, nSimulations))
   for (j = 1:nSimulations)
      f[:,j] += p
   end
   # Normalize vote shares by 1/(number of vacancies)
   f .*= nVacancies

   # C[i,L] = total number of seats won by candidate i
   # of candidate list L throughout all experiments
   C = zeros(max(listSize), length(listSize))

   # Count C[i,j]
   for (j = 1:nSimulations)
      # no. of seats won by each list by meeting Hare quota
      won = max(min(listSize, ifloor(f[:,j])), 0)

      # remaining (normalized) vote shares
      r = f[:,j] - won

      # final no. of seats won by a candidate list
      won += (r .>= sort(r)[nVacancies - sum(won)])

      # add 1 to the no. of seats won by each elected candidate
      for L=1:length(listSize)
         C[1:won[L], L] += 1
      end
   end

   # Probability of winning a seat for each candidate
   C = C./float(nSimulations)

   # Sort those probabilities
   result = [(0,0,0.0)]; pop(result)
   for j in 1:size(C,2)
      for i in 1:listSize[j]
         push(result, (j, i, C[i,j]))
      end
   end
   sort!((x,y)->x[3]>y[3], result)
end

# The main program
begin
t0 = time()
nVacancies = 9
# NOW News' rolling survey, Sep 2-6, NT West
p = [.04 .07 .05 .07 .01 .02 .16 .09 .03 .03 .08 .04 .01 .08 .1 .12]'
# HKUPOP's rolling survey, Aug 1 - Sep 7, NT West
#p = [4.6 6.8 4.6 7.9 0 2.3 15.7 9.5 2.5 2.7 7 4.8 1.6 8.4 9.8 11.8]'*0.01
listSize = [6 5 5 4 1 1 2 4 1 1 2 7 7 3 2 4]'  # no. of candidates in each list
nVoters = 503
nSimulations = 1000000
result = win_prob(nVacancies, p, listSize, nVoters, nSimulations)
for i=1:length(result)
   println(result[i])
end
println(time()-t0)
end