Random selection form a list with power law probability (closed form solution) .

SelectInverseRank.js

//
// Select the n-th item with probability `1/n`.
//  Return an 0 based index (`i = n-1`)
//
// If parameter `a` is given, we select from probability `1/(n^a)`
// In particular, if a negative value of `a` is given, we select with probability `n^y, y = -a`;
//
// The code use the inverse of the cumulative probability function, to map the uniform distribution to the power distribution.
// 
// Caveat: 
//    The math is done with continous probability function and rounded at the last minute.
//    This can be seen as an estimate of doing discreete probability math from the start.
//

function SelectInverseRank(len, a = 1.0){

  if( abs(a - 1.0) < 1e-9 ){
      // Special case for a == 1.0
      var rnd = Math.random()
      return Math.floor( Math.pow( len + 1.0, rnd ) ) - 1.0  // Alternatively: Math.floor( Math.exp( Math.ln(len + 1.0) * rnd) ) - 1.0
  }
  else{
      var factor = Math.pow( len + 1.0 , 1.0 - a ) - 1.0
      var power = 1.0/(1.0 - a)
      var rnd = Math.random()
      return Math.floor( Math.pow( 1.0 +  factor * rnd, power ) ) - 1.0;
  }
  
}

// Input Ranked Items
// Output list of pairs, with random pairing of similar rank

list = ranked.copy()
n = len(list)

if(n%2==1){
j = SelectInverseRank(n, 1.0)
list.append( list[n-j-1] )
}

for (i=0:i<n:i++){
j = SelectInverseRank(n-i, 1.0)
tmp = item[i]
item[i] = item[j]
item[j] = tmp
}

pairs = []
for (i=0:i<n-2:i+=2){
pairs.append( [item[i],item[i+1]] )
}

# Select item n with probability of about 1/(n+a)
# With a > 0, used to decrease greediness
# We approximate the discrete probabilities with continuous ones.
#
# Let s(x,a) = integrate 1/(u+a) du for u = 0..x, with a>=0, x>=1
# We find that s(x,a) = ln((x+a)/(1+a))
# We solve for x in the equation s(x,a) == rand()*s(size,a) 
# We convert x to integer in [0,size[ and return it.

def SelectInverseRank(size, a = 0.0):
    a1 = a + 1.0
    return math.floor(a1*math.pow((a1+size)/a1, rand())-a1)

# But this is simpler and close when a=0
math.floor(rand()*rand()*size)
math.floor(abs(rand()-rand())*size)