Implementation of random_permutation



Matthew Gwynne <mathew.gwynne at gmail.com> wrote:
> In our work, we make quite a lot of use of functions such as
> random_permutation for various experiments. In certain instances where
> performance is a big issue, we'd like to implement parts of the system
> in C++. However, we'd like to make sure that anything we write matches
> up with what we have at the Maxima level.
>
> To be able to do this, we need to know how the function
> random_permutation works to be able to replicate it in C++.
[...]
> Could anyone point me to some relevant documentation or maybe point me
> to some code for how these are implemented?

The code for random_permutation is in maxima/src/nset.lisp:

  (defun $random_permutation (a)
    (if ($listp a)
      (setq a (copy-list (cdr a)))
      (setq a (copy-list (require-set a "$random_permutation"))))
  
    (let ((n (length a)))
      (dotimes (i n)
        (let
          ((j (+ i ($random (- n i))))
           (tmp (nth i a)))
          (setf (nth i a) (nth j a))
          (setf (nth j a) tmp))))
  
    `((mlist) , at a))

Here's a translation into pseudo-code, with zero-based arrays:

  random_permutation(array a1)
  {
      a <-- copy_array(a1);
      n <-- length(a);
      for i from 0 to n-1 do
      {
          j <-- random integer (such that i <= j < n);
          swap(a[i],a[j]);
      }
      return(a);
  }

The random number generator used here is Mersenne Twister, MT19937, due
to Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
623-dimensionally equidistributed uniform pseudorandom number
generator.", ACM Transactions on Modeling and Computer Simulation, 1997.

      Mark