Re: [Bug-apl] RNG

[Juergen Sauermann]

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 <jinxiaoy...@gmail.com> 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 <juergen.sauerm...@t-online.de> 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 <juergen.sauerm...@t-online.de> >>>  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 >>>> >>>> >>>> >>>> >>>> >> >