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 >>>> >>>> >>>> >>>> >>>> >> >
