Source code in the tcl programming language

# A functional swap routine
proc swap {v idx1 idx2} {
    lset v $idx2 [lindex $v $idx1][lset v $idx1 [lindex $v $idx2];subst ""]
}

# Fill the integer array <vec> with the permutation at rank <rank>
proc computePermutation {vecName rank} {
    upvar 1 $vecName vec
    if {![info exist vec] || ![llength $vec]} return
    set N [llength $vec]
    for {set n 0} {$n < $N} {incr n} {lset vec $n $n}
    for {set n $N} {$n>=1} {incr n -1} {
	set r [expr {$rank % $n}]
	set rank [expr {$rank / $n}]
	set vec [swap $vec $r [expr {$n-1}]]
    }
}

# Return the rank of the current permutation.
proc computeRank {vec} {
    if {![llength $vec]} return
    set inv [lrepeat [llength $vec] 0]
    set i -1
    foreach v $vec {lset inv $v [incr i]}
    # First argument is lambda term
    set mrRank1 {{f n vec inv} {
	if {$n < 2} {return 0}
	set s [lindex $vec [set n1 [expr {$n - 1}]]]
	set vec [swap $vec $n1 [lindex $inv $n1]]
	set inv [swap $inv $s $n1]
	return [expr {$s + $n * [apply $f $f $n1 $vec $inv]}]
    }}
    return [apply $mrRank1 $mrRank1 [llength $vec] $vec $inv]
}


proc factorial {n} {
    for {set result 1; set i 2} {$i <= $n} {incr i} {
	set result [expr {$result * $i}]
    }
    return $result
}

set items1 [lrepeat 4 ""]
set rMax1 [factorial [llength $items1]]
for {set rank 0} {$rank < $rMax1} {incr rank} {
    computePermutation items1 $rank
    puts [format "%3d: \[%s\] = %d" \
	    $rank [join $items1 ", "] [computeRank $items1]]
}
puts ""
set items2 [lrepeat 21 ""]
set rMax2 [factorial [llength $items2]]
foreach _ {1 2 3 4} {
    # Note that we're casting to (potential) bignum, so entier() not int()
    set rank [expr {entier(rand() * $rMax2)}]
    computePermutation items2 $rank
    puts [format "%20lld: \[%s\] = %lld" \
	    $rank [join $items2 ", "] [computeRank $items2]]
}


proc uniform {limit} {
    set bits [expr {ceil(log($limit)/log(2))+10}]
    for {set b [set r 0]} {$b < $bits} {incr b 16} {
	incr r [expr {2**$b * int(rand() * 2**16)}]
    }
    return [expr {$r % $limit}]
}