epoc32/include/stdapis/boost/random/detail/const_mod.hpp
author William Roberts <williamr@symbian.org>
Wed, 31 Mar 2010 12:27:01 +0100
branchSymbian2
changeset 3 e1b950c65cb4
permissions -rw-r--r--
Attempt to represent the S^2->S^3 header reorganisation as a series of "hg rename" operations
williamr@2
     1
/* boost random/detail/const_mod.hpp header file
williamr@2
     2
 *
williamr@2
     3
 * Copyright Jens Maurer 2000-2001
williamr@2
     4
 * Distributed under the Boost Software License, Version 1.0. (See
williamr@2
     5
 * accompanying file LICENSE_1_0.txt or copy at
williamr@2
     6
 * http://www.boost.org/LICENSE_1_0.txt)
williamr@2
     7
 *
williamr@2
     8
 * See http://www.boost.org for most recent version including documentation.
williamr@2
     9
 *
williamr@2
    10
 * $Id: const_mod.hpp,v 1.8 2004/07/27 03:43:32 dgregor Exp $
williamr@2
    11
 *
williamr@2
    12
 * Revision history
williamr@2
    13
 *  2001-02-18  moved to individual header files
williamr@2
    14
 */
williamr@2
    15
williamr@2
    16
#ifndef BOOST_RANDOM_CONST_MOD_HPP
williamr@2
    17
#define BOOST_RANDOM_CONST_MOD_HPP
williamr@2
    18
williamr@2
    19
#include <cassert>
williamr@2
    20
#include <boost/static_assert.hpp>
williamr@2
    21
#include <boost/cstdint.hpp>
williamr@2
    22
#include <boost/integer_traits.hpp>
williamr@2
    23
#include <boost/detail/workaround.hpp>
williamr@2
    24
williamr@2
    25
