Scopira 20080306

constmod.h

00001 
00002 /*
00003  *  Copyright (c) 2002    National Research Council
00004  *
00005  *  All rights reserved.
00006  *
00007  *  This material is confidential and proprietary information of
00008  *  National Research Council Canada ("Confidential Information").
00009  *  This Confidential Information may only be used and reproduced
00010  *  in accordance with the terms of the license agreement.
00011  *
00012  */
00013 
00014 #ifndef __INCLUDED_SCOPIRA_TOOL_CONSTMOD_H__
00015 #define __INCLUDED_SCOPIRA_TOOL_CONSTMOD_H__
00016 
00017 //
00018 // taken from Boost, for random.h
00019 //
00020 
00021 #include <scopira/tool/limits.h>
00022 #include <scopira/tool/platform.h>
00023 
00024 namespace scopira {
00025 namespace tool {
00026 
00027 /*
00028  * Some random number generators require modular arithmetic.  Put
00029  * everything we need here.
00030  * IntType must be an integral type.
00031  */
00032 
00033 template<class IntType, IntType m>
00034 class const_mod
00035 {
00036 public:
00037   static IntType add(IntType x, IntType c)
00038   {
00039     if(c == 0)
00040       return x;
00041     else if(c <= traits::max() - m)    // i.e. m+c < max
00042       return add_small(x, c);
00043     else if(traits::is_signed)
00044       return add_signed(x, c);
00045     else {
00046       // difficult
00047       assert(!"const_mod::add with c too large");
00048       return 0;
00049     }
00050   }
00051 
00052   static IntType mult(IntType a, IntType x)
00053   {
00054     if(a == 1)
00055       return x;
00056     else if(m <= traits::max()/a)      // i.e. a*m <= max
00057       return mult_small(a, x);
00058     else if(traits::is_signed && (m%a < m/a))
00059       return mult_schrage(a, x);
00060     else {
00061       // difficult
00062       assert(!"const_mod::mult with a too large");
00063       return 0;
00064     }
00065   }
00066 
00067   static IntType mult_add(IntType a, IntType x, IntType c)
00068   { 
00069     if(m <= (traits::max()-c)/a)   // i.e. a*m+c <= max
00070       return (a*x+c) % m;
00071     else
00072       return add(mult(a, x), c);
00073   }
00074 
00075   static IntType invert(IntType x)
00076   { return x == 0 ? 0 : invert_euclidian(x); }
00077 
00078 private:
00079   typedef std::numeric_limits<IntType> traits;
00080 
00081   const_mod();      // don't instantiate
00082 
00083   static IntType add_small(IntType x, IntType c)
00084   {
00085     x += c;
00086     if(x >= m)
00087       x -= m;
00088     return x;
00089   }
00090 
00091   static IntType add_signed(IntType x, IntType c)
00092   {
00093     x += (c-m);
00094     if(x < 0)
00095       x += m;
00096     return x;
00097   }
00098 
00099   static IntType mult_small(IntType a, IntType x)
00100   {
00101     return a*x % m;
00102   }
00103 
00104   static IntType mult_schrage(IntType a, IntType value)
00105   {
00106     const IntType q = m / a;
00107     const IntType r = m % a;
00108 
00109     assert(r < q);        // check that overflow cannot happen
00110 
00111     value = a*(value%q) - r*(value/q);
00112     while(value <= 0)
00113       value += m;
00114     return value;
00115   }
00116 
00117   // invert c in the finite field (mod m) (m must be prime)
00118   static IntType invert_euclidian(IntType c)
00119   {
00120     // we are interested in the gcd factor for c, because this is our inverse
00121     BOOST_STATIC_ASSERT(m > 0);
00122 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
00123     BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_signed);
00124 #endif
00125     assert(c > 0);
00126     IntType l1 = 0;
00127     IntType l2 = 1;
00128     IntType n = c;
00129     IntType p = m;
00130     for(;;) {
00131       IntType q = p / n;
00132       l1 -= q * l2;           // this requires a signed IntType!
00133       p -= q * n;
00134       if(p == 0)
00135         return (l2 < 1 ? l2 + m : l2);
00136       IntType q2 = n / p;
00137       l2 -= q2 * l1;
00138       n -= q2 * p;
00139       if(n == 0)
00140         return (l1 < 1 ? l1 + m : l1);
00141     }
00142   }
00143 };
00144 
00145 // The modulus is exactly the word size: rely on machine overflow handling.
00146 // Due to a GCC bug, we cannot partially specialize in the presence of
00147 // template value parameters.
00148 template<>
00149 class const_mod<unsigned int, 0>
00150 {
00151   typedef unsigned int IntType;
00152 public:
00153   static IntType add(IntType x, IntType c) { return x+c; }
00154   static IntType mult(IntType a, IntType x) { return a*x; }
00155   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
00156 
00157   // m is not prime, thus invert is not useful
00158 private:                      // don't instantiate
00159   const_mod();
00160 };
00161 
00162 template<>
00163 class const_mod<unsigned long, 0>
00164 {
00165   typedef unsigned long IntType;
00166 public:
00167   static IntType add(IntType x, IntType c) { return x+c; }
00168   static IntType mult(IntType a, IntType x) { return a*x; }
00169   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
00170 
00171   // m is not prime, thus invert is not useful
00172 private:                      // don't instantiate
00173   const_mod();
00174 };
00175 
00176 template<>
00177 class const_mod<uint64_t, uint64_t(1) << 48>
00178 {
00179   typedef uint64_t IntType;
00180 public:
00181   static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
00182   static IntType mult(IntType a, IntType x) { return mod(a*x); }
00183   static IntType mult_add(IntType a, IntType x, IntType c)
00184     { return mod(a*x+c); }
00185   static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
00186 
00187   // m is not prime, thus invert is not useful
00188 private:                      // don't instantiate
00189   const_mod();
00190 };
00191 
00192 } // namespace random
00193 } // namespace boost
00194 
00195 #endif
00196