Re: [Bug-apl] RNG

[Xiao-Yong Jin]

We need reproducibility, so /dev/urandom isn’t an option.
⎕FIO is a wonderful thing in gnu apl.
I still want something easy to manipulate in apl so I can play with different 
ideas easily.

> On May 18, 2016, at 11:56 AM, Juergen Sauermann <address@hidden> wrote:
> 
> Hi Elias,
> 
> yes, even better! I was thinking of */dev/random *but that one blocks 
> immediately
> due to lack of entropy.
> 
> /// Jürgen
> 
> 
> On 05/18/2016 06:41 PM, Elias Mårtenson wrote:
>> 
>> In that case, one can just open /dev/urandom and read random data from 
>> there. That algorithm is supposed to be cryptographically secure.
>> 
>> On 19 May 2016 12:39 a.m., "Juergen Sauermann" <address@hidden 
>> <mailto:address@hidden>> wrote:
>> 
>>    Hi,
>> 
>>    does not work because Xiao-Yong is looking for a portable solution.
>> 
>>    A simple and fast - although not portable - solution is this:
>> 
>>    1. write your favourite RNG in C/C++ (ie. copy the source code
>>    from Numerical Recipes
>>       and write the random number, say 8 bytes each, to stdout,
>> 
>>    2. *popen()* the C/C++ program with ⎕FIO like:
>> 
>>    *      Handle←⎕FIO[24] 'path-to-the**-C-program'*
>> 
>>    Every chunk of 8 bytes *fread()* with, say,
>> 
>>    *      256⊥8 ⎕FIO[6] Handle**
>>    *
>>    will then be one (signed) random number.
>>    You can also read several random numbers in one go).
>> 
>>    The above essentially pipes the output of your C/C++ RND directly
>>    into GNU APL.
>> 
>>    The code below will probably run faster if you can avoid to
>>    convert between integers and bit vectors
>>    too often (like in bit∆mult) and pre-compute constants like
>>    *(⌽¯1+⍳64)* beforehand.
>> 
>>    /// Jürgen
>> 
>> 
>>    On 05/18/2016 05:44 PM, Elias Mårtenson wrote:
>>>    How about implementing this as a native quad-function?
>>> 
>>>    On 18 May 2016 at 23:09, Xiao-Yong Jin <address@hidden
>>>    <mailto:address@hidden>> wrote:
>>> 
>>>        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
>>>        <mailto: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
>>>        <mailto: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
>>>        >>>>
>>>        >>>>
>>>        >>>>
>>>        >>>>
>>>        >>>>
>>>        >>
>>>        >
>>> 
>>> 
>>> 
>> 
>