namespace boost {
williamr@2
    26
namespace random {
williamr@2
    27
williamr@2
    28
/*
williamr@2
    29
 * Some random number generators require modular arithmetic.  Put
williamr@2
    30
 * everything we need here.
williamr@2
    31
 * IntType must be an integral type.
williamr@2
    32
 */
williamr@2
    33
williamr@2
    34
namespace detail {
williamr@2
    35
williamr@2
    36
  template<bool is_signed>
williamr@2
    37
  struct do_add
williamr@2
    38
  { };
williamr@2
    39
williamr@2
    40
  template<>
williamr@2
    41
  struct do_add<true>
williamr@2
    42
  {
williamr@2
    43
    template<class IntType>
williamr@2
    44
    static IntType add(IntType m, IntType x, IntType c)
williamr@2
    45
    {
williamr@2
    46
      x += (c-m);
williamr@2
    47
      if(x < 0)
williamr@2
    48
        x += m;
williamr@2
    49
      return x;
williamr@2
    50
    }
williamr@2
    51
  };
williamr@2
    52
williamr@2
    53
  template<>
williamr@2
    54
  struct do_add<false>
williamr@2
    55
  {
williamr@2
    56
    template<class IntType>
williamr@2
    57
    static IntType add(IntType, IntType, IntType)
williamr@2
    58
    {
williamr@2
    59
      // difficult
williamr@2
    60
      assert(!"const_mod::add with c too large");
williamr@2
    61
      return 0;
williamr@2
    62
    }
williamr@2
    63
  };
williamr@2
    64
} // namespace detail
williamr@2
    65
williamr@2
    66
#if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560))
williamr@2
    67
williamr@2
    68
template<class IntType, IntType m>
williamr@2
    69
class const_mod
williamr@2
    70
{
williamr@2
    71
public:
williamr@2
    72
  static IntType add(IntType x, IntType c)
williamr@2
    73
  {
williamr@2
    74
    if(c == 0)
williamr@2
    75
      return x;
williamr@2
    76
    else if(c <= traits::const_max - m)    // i.e. m+c < max
williamr@2
    77
      return add_small(x, c);
williamr@2
    78
    else
williamr@2
    79
      return detail::do_add<traits::is_signed>::add(m, x, c);
williamr@2
    80
  }
williamr@2
    81
williamr@2
    82
  static IntType mult(IntType a, IntType x)
williamr@2
    83
  {
williamr@2
    84
    if(a == 1)
williamr@2
    85
      return x;
williamr@2
    86
    else if(m <= traits::const_max/a)      // i.e. a*m <= max
williamr@2
    87
      return mult_small(a, x);
williamr@2
    88
    else if(traits::is_signed && (m%a < m/a))
williamr@2
    89
      return mult_schrage(a, x);
williamr@2
    90
    else {
williamr@2
    91
      // difficult
williamr@2
    92
      assert(!"const_mod::mult with a too large");
williamr@2
    93
      return 0;
williamr@2
    94
    }
williamr@2
    95
  }
williamr@2
    96
williamr@2
    97
  static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2
    98
  {
williamr@2
    99
    if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
williamr@2
   100
      return (a*x+c) % m;
williamr@2
   101
    else
williamr@2
   102
      return add(mult(a, x), c);
williamr@2
   103
  }
williamr@2
   104
williamr@2
   105
  static IntType invert(IntType x)
williamr@2
   106
  { return x == 0 ? 0 : invert_euclidian(x); }
williamr@2
   107
williamr@2
   108
private:
williamr@2
   109
  typedef integer_traits<IntType> traits;
williamr@2
   110
williamr@2
   111
  const_mod();      // don't instantiate
williamr@2
   112
williamr@2
   113
  static IntType add_small(IntType x, IntType c)
williamr@2
   114
  {
williamr@2
   115
    x += c;
williamr@2
   116
    if(x >= m)
williamr@2
   117
      x -= m;
williamr@2
   118
    return x;
williamr@2
   119
  }
williamr@2
   120
williamr@2
   121
  static IntType mult_small(IntType a, IntType x)
williamr@2
   122
  {
williamr@2
   123
    return a*x % m;
williamr@2
   124
  }
williamr@2
   125
williamr@2
   126
  static IntType mult_schrage(IntType a, IntType value)
williamr@2
   127
  {
williamr@2
   128
    const IntType q = m / a;
williamr@2
   129
    const IntType r = m % a;
williamr@2
   130
williamr@2
   131
    assert(r < q);        // check that overflow cannot happen
williamr@2
   132
williamr@2
   133
    value = a*(value%q) - r*(value/q);
williamr@2
   134
    // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this
williamr@2
   135
    // convoluted formulation of the loop (Synge Todo)
williamr@2
   136
    for(;;) {
williamr@2
   137
      if (value > 0)
williamr@2
   138
        break;
williamr@2
   139
      value += m;
williamr@2
   140
    }
williamr@2
   141
    return value;
williamr@2
   142
  }
williamr@2
   143
williamr@2
   144
  // invert c in the finite field (mod m) (m must be prime)
williamr@2
   145
  static IntType invert_euclidian(IntType c)
williamr@2
   146
  {
williamr@2
   147
    // we are interested in the gcd factor for c, because this is our inverse
williamr@2
   148
    BOOST_STATIC_ASSERT(m > 0);
williamr@2
   149
#if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003))
williamr@2
   150
    assert(boost::integer_traits<IntType>::is_signed);
williamr@2
   151
#elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS)
williamr@2
   152
    BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
williamr@2
   153
#endif
williamr@2
   154
    assert(c > 0);
williamr@2
   155
    IntType l1 = 0;
williamr@2
   156
    IntType l2 = 1;
williamr@2
   157
    IntType n = c;
williamr@2
   158
    IntType p = m;
williamr@2
   159
    for(;;) {
williamr@2
   160
      IntType q = p / n;
williamr@2
   161
      l1 -= q * l2;           // this requires a signed IntType!
williamr@2
   162
      p -= q * n;
williamr@2
   163
      if(p == 0)
williamr@2
   164
        return (l2 < 1 ? l2 + m : l2);
williamr@2
   165
      IntType q2 = n / p;
williamr@2
   166
      l2 -= q2 * l1;
williamr@2
   167
      n -= q2 * p;
williamr@2
   168
      if(n == 0)
williamr@2
   169
        return (l1 < 1 ? l1 + m : l1);
williamr@2
   170
    }
williamr@2
   171
  }
williamr@2
   172
};
williamr@2
   173
williamr@2
   174
// The modulus is exactly the word size: rely on machine overflow handling.
williamr@2
   175
// Due to a GCC bug, we cannot partially specialize in the presence of
williamr@2
   176
// template value parameters.
williamr@2
   177
template<>
williamr@2
   178
class const_mod<unsigned int, 0>
williamr@2
   179
{
williamr@2
   180
  typedef unsigned int IntType;
williamr@2
   181
public:
williamr@2
   182
  static IntType add(IntType x, IntType c) { return x+c; }
williamr@2
   183
  static IntType mult(IntType a, IntType x) { return a*x; }
williamr@2
   184
  static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
williamr@2
   185
williamr@2
   186
  // m is not prime, thus invert is not useful
williamr@2
   187
private:                      // don't instantiate
williamr@2
   188
  const_mod();
williamr@2
   189
};
williamr@2
   190
williamr@2
   191
template<>
williamr@2
   192
class const_mod<unsigned long, 0>
williamr@2
   193
{
williamr@2
   194
  typedef unsigned long IntType;
williamr@2
   195
public:
williamr@2
   196
  static IntType add(IntType x, IntType c) { return x+c; }
williamr@2
   197
  static IntType mult(IntType a, IntType x) { return a*x; }
williamr@2
   198
  static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
williamr@2
   199
williamr@2
   200
  // m is not prime, thus invert is not useful
williamr@2
   201
private:                      // don't instantiate
williamr@2
   202
  const_mod();
williamr@2
   203
};
williamr@2
   204
williamr@2
   205
// the modulus is some power of 2: rely partly on machine overflow handling
williamr@2
   206
// we only specialize for rand48 at the moment
williamr@2
   207
#ifndef BOOST_NO_INT64_T
williamr@2
   208
template<>
williamr@2
   209
class const_mod<uint64_t, uint64_t(1) << 48>
williamr@2
   210
{
williamr@2
   211
  typedef uint64_t IntType;
williamr@2
   212
public:
williamr@2
   213
  static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
williamr@2
   214
  static IntType mult(IntType a, IntType x) { return mod(a*x); }
williamr@2
   215
  static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2
   216
    { return mod(a*x+c); }
williamr@2
   217
  static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
williamr@2
   218
williamr@2
   219
  // m is not prime, thus invert is not useful
williamr@2
   220
private:                      // don't instantiate
williamr@2
   221
  const_mod();
williamr@2
   222
};
williamr@2
   223
#endif /* !BOOST_NO_INT64_T */
williamr@2
   224
williamr@2
   225
#else
williamr@2
   226
williamr@2
   227
//
williamr@2
   228
// for some reason Borland C++ Builder 6 has problems with
williamr@2
   229
// the full specialisations of const_mod, define a generic version
williamr@2
   230
// instead, the compiler will optimise away the const-if statements:
williamr@2
   231
//
williamr@2
   232
williamr@2
   233
template<class IntType, IntType m>
williamr@2
   234
class const_mod
williamr@2
   235
{
williamr@2
   236
public:
williamr@2
   237
  static IntType add(IntType x, IntType c)
williamr@2
   238
  {
williamr@2
   239
    if(0 == m)
williamr@2
   240
    {
williamr@2
   241
       return x+c;
williamr@2
   242
    }
williamr@2
   243
    else
williamr@2
   244
    {
williamr@2
   245
       if(c == 0)
williamr@2
   246
         return x;
williamr@2
   247
       else if(c <= traits::const_max - m)    // i.e. m+c < max
williamr@2
   248
         return add_small(x, c);
williamr@2
   249
       else
williamr@2
   250
         return detail::do_add<traits::is_signed>::add(m, x, c);
williamr@2
   251
    }
williamr@2
   252
  }
williamr@2
   253
williamr@2
   254
  static IntType mult(IntType a, IntType x)
williamr@2
   255
  {
williamr@2
   256
    if(x == 0)
williamr@2
   257
    {
williamr@2
   258
       return a*x;
williamr@2
   259
    }
williamr@2
   260
    else
williamr@2
   261
    {
williamr@2
   262
       if(a == 1)
williamr@2
   263
         return x;
williamr@2
   264
       else if(m <= traits::const_max/a)      // i.e. a*m <= max
williamr@2
   265
         return mult_small(a, x);
williamr@2
   266
       else if(traits::is_signed && (m%a < m/a))
williamr@2
   267
         return mult_schrage(a, x);
williamr@2
   268
       else {
williamr@2
   269
         // difficult
williamr@2
   270
         assert(!"const_mod::mult with a too large");
williamr@2
   271
         return 0;
williamr@2
   272
       }
williamr@2
   273
    }
williamr@2
   274
  }
williamr@2
   275
williamr@2
   276
  static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2
   277
  {
williamr@2
   278
    if(m == 0)
williamr@2
   279
    {
williamr@2
   280
       return a*x+c;
williamr@2
   281
    }
williamr@2
   282
    else
williamr@2
   283
    {
williamr@2
   284
       if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
williamr@2
   285
         return (a*x+c) % m;
williamr@2
   286
       else
williamr@2
   287
         return add(mult(a, x), c);
williamr@2
   288
    }
williamr@2
   289
  }
williamr@2
   290
williamr@2
   291
  static IntType invert(IntType x)
williamr@2
   292
  { return x == 0 ? 0 : invert_euclidian(x); }
williamr@2
   293
williamr@2
   294
private:
williamr@2
   295
  typedef integer_traits<IntType> traits;
williamr@2
   296
williamr@2
   297
  const_mod();      // don't instantiate
williamr@2
   298
williamr@2
   299
  static IntType add_small(IntType x, IntType c)
williamr@2
   300
  {
williamr@2
   301
    x += c;
williamr@2
   302
    if(x >= m)
williamr@2
   303
      x -= m;
williamr@2
   304
    return x;
williamr@2
   305
  }
williamr@2
   306
williamr@2
   307
  static IntType mult_small(IntType a, IntType x)
williamr@2
   308
  {
williamr@2
   309
    return a*x % m;
williamr@2
   310
  }
williamr@2
   311
williamr@2
   312
  static IntType mult_schrage(IntType a, IntType value)
williamr@2
   313
  {
williamr@2
   314
    const IntType q = m / a;
williamr@2
   315
    const IntType r = m % a;
williamr@2
   316
williamr@2
   317
    assert(r < q);        // check that overflow cannot happen
williamr@2
   318
williamr@2
   319
    value = a*(value%q) - r*(value/q);
williamr@2
   320
    while(value <= 0)
williamr@2
   321
      value += m;
williamr@2
   322
    return value;
williamr@2
   323
  }
williamr@2
   324
williamr@2
   325
  // invert c in the finite field (mod m) (m must be prime)
williamr@2
   326
  static IntType invert_euclidian(IntType c)
williamr@2
   327
  {
williamr@2
   328
    // we are interested in the gcd factor for c, because this is our inverse
williamr@2
   329
    BOOST_STATIC_ASSERT(m > 0);
williamr@2
   330
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
williamr@2
   331
    BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
williamr@2
   332
#endif
williamr@2
   333
    assert(c > 0);
williamr@2
   334
    IntType l1 = 0;
williamr@2
   335
    IntType l2 = 1;
williamr@2
   336
    IntType n = c;
williamr@2
   337
    IntType p = m;
williamr@2
   338
    for(;;) {
williamr@2
   339
      IntType q = p / n;
williamr@2
   340
      l1 -= q * l2;           // this requires a signed IntType!
williamr@2
   341
      p -= q * n;
williamr@2
   342
      if(p == 0)
williamr@2
   343
        return (l2 < 1 ? l2 + m : l2);
williamr@2
   344
      IntType q2 = n / p;
williamr@2
   345
      l2 -= q2 * l1;
williamr@2
   346
      n -= q2 * p;
williamr@2
   347
      if(n == 0)
williamr@2
   348
        return (l1 < 1 ? l1 + m : l1);
williamr@2
   349
    }
williamr@2
   350
  }
williamr@2
   351
};
williamr@2
   352
williamr@2
   353
williamr@2
   354
#endif
williamr@2
   355
williamr@2
   356
} // namespace random
williamr@2
   357
} // namespace boost
williamr@2
   358
williamr@2
   359
#endif // BOOST_RANDOM_CONST_MOD_HPP