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