a simple rng in scheme

[Andy Wingo]

Hi,

Here is a simple implementation of what I think is xorshift128+, a
simple random number generator.  I have verified the first 10 outputs
against
https://github.com/AndreasMadsen/xorshift/blob/master/reference.c.

(define (xorshift128+ s0 s1)
  (define (u64 x) (logand #xffffFFFFffffFFFF x))
  (define (ulogxor x y) (u64 (logxor x y)))
  (define (ulsh x y) (u64 (ash x y)))
  (define (ursh x y) (u64 (ash x (- y))))
  (define (uadd x y) (u64 (+ x y)))
  (let ((state (u64vector s0 s1)))
    (lambda ()
      (let* ((x (u64vector-ref state 0))
             (y (u64vector-ref state 1))
             (x* y)
             (y* (ulogxor x (ulsh x 23)))
             (y* (ulogxor y* (ursh y* 17)))
             (y* (ulogxor y* (ulogxor y (ursh y 26)))))
        (u64vector-set! state 0 x*)
        (u64vector-set! state 1 y*)
        (let ((res (uadd y* y)))
          (values (logand res #xffffFFFF) (ash res -32)))))))

This function returns 2 32-bit values instead of one 64-bit value, to
avoid allocating bignums.  In fact with the latest Guile master, this
function does not allocate at all.  It can produce around 20M random
numbers per second on this old desktop.

The disassembly looks like:

   ----------------------------------------
   Disassembly of anonymous procedure at #x7fdb8320f168:

      0    (assert-nargs-ee/locals 1 5)    ;; 6 slots (0 args)   at 
/tmp/xorshift128+.scm:8:4
      1    (load-u64 4 0 0)                                      at 
/tmp/xorshift128+.scm:9:16
      4    (free-ref 5 5 0)                ;; free var 0
      6    (bv-u64-ref 3 5 4)              
      7    (load-u64 2 0 8)                                      at 
/tmp/xorshift128+.scm:10:16
     10    (bv-u64-ref 1 5 2)              
     11    (ulsh/immediate 0 3 23)                               at 
/tmp/xorshift128+.scm:4:26
     12    (ulogxor 3 3 0)                                       at 
/tmp/xorshift128+.scm:3:29
     13    (ursh/immediate 0 3 17)                               at 
/tmp/xorshift128+.scm:5:26
     14    (ulogxor 3 3 0)                                       at 
/tmp/xorshift128+.scm:3:29
     15    (ursh/immediate 0 1 26)                               at 
/tmp/xorshift128+.scm:5:26
     16    (ulogxor 0 1 0)                                       at 
/tmp/xorshift128+.scm:3:29
     17    (ulogxor 3 3 0)                 
     18    (bv-u64-set! 5 4 1)                                   at 
/tmp/xorshift128+.scm:15:8
     19    (bv-u64-set! 5 2 3)                                   at 
/tmp/xorshift128+.scm:16:8
     20    (uadd 5 3 1)                                          at 
/tmp/xorshift128+.scm:6:26
     21    (load-u64 4 0 4294967295)                             at 
/tmp/xorshift128+.scm:18:18
     24    (ulogand 4 5 4)                 
     25    (u64->scm 4 4)                  
     26    (ursh/immediate 5 5 32)                               at 
/tmp/xorshift128+.scm:18:42
     27    (u64->scm 3 5)                  
     28    (return-values 3)               ;; 2 values           at 
/tmp/xorshift128+.scm:18:10

As you can see, all unboxed operations.  To get this, I had to do a few
things to Guile:

 * Introduce unboxed logxor.

 * Compute the set of bits needed in any integer result, as a backwards
   flow problem.

 * Use that information to back-propagate constraints from
   (logand #xffffffffffffffff EXP) back to EXP, allowing EXP to unbox
   even if its output is out of the u64 range.

 * Elide (logand #xffffffffffffffff EXP) once EXP has been unboxed.

Regards,

Andy