summaryrefslogtreecommitdiff
path: root/boost/random/inversive_congruential.hpp
blob: 2329c0db3a909b9f7551692184547ee8951fd53d (plain)
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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
/* boost random/inversive_congruential.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP
#define BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP

#include <iosfwd>
#include <stdexcept>
#include <boost/assert.hpp>
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/integer/static_log2.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/const_mod.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/seed_impl.hpp>

#include <boost/random/detail/disable_warnings.hpp>

namespace boost {
namespace random {

// Eichenauer and Lehn 1986
/**
 * Instantiations of class template @c inversive_congruential_engine model a
 * \pseudo_random_number_generator. It uses the inversive congruential
 * algorithm (ICG) described in
 *
 *  @blockquote
 *  "Inversive pseudorandom number generators: concepts, results and links",
 *  Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation
 *  Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman
 *  (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps
 *  @endblockquote
 *
 * The output sequence is defined by x(n+1) = (a*inv(x(n)) - b) (mod p),
 * where x(0), a, b, and the prime number p are parameters of the generator.
 * The expression inv(k) denotes the multiplicative inverse of k in the
 * field of integer numbers modulo p, with inv(0) := 0.
 *
 * The template parameter IntType shall denote a signed integral type large
 * enough to hold p; a, b, and p are the parameters of the generators. The
 * template parameter val is the validation value checked by validation.
 *
 * @xmlnote
 * The implementation currently uses the Euclidian Algorithm to compute
 * the multiplicative inverse. Therefore, the inversive generators are about
 * 10-20 times slower than the others (see section"performance"). However,
 * the paper talks of only 3x slowdown, so the Euclidian Algorithm is probably
 * not optimal for calculating the multiplicative inverse.
 * @endxmlnote
 */
template<class IntType, IntType a, IntType b, IntType p>
class inversive_congruential_engine
{
public:
    typedef IntType result_type;
    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

