Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Lukas Weber
load_leveller
Commits
0454c645
Commit
0454c645
authored
Feb 14, 2019
by
Lukas Weber
Browse files
clang formatted the code and in transition to fixing the build system
parent
4828f95b
Changes
30
Hide whitespace changes
Inline
Side-by-side
CMakeLists.txt
View file @
0454c645
cmake_minimum_required
(
VERSION
2.6
)
cmake_minimum_required
(
VERSION
3.1
)
project
(
load_leveller
)
include
(
${
CMAKE_BINARY_DIR
}
/conanbuildinfo.cmake
)
conan_basic_setup
()
find_package
(
MPI REQUIRED
)
set
(
CMAKE_CXX_COMPILER
${
MPI_CXX_COMPILER
}
)
set
(
CMAKE_CXX_FLAGS
"
${
CMAKE_CXX_FLAGS
}
-g -Wall -O3 -std=c++17 -pedantic
${
MPI_LINK_FLAGS
}
${
MPI_COMPILE_FLAGS
}
"
)
include_directories
(
${
MPI_INCLUDE_PATH
}
)
find_package
(
fmt REQUIRED
)
find_package
(
HDF5 REQUIRED
)
set
(
CMAKE_CXX_STANDARD 17
)
set
(
SRCs
dump.cpp
...
...
@@ -25,5 +26,5 @@ set(SRCs
)
add_library
(
load_leveller
${
SRCs
}
)
target_link_libraries
(
load_leveller
${
CONAN_LIBS
}
${
MPI_LIBRARIES
}
)
target_link_libraries
(
load_leveller
PUBLIC MPI::MPI_CXX fmt::fmt
)
install
(
FILES load_leveller.h DESTINATION include
)
MersenneTwister.h
View file @
0454c645
...
...
@@ -22,7 +22,7 @@
// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
// Copyright (C) 2000 - 2003, Richard J. Wagner
// All rights reserved.
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
...
...
@@ -35,8 +35,8 @@
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. The names of its contributors may not be used to endorse or promote
// products derived from this software without specific prior written
// 3. The names of its contributors may not be used to endorse or promote
// products derived from this software without specific prior written
// permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
...
...
@@ -73,139 +73,156 @@
#include
<iostream>
#include
<limits.h>
#include
<math.h>
#include
<stdio.h>
#include
<time.h>
#include
<math.h>
class
MTRand
{
// Data
// Data
public:
typedef
unsigned
long
uint32
;
// unsigned integer type, at least 32 bits
static
const
int
N
=
624
;
// length of state vector
typedef
unsigned
long
uint32
;
// unsigned integer type, at least 32 bits
static
const
int
N
=
624
;
// length of state vector
protected:
enum
{
M
=
397
};
// period parameter
uint32
state
[
N
];
// internal state
uint32
*
pNext
;
// next value to get from state
int
left
;
// number of values left before reload needed
enum
{
M
=
397
};
// period parameter
uint32
state
[
N
];
// internal state
uint32
*
pNext
;
// next value to get from state
int
left
;
// number of values left before reload needed
uint32
myoneseed
;
//Methods
//
Methods
public:
MTRand
(
const
uint32
&
oneSeed
);
// initialize with a simple uint32
MTRand
(
uint32
*
const
bigSeed
,
uint32
const
seedLength
=
N
);
// or an array
MTRand
();
// auto-initialize with /dev/urandom or time() and clock()
MTRand
(
const
uint32
&
oneSeed
);
// initialize with a simple uint32
MTRand
(
uint32
*
const
bigSeed
,
uint32
const
seedLength
=
N
);
// or an array
MTRand
();
// auto-initialize with /dev/urandom or time() and clock()
// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
// values together, otherwise the generator state can be learned after
// reading 624 consecutive values.
// Access to 32-bit random numbers
double
rand
();
// real number in [0,1]
double
rand
(
const
double
&
n
);
// real number in [0,n]
double
randExc
();
// real number in [0,1)
double
randExc
(
const
double
&
n
);
// real number in [0,n)
double
randDblExc
();
// real number in (0,1)
double
randDblExc
(
const
double
&
n
);
// real number in (0,n)
uint32
randInt
();
// integer in [0,2^32-1]
uint32
randInt
(
const
uint32
&
n
);
// integer in [0,n] for n < 2^32
double
operator
()()
{
return
rand
();
}
// same as rand()
double
rand
();
// real number in [0,1]
double
rand
(
const
double
&
n
);
// real number in [0,n]
double
randExc
();
// real number in [0,1)
double
randExc
(
const
double
&
n
);
// real number in [0,n)
double
randDblExc
();
// real number in (0,1)
double
randDblExc
(
const
double
&
n
);
// real number in (0,n)
uint32
randInt
();
// integer in [0,2^32-1]
uint32
randInt
(
const
uint32
&
n
);
// integer in [0,n] for n < 2^32
double
operator
()()
{
return
rand
();
}
// same as rand()
// Access to 53-bit random numbers (capacity of IEEE double precision)
double
rand53
();
// real number in [0,1)
double
rand53
();
// real number in [0,1)
// Access to nonuniform random number distributions
double
randNorm
(
const
double
&
mean
=
0.0
,
const
double
&
variance
=
0.0
);
double
randNorm
(
const
double
&
mean
=
0.0
,
const
double
&
variance
=
0.0
);
// Re-seeding functions with same behavior as initializers
void
seed
(
const
uint32
oneSeed
);
void
seed
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
=
N
);
void
seed
(
const
uint32
oneSeed
);
void
seed
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
=
N
);
void
seed
();
uint32
get_seed
()
{
return
myoneseed
;}
uint32
get_seed
()
{
return
myoneseed
;
}
// Saving and loading generator state
void
save
(
std
::
vector
<
uint32
>
&
saveArray
)
const
;
void
load
(
const
std
::
vector
<
uint32
>
&
loadArray
);
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
MTRand
&
mtrand
);
friend
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
MTRand
&
mtrand
);
void
save
(
std
::
vector
<
uint32
>
&
saveArray
)
const
;
void
load
(
const
std
::
vector
<
uint32
>
&
loadArray
);
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
MTRand
&
mtrand
);
friend
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
MTRand
&
mtrand
);
protected:
void
initialize
(
const
uint32
oneSeed
);
void
initialize
(
const
uint32
oneSeed
);
void
reload
();
uint32
hiBit
(
const
uint32
&
u
)
const
{
return
u
&
0x80000000UL
;
}
uint32
loBit
(
const
uint32
&
u
)
const
{
return
u
&
0x00000001UL
;
}
uint32
loBits
(
const
uint32
&
u
)
const
{
return
u
&
0x7fffffffUL
;
}
uint32
mixBits
(
const
uint32
&
u
,
const
uint32
&
v
)
const
{
return
hiBit
(
u
)
|
loBits
(
v
);
}
uint32
twist
(
const
uint32
&
m
,
const
uint32
&
s0
,
const
uint32
&
s1
)
const
{
return
m
^
(
mixBits
(
s0
,
s1
)
>>
1
)
^
(
-
loBit
(
s1
)
&
0x9908b0dfUL
);
}
static
uint32
hash
(
time_t
t
,
clock_t
c
);
uint32
hiBit
(
const
uint32
&
u
)
const
{
return
u
&
0x80000000UL
;
}
uint32
loBit
(
const
uint32
&
u
)
const
{
return
u
&
0x00000001UL
;
}
uint32
loBits
(
const
uint32
&
u
)
const
{
return
u
&
0x7fffffffUL
;
}
uint32
mixBits
(
const
uint32
&
u
,
const
uint32
&
v
)
const
{
return
hiBit
(
u
)
|
loBits
(
v
);
}
uint32
twist
(
const
uint32
&
m
,
const
uint32
&
s0
,
const
uint32
&
s1
)
const
{
return
m
^
(
mixBits
(
s0
,
s1
)
>>
1
)
^
(
-
loBit
(
s1
)
&
0x9908b0dfUL
);
}
static
uint32
hash
(
time_t
t
,
clock_t
c
);
};
inline
MTRand
::
MTRand
(
const
uint32
&
oneSeed
)
{
seed
(
oneSeed
);
}
inline
MTRand
::
MTRand
(
const
uint32
&
oneSeed
)
{
seed
(
oneSeed
);
}
inline
MTRand
::
MTRand
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
)
{
seed
(
bigSeed
,
seedLength
);
}
inline
MTRand
::
MTRand
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
)
{
seed
(
bigSeed
,
seedLength
);
}
inline
MTRand
::
MTRand
()
{
seed
();
}
inline
MTRand
::
MTRand
()
{
seed
();
}
inline
double
MTRand
::
rand
()
{
return
double
(
randInt
())
*
(
1.0
/
4294967295.0
);
}
inline
double
MTRand
::
rand
()
{
return
double
(
randInt
())
*
(
1.0
/
4294967295.0
);
}
inline
double
MTRand
::
rand
(
const
double
&
n
)
{
return
rand
()
*
n
;
}
inline
double
MTRand
::
rand
(
const
double
&
n
)
{
return
rand
()
*
n
;
}
inline
double
MTRand
::
randExc
()
{
return
double
(
randInt
())
*
(
1.0
/
4294967296.0
);
}
inline
double
MTRand
::
randExc
()
{
return
double
(
randInt
())
*
(
1.0
/
4294967296.0
);
}
inline
double
MTRand
::
randExc
(
const
double
&
n
)
{
return
randExc
()
*
n
;
}
inline
double
MTRand
::
randExc
(
const
double
&
n
)
{
return
randExc
()
*
n
;
}
inline
double
MTRand
::
randDblExc
()
{
return
(
double
(
randInt
())
+
0.5
)
*
(
1.0
/
4294967296.0
);
}
inline
double
MTRand
::
randDblExc
()
{
return
(
double
(
randInt
())
+
0.5
)
*
(
1.0
/
4294967296.0
);
}
inline
double
MTRand
::
randDblExc
(
const
double
&
n
)
{
return
randDblExc
()
*
n
;
}
inline
double
MTRand
::
randDblExc
(
const
double
&
n
)
{
return
randDblExc
()
*
n
;
}
inline
double
MTRand
::
rand53
()
{
inline
double
MTRand
::
rand53
()
{
uint32
a
=
randInt
()
>>
5
,
b
=
randInt
()
>>
6
;
return
(
a
*
67108864.0
+
b
)
*
(
1.0
/
9007199254740992.0
);
// by Isaku Wada
return
(
a
*
67108864.0
+
b
)
*
(
1.0
/
9007199254740992.0
);
// by Isaku Wada
}
inline
double
MTRand
::
randNorm
(
const
double
&
mean
,
const
double
&
variance
)
{
inline
double
MTRand
::
randNorm
(
const
double
&
mean
,
const
double
&
variance
)
{
// Return a real number from a normal (Gaussian) distribution with given
// mean and variance by Box-Muller method
double
r
=
sqrt
(
-
2.0
*
log
(
1.0
-
randDblExc
())
)
*
variance
;
double
r
=
sqrt
(
-
2.0
*
log
(
1.0
-
randDblExc
()))
*
variance
;
double
phi
=
2.0
*
3.14159265358979323846264338328
*
randExc
();
return
mean
+
r
*
cos
(
phi
);
}
inline
MTRand
::
uint32
MTRand
::
randInt
()
{
inline
MTRand
::
uint32
MTRand
::
randInt
()
{
// Pull a 32-bit integer from the generator state
// Every other access function simply transforms the numbers extracted here
if
(
left
==
0
)
reload
();
if
(
left
==
0
)
reload
();
--
left
;
uint32
s1
;
s1
=
*
pNext
++
;
s1
^=
(
s1
>>
11
);
s1
^=
(
s1
<<
7
)
&
0x9d2c5680UL
;
s1
^=
(
s1
<<
7
)
&
0x9d2c5680UL
;
s1
^=
(
s1
<<
15
)
&
0xefc60000UL
;
return
(
s1
^
(
s1
>>
18
)
);
return
(
s1
^
(
s1
>>
18
));
}
inline
MTRand
::
uint32
MTRand
::
randInt
(
const
uint32
&
n
)
{
inline
MTRand
::
uint32
MTRand
::
randInt
(
const
uint32
&
n
)
{
// Find which bits are used in n
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
uint32
used
=
n
;
...
...
@@ -214,27 +231,23 @@ inline MTRand::uint32 MTRand::randInt( const uint32& n )
used
|=
used
>>
4
;
used
|=
used
>>
8
;
used
|=
used
>>
16
;
// Draw numbers until one is found in [0,n]
uint32
i
;
do
i
=
randInt
()
&
used
;
// toss unused bits to shorten search
while
(
i
>
n
);
i
=
randInt
()
&
used
;
// toss unused bits to shorten search
while
(
i
>
n
);
return
i
;
}
inline
void
MTRand
::
seed
(
const
uint32
oneSeed
)
{
inline
void
MTRand
::
seed
(
const
uint32
oneSeed
)
{
// Seed the generator with a simple uint32
myoneseed
=
oneSeed
;
myoneseed
=
oneSeed
;
initialize
(
oneSeed
);
reload
();
}
inline
void
MTRand
::
seed
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
)
{
inline
void
MTRand
::
seed
(
uint32
*
const
bigSeed
,
const
uint32
seedLength
)
{
// Seed the generator with an array of uint32's
// There are 2^19937-1 possible initial states. This function allows
// all of those to be accessed by providing at least 19937 bits (with a
...
...
@@ -244,40 +257,42 @@ inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
initialize
(
19650218UL
);
int
i
=
1
;
uint32
j
=
0
;
int
k
=
(
N
>
seedLength
?
N
:
seedLength
);
for
(
;
k
;
--
k
)
{
state
[
i
]
=
state
[
i
]
^
(
(
state
[
i
-
1
]
^
(
state
[
i
-
1
]
>>
30
))
*
1664525UL
);
state
[
i
]
+=
(
bigSeed
[
j
]
&
0xffffffffUL
)
+
j
;
int
k
=
(
N
>
seedLength
?
N
:
seedLength
);
for
(;
k
;
--
k
)
{
state
[
i
]
=
state
[
i
]
^
((
state
[
i
-
1
]
^
(
state
[
i
-
1
]
>>
30
))
*
1664525UL
);
state
[
i
]
+=
(
bigSeed
[
j
]
&
0xffffffffUL
)
+
j
;
state
[
i
]
&=
0xffffffffUL
;
++
i
;
++
j
;
if
(
i
>=
N
)
{
state
[
0
]
=
state
[
N
-
1
];
i
=
1
;
}
if
(
j
>=
seedLength
)
j
=
0
;
++
i
;
++
j
;
if
(
i
>=
N
)
{
state
[
0
]
=
state
[
N
-
1
];
i
=
1
;
}
if
(
j
>=
seedLength
)
j
=
0
;
}
for
(
k
=
N
-
1
;
k
;
--
k
)
{
state
[
i
]
=
state
[
i
]
^
(
(
state
[
i
-
1
]
^
(
state
[
i
-
1
]
>>
30
))
*
1566083941UL
);
for
(
k
=
N
-
1
;
k
;
--
k
)
{
state
[
i
]
=
state
[
i
]
^
((
state
[
i
-
1
]
^
(
state
[
i
-
1
]
>>
30
))
*
1566083941UL
);
state
[
i
]
-=
i
;
state
[
i
]
&=
0xffffffffUL
;
++
i
;
if
(
i
>=
N
)
{
state
[
0
]
=
state
[
N
-
1
];
i
=
1
;
}
if
(
i
>=
N
)
{
state
[
0
]
=
state
[
N
-
1
];
i
=
1
;
}
}
state
[
0
]
=
0x80000000UL
;
// MSB is 1, assuring non-zero initial array
state
[
0
]
=
0x80000000UL
;
// MSB is 1, assuring non-zero initial array
reload
();
//myoneseed=-1;
//
myoneseed=-1;
}
inline
void
MTRand
::
seed
()
{
inline
void
MTRand
::
seed
()
{
// Seed the generator with an array from /dev/urandom if available
// Otherwise use a hash of time() and clock() values
// First try getting an array from /dev/urandom
//FILE* urandom = fopen( "/dev/urandom", "rb" );
//if( urandom )
//
FILE* urandom = fopen( "/dev/urandom", "rb" );
//
if( urandom )
//{
// uint32 bigSeed[N];
// uint32 *s = bigSeed;
...
...
@@ -288,15 +303,13 @@ inline void MTRand::seed()
// fclose(urandom);
// if( success ) { seed( bigSeed, N ); return; }
//}
// Was not successful, so use time() and clock() instead
myoneseed
=
hash
(
time
(
NULL
),
clock
()
);
seed
(
myoneseed
);
myoneseed
=
hash
(
time
(
NULL
),
clock
());
seed
(
myoneseed
);
}
inline
void
MTRand
::
initialize
(
const
uint32
seed
)
{
inline
void
MTRand
::
initialize
(
const
uint32
seed
)
{
// Initialize generator state with seed
// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
// In previous versions, most significant bits (MSBs) of the seed affect
...
...
@@ -305,25 +318,22 @@ inline void MTRand::initialize( const uint32 seed )
uint32
*
r
=
state
;
int
i
=
1
;
*
s
++
=
seed
&
0xffffffffUL
;
for
(
;
i
<
N
;
++
i
)
{
*
s
++
=
(
1812433253UL
*
(
*
r
^
(
*
r
>>
30
)
)
+
i
)
&
0xffffffffUL
;
for
(;
i
<
N
;
++
i
)
{
*
s
++
=
(
1812433253UL
*
(
*
r
^
(
*
r
>>
30
))
+
i
)
&
0xffffffffUL
;
r
++
;
}
}
inline
void
MTRand
::
reload
()
{
inline
void
MTRand
::
reload
()
{
// Generate N new values in state
// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
uint32
*
p
=
state
;
int
i
;
for
(
i
=
N
-
M
;
i
--
;
++
p
)
*
p
=
twist
(
p
[
M
],
p
[
0
],
p
[
1
]
);
for
(
i
=
M
;
--
i
;
++
p
)
*
p
=
twist
(
p
[
M
-
N
],
p
[
0
],
p
[
1
]
);
*
p
=
twist
(
p
[
M
-
N
],
p
[
0
],
state
[
0
]
);
for
(
i
=
N
-
M
;
i
--
;
++
p
)
*
p
=
twist
(
p
[
M
],
p
[
0
],
p
[
1
]);
for
(
i
=
M
;
--
i
;
++
p
)
*
p
=
twist
(
p
[
M
-
N
],
p
[
0
],
p
[
1
]);
*
p
=
twist
(
p
[
M
-
N
],
p
[
0
],
state
[
0
]);
left
=
N
,
pNext
=
state
;
}
...
...
@@ -341,70 +351,62 @@ static uint32_t get_rank() {
return
0
;
}
inline
MTRand
::
uint32
MTRand
::
hash
(
time_t
t
,
clock_t
c
)
{
inline
MTRand
::
uint32
MTRand
::
hash
(
time_t
t
,
clock_t
c
)
{
// Get a uint32 from t and c
// Better than uint32(x) in case x is floating point in [0,1]
// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
static
uint32
differ
=
0
;
// guarantee time-based seeds will change
static
uint32
differ
=
0
;
// guarantee time-based seeds will change
differ
=
get_rank
();
uint32
h1
=
0
;
unsigned
char
*
p
=
(
unsigned
char
*
)
&
t
;
for
(
size_t
i
=
0
;
i
<
sizeof
(
t
);
++
i
)
{
unsigned
char
*
p
=
(
unsigned
char
*
)
&
t
;
for
(
size_t
i
=
0
;
i
<
sizeof
(
t
);
++
i
)
{
h1
*=
UCHAR_MAX
+
2U
;
h1
+=
p
[
i
];
}
uint32
h2
=
0
;
p
=
(
unsigned
char
*
)
&
c
;
for
(
size_t
j
=
0
;
j
<
sizeof
(
c
);
++
j
)
{
p
=
(
unsigned
char
*
)
&
c
;
for
(
size_t
j
=
0
;
j
<
sizeof
(
c
);
++
j
)
{
h2
*=
UCHAR_MAX
+
2U
;
h2
+=
p
[
j
];
}
return
(
h1
+
differ
++
)
^
h2
;
return
(
h1
+
differ
++
)
^
h2
;
}
inline
void
MTRand
::
save
(
std
::
vector
<
uint32
>&
saveArray
)
const
{
saveArray
=
std
::
vector
<
uint32
>
{
state
,
state
+
N
};
inline
void
MTRand
::
save
(
std
::
vector
<
uint32
>
&
saveArray
)
const
{
saveArray
=
std
::
vector
<
uint32
>
{
state
,
state
+
N
};
saveArray
.
push_back
(
left
);
}
inline
void
MTRand
::
load
(
const
std
::
vector
<
uint32
>&
loadArray
)
{
inline
void
MTRand
::
load
(
const
std
::
vector
<
uint32
>
&
loadArray
)
{
left
=
loadArray
.
back
();
assert
(
loadArray
.
size
()
==
N
+
1
);
memcpy
(
state
,
loadArray
.
data
(),
sizeof
(
uint32
)
*
N
);
assert
(
loadArray
.
size
()
==
N
+
1
);
memcpy
(
state
,
loadArray
.
data
(),
sizeof
(
uint32
)
*
N
);
pNext
=
&
state
[
N
-
left
];
pNext
=
&
state
[
N
-
left
];
}
inline
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
MTRand
&
mtrand
)
{
inline
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
MTRand
&
mtrand
)
{
const
MTRand
::
uint32
*
s
=
mtrand
.
state
;
int
i
=
mtrand
.
N
;
for
(
;
i
--
;
os
<<
*
s
++
<<
"
\t
"
)
{}
for
(;
i
--
;
os
<<
*
s
++
<<
"
\t
"
)
{
}
return
os
<<
mtrand
.
left
;
}
inline
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
MTRand
&
mtrand
)
{
inline
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
MTRand
&
mtrand
)
{
MTRand
::
uint32
*
s
=
mtrand
.
state
;
int
i
=
mtrand
.
N
;
for
(
;
i
--
;
is
>>
*
s
++
)
{}
for
(;
i
--
;
is
>>
*
s
++
)
{
}
is
>>
mtrand
.
left
;
mtrand
.
pNext
=
&
mtrand
.
state
[
mtrand
.
N
-
mtrand
.
left
];
mtrand
.
pNext
=
&
mtrand
.
state
[
mtrand
.
N
-
mtrand
.
left
];
return
is
;
}
#endif
// MERSENNETWISTER_H
#endif // MERSENNETWISTER_H
// Change log:
//
...
...
@@ -447,5 +449,3 @@ inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
// - Changed license from GNU LGPL to BSD
//
// Added function get_seed, that is needed for our Monte-Carlo class (S.W.)
acml_rand.h
View file @
0454c645
#ifndef MCL_ACML_RAND_H
#ifndef MCL_ACML_RAND_H
#define MCL_ACML_RAND_H
extern
"C"
void
dranduniform_
(
int
*
N
,
double
*
A
,
double
*
B
,
int
*
state
,
double
*
X
,
int
*
info
);
extern
"C"
void
dranddiscreteuniform_
(
int
*
N
,
int
*
A
,
int
*
B
,
int
*
state
,
int
*
X
,
int
*
info
);
extern
"C"
void
drandinitialize_
(
int
*
GENID
,
int
*
SUBID
,
int
*
SEED
,
int
*
LSEED
,
int
*
STATE
,
int
*
LSTATE
,
int
*
INFO
);
extern
"C"
void
dranduniform_
(
int
*
N
,
double
*
A
,
double
*
B
,
int
*
state
,
double
*
X
,
int
*
info
);
extern
"C"
void
dranddiscreteuniform_
(
int
*
N
,
int
*
A
,
int
*
B
,
int
*
state
,
int
*
X
,
int
*
info
);
extern
"C"
void
drandinitialize_
(
int
*
GENID
,
int
*
SUBID
,
int
*
SEED
,
int
*
LSEED
,
int
*
STATE
,
int
*
LSTATE
,
int
*
INFO
);
#include
<cstdlib>
class
acml_rand
{
public:
acml_rand
(
uint
seed
,
double
A
,
double
B
,
int
len
=
1000
)
{
len_
=
len
;
vec
=
new
double
[
len_
];
A_
=
A
;
B_
=
B
;
//fill seed array
seed_len
=
624
;
int
*
seed_array
=
new
int
[
seed_len
];
srand48
(
seed
);
for
(
uint
i
=
0
;
i
<
(
uint
)
seed_len
;
++
i
)
{
seed_array
[
i
]
=
(
int
)
lrand48
();
public:
acml_rand
(
uint
seed
,
double
A
,
double
B
,
int
len
=
1000
)
{
len_
=
len
;
vec
=
new
double
[
len_
];
A_
=
A
;
B_
=
B
;
// fill seed array
seed_len
=
624
;
int
*
seed_array
=
new
int
[
seed_len
];
srand48
(
seed
);
for
(
uint
i
=
0
;
i
<
(
uint
)
seed_len
;
++
i
)
{
seed_array
[
i
]
=
(
int
)
lrand48
();
}
// initialize Mersenne PRNG int genid = 3;
int
info
;
int
ignored
=
0
;
int
genid
=
3
;
state_len
=
700
;
state
=
new
int
[
state_len
];
drandinitialize_
(
&
genid
,
&
ignored
,
seed_array
,
&
seed_len
,
state
,
&
state_len
,
&
info
);
generate
();
}
//initialize Mersenne PRNG int genid = 3;
int
info
;
int
ignored
=
0
;
int
genid
=
3
;
state_len
=
700
;
state
=
new
int
[
state_len
];
drandinitialize_
(
&
genid
,
&
ignored
,
seed_array
,
&
seed_len
,
state
,
&
state_len
,
&
info
);
generate
();
}
~
acml_rand
()
{
delete
[]
vec
;
delete
[]
state
;
}
inline
const
double
&
operator
()
(
void
)
{
if
(
pos
<
len_
)
return
vec
[
pos
++
];
generate
();
return
vec
[
pos
++
];
}
private: