Scopira  20080306
constmod.h
1 
2 /*
3  * Copyright (c) 2002 National Research Council
4  *
5  * All rights reserved.
6  *
7  * This material is confidential and proprietary information of
8  * National Research Council Canada ("Confidential Information").
9  * This Confidential Information may only be used and reproduced
10  * in accordance with the terms of the license agreement.
11  *
12  */
13 
14 #ifndef __INCLUDED_SCOPIRA_TOOL_CONSTMOD_H__
15 #define __INCLUDED_SCOPIRA_TOOL_CONSTMOD_H__
16 
17 //
18 // taken from Boost, for random.h
19 //
20 
21 #include <scopira/tool/limits.h>
22 #include <scopira/tool/platform.h>
23 
24 namespace scopira {
25 namespace tool {
26 
27 /*
28  * Some random number generators require modular arithmetic. Put
29  * everything we need here.
30  * IntType must be an integral type.
31  */
32 
33 template<class IntType, IntType m>
34 class const_mod
35 {
36 public:
37  static IntType add(IntType x, IntType c)
38  {
39  if(c == 0)
40  return x;
41  else if(c <= traits::max() - m) // i.e. m+c < max
42  return add_small(x, c);
43  else if(traits::is_signed)
44  return add_signed(x, c);
45  else {
46  // difficult
47  assert(!"const_mod::add with c too large");
48  return 0;
49  }
50  }
51 
52  static IntType mult(IntType a, IntType x)
53  {
54  if(a == 1)
55  return x;
56  else if(m <= traits::max()/a) // i.e. a*m <= max
57  return mult_small(a, x);
58  else if(traits::is_signed && (m%a < m/a))
59  return mult_schrage(a, x);
60  else {
61  // difficult
62  assert(!"const_mod::mult with a too large");
63  return 0;
64  }
65  }
66 
67  static IntType mult_add(IntType a, IntType x, IntType c)
68  {
69  if(m <= (traits::max()-c)/a) // i.e. a*m+c <= max
70  return (a*x+c) % m;
71  else
72  return add(mult(a, x), c);
73  }
74 
75  static IntType invert(IntType x)
76  { return x == 0 ? 0 : invert_euclidian(x); }
77 
78 private:
79  typedef std::numeric_limits<IntType> traits;
80 
81  const_mod(); // don't instantiate
82 
83  static IntType add_small(IntType x, IntType c)
84  {
85  x += c;
86  if(x >= m)
87  x -= m;
88  return x;
89  }
90 
91  static IntType add_signed(IntType x, IntType c)
92  {
93  x += (c-m);
94  if(x < 0)
95  x += m;
96  return x;
97  }
98 
99  static IntType mult_small(IntType a, IntType x)
100  {
101  return a*x % m;
102  }
103 
104  static IntType mult_schrage(IntType a, IntType value)
105  {
106  const IntType q = m / a;
107  const IntType r = m % a;
108 
109  assert(r < q); // check that overflow cannot happen
110 
111  value = a*(value%q) - r*(value/q);
112  while(value <= 0)
113  value += m;
114  return value;
115  }
116 
117  // invert c in the finite field (mod m) (m must be prime)
118  static IntType invert_euclidian(IntType c)
119  {
120  // we are interested in the gcd factor for c, because this is our inverse
121  BOOST_STATIC_ASSERT(m > 0);
122 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
123  BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_signed);
124 #endif
125  assert(c > 0);
126  IntType l1 = 0;
127  IntType l2 = 1;
128  IntType n = c;
129  IntType p = m;
130  for(;;) {
131  IntType q = p / n;
132  l1 -= q * l2; // this requires a signed IntType!
133  p -= q * n;
134  if(p == 0)
135  return (l2 < 1 ? l2 + m : l2);
136  IntType q2 = n / p;
137  l2 -= q2 * l1;
138  n -= q2 * p;
139  if(n == 0)
140  return (l1 < 1 ? l1 + m : l1);
141  }
142  }
143 };
144 
145 // The modulus is exactly the word size: rely on machine overflow handling.
146 // Due to a GCC bug, we cannot partially specialize in the presence of
147 // template value parameters.
148 template<>
149 class const_mod<unsigned int, 0>
150 {
151  typedef unsigned int IntType;
152 public:
153  static IntType add(IntType x, IntType c) { return x+c; }
154  static IntType mult(IntType a, IntType x) { return a*x; }
155  static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
156 
157  // m is not prime, thus invert is not useful
158 private: // don't instantiate
159  const_mod();
160 };
161 
162 template<>
163 class const_mod<unsigned long, 0>
164 {
165  typedef unsigned long IntType;
166 public:
167  static IntType add(IntType x, IntType c) { return x+c; }
168  static IntType mult(IntType a, IntType x) { return a*x; }
169  static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
170 
171  // m is not prime, thus invert is not useful
172 private: // don't instantiate
173  const_mod();
174 };
175 
176 template<>
177 class const_mod<uint64_t, uint64_t(1) << 48>
178 {
179  typedef uint64_t IntType;
180 public:
181  static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
182  static IntType mult(IntType a, IntType x) { return mod(a*x); }
183  static IntType mult_add(IntType a, IntType x, IntType c)
184  { return mod(a*x+c); }
185  static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
186 
187  // m is not prime, thus invert is not useful
188 private: // don't instantiate
189  const_mod();
190 };
191 
192 } // namespace random
193 } // namespace boost
194 
195 #endif
196 
Definition: archiveflow.h:20
Definition: constmod.h:34