1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
// -*- C++ -*-
// Copyright (C) 2007, 2008 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library. This library is free
// software; you can redistribute it and/or modify it under the terms
// of the GNU General Public License as published by the Free Software
// Foundation; either version 2, or (at your option) any later
// version.
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this library; see the file COPYING. If not, write to
// the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.
// As a special exception, you may use this file as part of a free
// software library without restriction. Specifically, if other files
// instantiate templates or use macros or inline functions from this
// file, or you compile this file and link it with other files to
// produce an executable, this file does not by itself cause the
// resulting executable to be covered by the GNU General Public
// License. This exception does not however invalidate any other
// reasons why the executable file might be covered by the GNU General
// Public License.
/** @file parallel/random_number.h
* @brief Random number generator based on the Mersenne twister.
* This file is a GNU parallel extension to the Standard C++ Library.
*/
// Written by Johannes Singler.
#ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
#define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
#include <parallel/types.h>
#include <tr1/random>
namespace __gnu_parallel
{
/** @brief Random number generator, based on the Mersenne twister. */
class random_number
{
private:
std::tr1::mt19937 mt;
uint64 supremum;
uint64 RAND_SUP;
double supremum_reciprocal;
double RAND_SUP_REC;
// Assumed to be twice as long as the usual random number.
uint64 cache;
// Bit results.
int bits_left;
static uint32
scale_down(uint64 x,
#if _GLIBCXX_SCALE_DOWN_FPU
uint64 /*supremum*/, double supremum_reciprocal)
#else
uint64 supremum, double /*supremum_reciprocal*/)
#endif
{
#if _GLIBCXX_SCALE_DOWN_FPU
return uint32(x * supremum_reciprocal);
#else
return static_cast<uint32>(x % supremum);
#endif
}
public:
/** @brief Default constructor. Seed with 0. */
random_number()
: mt(0), supremum(0x100000000ULL),
RAND_SUP(1ULL << (sizeof(uint32) * 8)),
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
RAND_SUP_REC(1.0 / double(RAND_SUP)),
cache(0), bits_left(0) { }
/** @brief Constructor.
* @param seed Random seed.
* @param supremum Generate integer random numbers in the
* interval @c [0,supremum). */
random_number(uint32 seed, uint64 supremum = 0x100000000ULL)
: mt(seed), supremum(supremum),
RAND_SUP(1ULL << (sizeof(uint32) * 8)),
supremum_reciprocal(double(supremum) / double(RAND_SUP)),
RAND_SUP_REC(1.0 / double(RAND_SUP)),
cache(0), bits_left(0) { }
/** @brief Generate unsigned random 32-bit integer. */
uint32
operator()()
{ return scale_down(mt(), supremum, supremum_reciprocal); }
/** @brief Generate unsigned random 32-bit integer in the
interval @c [0,local_supremum). */
uint32
operator()(uint64 local_supremum)
{
return scale_down(mt(), local_supremum,
double(local_supremum * RAND_SUP_REC));
}
/** @brief Generate a number of random bits, run-time parameter.
* @param bits Number of bits to generate. */
unsigned long
genrand_bits(int bits)
{
unsigned long res = cache & ((1 << bits) - 1);
cache = cache >> bits;
bits_left -= bits;
if (bits_left < 32)
{
cache |= ((uint64(mt())) << bits_left);
bits_left += 32;
}
return res;
}
};
} // namespace __gnu_parallel
#endif /* _GLIBCXX_PARALLEL_RANDOM_NUMBER_H */
|
