From 0f996e934274eb3985ffc5a732120709be1c1857 Mon Sep 17 00:00:00 2001
From: Doug Coleman <doug.coleman@gmail.com>
Date: Sat, 29 Dec 2018 13:27:41 -0600
Subject: [PATCH] random.xoshiro: Implement xoshiro256starstar.

It's rather slow because Factor doesn't have 64bit integers.

[
	1,000,000 <iota>
	1 2 3 4 <xoshiro-256-star-star> '[ drop _ next-256 ] map
] profile

Reference implementation: http://xoshiro.di.unimi.it/xoshiro256starstar.c
Debugging C code: https://paste.factorcode.org/paste?id=3949

Related issue: #1518
---
 extra/random/xoshiro/authors.txt          |  1 +
 extra/random/xoshiro/tags.txt             |  1 +
 extra/random/xoshiro/xoshiro-tests.factor | 38 ++++++++++++
 extra/random/xoshiro/xoshiro.factor       | 70 +++++++++++++++++++++++
 4 files changed, 110 insertions(+)
 create mode 100644 extra/random/xoshiro/authors.txt
 create mode 100644 extra/random/xoshiro/tags.txt
 create mode 100644 extra/random/xoshiro/xoshiro-tests.factor
 create mode 100644 extra/random/xoshiro/xoshiro.factor

diff --git a/extra/random/xoshiro/authors.txt b/extra/random/xoshiro/authors.txt
new file mode 100644
index 0000000000..7c1b2f2279
--- /dev/null
+++ b/extra/random/xoshiro/authors.txt
@@ -0,0 +1 @@
+Doug Coleman
diff --git a/extra/random/xoshiro/tags.txt b/extra/random/xoshiro/tags.txt
new file mode 100644
index 0000000000..fb5bea3ff6
--- /dev/null
+++ b/extra/random/xoshiro/tags.txt
@@ -0,0 +1 @@
+rng
diff --git a/extra/random/xoshiro/xoshiro-tests.factor b/extra/random/xoshiro/xoshiro-tests.factor
new file mode 100644
index 0000000000..ef6a25667a
--- /dev/null
+++ b/extra/random/xoshiro/xoshiro-tests.factor
@@ -0,0 +1,38 @@
+! Copyright (C) 2018 Doug Coleman.
+! See http://factorcode.org/license.txt for BSD license.
+USING: kernel math random.xoshiro tools.test ;
+IN: random.xoshiro.tests
+
+{
+  2 3 131074 70368744177664
+  5760
+} [
+  0 1 2 3
+  (next-256)
+] unit-test
+
+{
+  4046638647718970624
+  2015137892536077249
+  6184416992699500823
+  16308606917844226410
+} [
+  0 1 2 3
+  10,000 [
+    (next-256) drop
+  ] times
+] unit-test
+
+{
+  14662298501051415801
+  12883398035623381500
+  17052052954271276209
+  1546841944388125985
+} [ 0 1 2 3 jump-256 ] unit-test
+
+{
+  15716266295256758020
+  17232205271518152816
+  9857397594961175947
+  8327361040835137714
+} [ 0 1 2 3 long-jump-256 ] unit-test
\ No newline at end of file
diff --git a/extra/random/xoshiro/xoshiro.factor b/extra/random/xoshiro/xoshiro.factor
new file mode 100644
index 0000000000..acfff6c26f
--- /dev/null
+++ b/extra/random/xoshiro/xoshiro.factor
@@ -0,0 +1,70 @@
+! Copyright (C) 2018 Doug Coleman.
+! See http://factorcode.org/license.txt for BSD license.
+USING: accessors alien.c-types classes.struct kernel locals math
+math.bitwise random sequences slots.syntax ;
+IN: random.xoshiro
+
+! http://xoshiro.di.unimi.it/xoshiro256starstar.c
+
+CONSTANT: JUMP-256 {
+    0x180ec6d33cfd0aba
+    0xd5a61266f0c9392c
+    0xa9582618e03fc9aa
+    0x39abdc4529b1661c
+}
+
+CONSTANT: LONG-JUMP-256 {
+    0x76e15d3efefdcbbf
+    0xc5004e441c522fb3
+    0x77710069854ee241
+    0x39109bb02acbe635
+}
+
+STRUCT: xoshiro-256-star-star { s0 ulonglong } { s1 ulonglong } { s2 ulonglong } { s3 ulonglong } ;
+
+: <xoshiro-256-star-star> ( s0 s1 s2 s3 -- obj )
+    xoshiro-256-star-star <struct>
+        swap >>s3
+        swap >>s2
+        swap >>s1
+        swap >>s0 ; inline
+
+: rotl-256 ( x: uint64_t k: int -- out: uint64_t )
+    [ shift ]
+    [ 64 swap - neg shift ] 2bi bitor 64 bits ; inline
+
+:: (next-256) ( s0! s1! s2! s3! -- s0 s1 s2 s3 64-random-bits )
+    s1 5 * 7 rotl-256 9 * 64 bits :> 64-random-bits
+    s1 17 shift 64 bits :> t
+    s0 s2 bitxor s2!
+    s1 s3 bitxor s3!
+    s2 s1 bitxor s1!
+    s3 s0 bitxor s0!
+    s2 t bitxor s2!
+    s3 45 rotl-256 s3!
+    s0 s1 s2 s3 64-random-bits ; inline
+
+: next-256 ( xoshiro-256-star-star -- r64 )
+    dup get[ s0 s1 s2 s3 ] (next-256)
+    [ set[ s0 s1 s2 s3 ] drop ] dip ; 
+
+:: jump ( s0! s1! s2! s3! jump-table -- s0' s1' s2' s3' )
+    0 0 0 0 :> ( t0! t1! t2! t3! )
+    4 <iota> [
+        64 <iota> [
+        [ jump-table nth ] [ 1 swap shift ] bi* bitand 0 > [
+            s0 t0 bitxor t0!
+            s1 t1 bitxor t1!
+            s2 t2 bitxor t2!
+            s3 t3 bitxor t3!
+        ] when
+        s0 s1 s2 s3 (next-256) drop s3! s2! s1! s0!
+        ] with each
+    ] each
+    t0 t1 t2 t3 ;
+
+: jump-256 ( s0 s1 s2 s3 -- s0' s1' s2' s3' ) JUMP-256 jump ;
+: long-jump-256 ( s0 s1 s2 s3 -- s0' s1' s2' s3' ) LONG-JUMP-256 jump ;
+
+M: xoshiro-256-star-star random-32*
+    next-256 ;