Re: [Bug-apl] RNG

[Xiao-Yong Jin]

I translated a simple RNG from the book, Numerical Recipes.  APL code is slow, 
probably because of the bit operations.  The bit operations are not portable, 
relying on ¯1=2⊥64⍴1 (always a 64-bit signed integer), for which Dyalog does 
differently.  I’d welcome some suggestions.  Full code follows.

∇z←x bit∆add y
 ⍝ 64-bit-vector addition as unsigned int.
 z←(64⍴2)⊤2⊥x+y
∇

∇z←x bit∆mul y
 ⍝ 64-bit-vector multiplication as unsigned int.
 z←⊃bit∆add/(⌽¯1+⍳64)bit∆shift¨⊂[1]x∘.∧y
∇

∇z←n bit∆shift x
 ⍝ Shift vector x by n positions with filler ↑0↑x.
 ⍝ Shift direction is (1 ¯1=×n)/'left' 'right'.
 →(n≠0)/nonzero
 z←x
 →0
 nonzero: z←((×n)×⍴x)↑n↓x
∇

∇r←ran∆bits;u;v;w
 (u v w)←ran∆state
 u←((64⍴2)⊤7046029254386353087) bit∆add u bit∆mul ((64⍴2)⊤2862933555777941757)
 v←v≠¯17 bit∆shift v
 v←v≠31 bit∆shift v
 v←v≠¯8 bit∆shift v
 w←(¯32 bit∆shift w) bit∆add ((64⍴2)⊤4294957665) bit∆mul w∧¯64↑32⍴1
 ran∆state←u v w
 r←u≠21 bit∆shift u
 r←r≠¯35 bit∆shift r
 r←r≠4 bit∆shift r
 r←w≠r bit∆add v
∇

∇r←ran∆double
 ⍝ Returns a random 64-bit floating point number in the range of [0,1).
 r←5.42101086242752217E¯20×ran∆int64
 →(r≥0)/0
 r←r+1
∇

∇ran∆init j;u;v;w
 v←(64⍴2)⊤4101842887655102017
 w←(64⍴2)⊤1
 u←v≠(64⍴2)⊤j
 ran∆state←u v w
 ⊣ran∆int64
 ran∆state[2]←ran∆state[1]
 ⊣ran∆int64
 ran∆state[3]←ran∆state[2]
 ⊣ran∆int64
∇

∇r←ran∆int32
 r←2⊥¯32↑ran∆bits
∇

∇r←ran∆int64
 r←2⊥ran∆bits
∇

∇pass←ran∆test;n;r;x;y;z;ran∆state;expected;⎕CT
 ran∆init 12345
 n←0
 r←⍳0
 loop:x←ran∆int64
 y←ran∆double
 z←ran∆int32
 →(∧/n≠0 3 99)/nosave
 r←r,x,y,z
 nosave:→(10000>n←n+1)/loop
 nosave:→(100>n←n+1)/loop
 expected←¯2246894694364600497 0.9394142395716724714 1367803369 
2961111174787631927 0.1878554618793005226 3059533365 ¯1847334932704710330 
0.7547241536014889229 1532567919
 ⎕CT←1E¯15
 pass←(r=expected),0 1 2 9 10 11 297 298 299,r,[1.5]expected
∇

> On May 17, 2016, at 1:41 PM, Juergen Sauermann <address@hidden> wrote:
> 
> Hi Xiao-Yong,
> 
> I have fixed the redundant %, see SVN 728.
> 
> Regarding length, the GNU APL generator is 64-bit long (and so are
> GNU APL integers), which should suffice for most purposes.
> 
> Regarding bit vectors in APL, most people use integer 0/1 vectors
> for that (and then you have all boolean functions available) and
> 2⊥ resp. 2 2 2 ... 2⊤ for conversions between 0/1 vectors and integers
> You can also call into C/C++ libraries from GNU APL using native functions.
> 
> /// Jürgen
> 
> 
> 
> 
> On 05/17/2016 07:44 PM, Xiao-Yong Jin wrote:
>> Hi,
>> 
>> 
>>> On May 17, 2016, at 12:06 PM, Juergen Sauermann <address@hidden>
>>>  wrote:
>>> 
>>> Hi Xiao-Yong,
>>> 
>>> the reason is that ⎕RL is defined as a single integer in the ISO standard.
>>> That prevents generators with a too large state.
>>> 
>> Okay.  I was thinking whether ⎕RL can be an integer vector.  Even a combined 
>> generator with a 3-int-state would be quite useful for numerical simulations 
>> if it could be.
>> 
>>> For somebody seriously into simulations a general purpose generator
>>> will never suffice, but it is fairly easy to program one in APL.
>>> 
>> We definitely don’t want to make it cryptographically strong, but as a 
>> general purpose one, I would hope for decent high quality for relatively 
>> long monte carlo simulations.  
>> 
>> I don’t see an easy implementation because we don’t have exact 64bit 
>> unsigned integers and bit operations in APL.  If any of you on this list 
>> have suggestions in implementing a good RNG in APL, please let me know.
>> 
>>>  
>>> c++11 is currently not an option because I would like to not reduce the
>>> portability of GNU APL onto different platforms.
>>> 
>>> I'll have a look at the % usage.
>>> 
>>> /// Jürgen
>>> 
>>> 
>>> 
>>>  
>>> On 05/17/2016 06:16 PM, Xiao-Yong Jin wrote:
>>> 
>>>> Hi,
>>>> 
>>>> The LCG used for roll may be fine for some casual uses, but I would really 
>>>> like to see a higher quality RNG adopted here.  Since speed may not be the 
>>>> main concern here, employing the use of c++11 <random> and preferably 
>>>> using the std::mt19937_64 seems to be much better for any monte carlo type 
>>>> calculations.  It could be a trivial change to Quad_RL class, although 
>>>> saving the RNG state in a workspace may require a bit more code.  What do 
>>>> you say?
>>>> 
>>>> By the way, since in Workspace::get_RL 'return rand % mod;', you can 
>>>> remove the same ‘%’ in all the bif_roll definitions.
>>>> 
>>>> Best,
>>>> Xiao-Yong
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>> 
>