    BOOST_STATIC_CONSTANT(result_type, multiplier = a);
    BOOST_STATIC_CONSTANT(result_type, increment = b);
    BOOST_STATIC_CONSTANT(result_type, modulus = p);
    BOOST_STATIC_CONSTANT(IntType, default_seed = 1);

    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return b == 0 ? 1 : 0; }
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return p-1; }
    
    /**
     * Constructs an @c inversive_congruential_engine, seeding it with
     * the default seed.
     */
    inversive_congruential_engine() { seed(); }

    /**
     * Constructs an @c inversive_congruential_engine, seeding it with @c x0.
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(inversive_congruential_engine,
                                               IntType, x0)
    { seed(x0); }
    
    /**
     * Constructs an @c inversive_congruential_engine, seeding it with values
     * produced by a call to @c seq.generate().
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(inversive_congruential_engine,
                                             SeedSeq, seq)
    { seed(seq); }
    
    /**
     * Constructs an @c inversive_congruential_engine, seeds it
     * with values taken from the itrator range [first, last),
     * and adjusts first to point to the element after the last one
     * used.  If there are not enough elements, throws @c std::invalid_argument.
     *
     * first and last must be input iterators.
     */
    template<class It> inversive_congruential_engine(It& first, It last)
    { seed(first, last); }

    /**
     * Calls seed(default_seed)
     */
    void seed() { seed(default_seed); }
  
    /**
     * If c mod m is zero and x0 mod m is zero, changes the current value of
     * the generator to 1. Otherwise, changes it to x0 mod m. If c is zero,
     * distinct seeds in the range [1,m) will leave the generator in distinct
     * states. If c is not zero, the range is [0,m).
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(inversive_congruential_engine, IntType, x0)
    {
        // wrap _x if it doesn't fit in the destination
        if(modulus == 0) {
            _value = x0;
        } else {
            _value = x0 % modulus;
        }
        // handle negative seeds
        if(_value <= 0 && _value != 0) {
            _value += modulus;
        }
        // adjust to the correct range
        if(increment == 0 && _value == 0) {
            _value = 1;
        }
        BOOST_ASSERT(_value >= (min)());
        BOOST_ASSERT(_value <= (max)());
    }

    /**
     * Seeds an @c inversive_congruential_engine using values from a SeedSeq.
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(inversive_congruential_engine, SeedSeq, seq)
    { seed(detail::seed_one_int<IntType, modulus>(seq)); }
    
    /**
     * seeds an @c inversive_congruential_engine with values taken
     * from the itrator range [first, last) and adjusts @c first to
     * point to the element after the last one used.  If there are
     * not enough elements, throws @c std::invalid_argument.
     *
     * @c first and @c last must be input iterators.
     */
    template<class It> void seed(It& first, It last)
    { seed(detail::get_one_int<IntType, modulus>(first, last)); }

    /** Returns the next output of the generator. */
    IntType operator()()
    {
        typedef const_mod<IntType, p> do_mod;
        _value = do_mod::mult_add(a, do_mod::invert(_value), b);
        return _value;
    }
  
    /** Fills a range with random values */
    template<class Iter>
    void generate(Iter first, Iter last)
    { detail::generate_from_int(*this, first, last); }

    /** Advances the state of the generator by @c z. */
    void discard(boost::uintmax_t z)
    {
        for(boost::uintmax_t j = 0; j < z; ++j) {
            (*this)();
        }
    }

    /**
     * Writes the textual representation of the generator to a @c std::ostream.
     */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, inversive_congruential_engine, x)
    {
        os << x._value;
        return os;
    }

    /**
     * Reads the textual representation of the generator from a @c std::istream.
     */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, inversive_congruential_engine, x)
    {
        is >> x._value;
        return is;
    }

    /**
     * Returns true if the two generators will produce identical
     * sequences of outputs.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(inversive_congruential_engine, x, y)
    { return x._value == y._value; }

    /**
     * Returns true if the two generators will produce different
     * sequences of outputs.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(inversive_congruential_engine)

private:
    IntType _value;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class IntType, IntType a, IntType b, IntType p>
const bool inversive_congruential_engine<IntType, a, b, p>::has_fixed_range;
template<class IntType, IntType a, IntType b, IntType p>
const typename inversive_congruential_engine<IntType, a, b, p>::result_type inversive_congruential_engine<IntType, a, b, p>::multiplier;
template<class IntType, IntType a, IntType b, IntType p>
const typename inversive_congruential_engine<IntType, a, b, p>::result_type inversive_congruential_engine<IntType, a, b, p>::increment;
template<class IntType, IntType a, IntType b, IntType p>
const typename inversive_congruential_engine<IntType, a, b, p>::result_type inversive_congruential_engine<IntType, a, b, p>::modulus;
template<class IntType, IntType a, IntType b, IntType p>
const typename inversive_congruential_engine<IntType, a, b, p>::result_type inversive_congruential_engine<IntType, a, b, p>::default_seed;
#endif

/// \cond show_deprecated

// provided for backwards compatibility
template<class IntType, IntType a, IntType b, IntType p, IntType val = 0>
class inversive_congruential : public inversive_congruential_engine<IntType, a, b, p>
{
    typedef inversive_congruential_engine<IntType, a, b, p> base_type;
public:
    inversive_congruential(IntType x0 = 1) : base_type(x0) {}
    template<class It>
    inversive_congruential(It& first, It last) : base_type(first, last) {}
};

/// \endcond

/**
 * The specialization hellekalek1995 was suggested in
 *
 *  @blockquote
 *  "Inversive pseudorandom number generators: concepts, results and links",
 *  Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation
 *  Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman
 *  (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps
 *  @endblockquote
 */
typedef inversive_congruential_engine<uint32_t, 9102, 2147483647-36884165,
  2147483647> hellekalek1995;

} // namespace random

using random::hellekalek1995;

} // namespace boost

#include <boost/random/detail/enable_warnings.hpp>

#endif // BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP