basix: include/basix/int.hpp Source File

00001 
00002 /******************************************************************************
00003 * MODULE     : int.hpp
00004 * DESCRIPTION: additional operations for hardware integers
00005 * COPYRIGHT  : (C) 2008  Gregoire Lecerf
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX__INT__HPP
00014 #define __MMX__INT__HPP
00015 #include "basix/port.hpp"
00016 #include "basix/basix-config.hpp"
00017 #if BASIX_HAVE_STDINT_H
00018 #include "stdint.h"
00019 #endif
00020 namespace mmx {
00021 
00022 const int int_plus_infinity = (int) (((nat) (-1)) >> 1);
00023 const int int_minus_infinity= (int) ((((nat) (-1)) >> 1) + 1);
00024 
00025 const int int_size= sizeof (int);
00026 const int long_int_size= sizeof (long int);
00027 
00028 #if long_int_size < (int_size << 1)
00029   typedef long long int double_int;
00030 #else
00031   typedef long int double_int;
00032 #endif
00033 
00034 inline bool is_int (const double_int& x) {
00035   return x == ((double_int) (int) x); }
00036 inline int as_int (const double_int& x) {
00037   return (int) x; }
00038 
00039 inline void set_default (int& x) { x= 0; }
00040 inline void set_default (nat& x) { x= 0; }
00041 inline void set_default (long int& x) { x= 0; }
00042 inline void set_default (long unsigned int& x) { x= 0; }
00043 
00044 inline void set_maximal (int& x) {
00045   x= (int) (((nat) (-1)) >> 1); }
00046 inline void set_minimal (int& x) {
00047   x= (int) ((((nat) (-1)) >> 1) + 1); }
00048 inline void set_maximal (nat& x) {
00049   x= (nat) (-1); }
00050 inline void set_minimal (nat& x) {
00051   x= (nat) 0; }
00052 inline void set_maximal (long int& x) {
00053   x= (long int) (((long unsigned int) (-1)) >> 1); }
00054 inline void set_minimal (long int& x) {
00055   x= (long int) ((((long unsigned int) (-1)) >> 1) + 1); }
00056 inline void set_maximal (long unsigned int& x) {
00057   x= (long unsigned int) (-1); }
00058 inline void set_minimal (long unsigned int& x) {
00059   x= (long unsigned int) 0; }
00060 inline void set_maximal (long long int& x) {
00061   x= (long long int) (((long long unsigned int) (-1)) >> 1); }
00062 inline void set_minimal (long long int& x) {
00063   x= (long long int) ((((long long unsigned int) (-1)) >> 1) + 1); }
00064 inline void set_maximal (long long unsigned int& x) {
00065   x= (long long unsigned int) (-1); }
00066 inline void set_minimal (long long unsigned int& x) {
00067   x= (long long unsigned int) 0; }
00068 
00069 inline nat
00070 log_2 (nat p) {
00071   nat i=0;
00072   while (p>1) {
00073     p >>= 1;
00074     i++;
00075   }
00076   return i;
00077 }
00078 
00079 inline nat
00080 next_power_of_two (nat p) {
00081   // return 0 in case of overflow
00082   nat n= 1;
00083   while (n < p && n != 0)
00084     n = n << 1;
00085   return n;
00086 }
00087 
00088 inline bool
00089 is_power_of_two (nat p) {
00090   if (p == 0) return false;
00091   while ((p & 1) == 0) p >>= 1;
00092   return p == 1;
00093 }
00094 
00095 inline nat // inline instead of static in order to avoid warnings
00096 log_3 (nat p) {
00097   nat i=0;
00098   while (p>1) {
00099     p /= 3;
00100     i++;
00101   }
00102   return i;
00103 }
00104 
00105 inline nat
00106 next_power_of_three (nat p) {
00107   // return 0 in case of overflow
00108   nat n= 1;
00109   while (n < p)
00110     n *= 3;
00111   return (n<p) ? 0 : n;
00112 }
00113 
00115 
00116 // These two macro prevent the compiler from complaining
00117 // when n = 8 * sizeof(C)
00118 #define MMX_SAFE_LEFT_SHIFT_INT(C,a,n)                          \
00119   (((n) >= 8*sizeof(C)) ? ((C) 0) : (((C) a) << ((n)/2)) << ((n)-(n)/2))
00120 
00121 #define MMX_SAFE_RIGHT_SHIFT_INT(C,a,n)                         \
00122   (((n) >= 8*sizeof(C)) ? ((C) 0) : (((C) a) >> ((n)/2)) >> ((n)-(n)/2))
00123 
00125 template<nat s>
00126 struct unsigned_int_with_size_at_least_helper {
00127 private:
00128   template<bool b, typename T, typename E>
00129   struct if_helper { typedef T ret; };
00130   template<typename T, typename E>
00131   struct if_helper<false,T,E> { typedef E ret; };
00132 public:
00133   typedef
00134     typename if_helper<sizeof(unsigned char)>=s, unsigned char,
00135     typename if_helper<sizeof(short unsigned int)>=s, short unsigned int,
00136     typename if_helper<sizeof(unsigned int)>=s, unsigned int,
00137     typename if_helper<sizeof(long unsigned int)>=s, long unsigned int,
00138     typename if_helper<sizeof(long long unsigned int)>=s,
00139                               long long unsigned int,
00140       void >::ret >::ret >::ret >::ret >::ret type ; };
00141 
00143 template<typename C>
00144 struct unsigned_int_with_double_size_helper {
00145   typedef typename
00146   unsigned_int_with_size_at_least_helper<2*sizeof(C)>::type
00147   type; };
00148 
00150 template<typename C>
00151 struct signed_of_helper {
00152   typedef void type; };
00153 
00154 template<>
00155 struct signed_of_helper<double> {
00156   typedef double type; };
00157 
00158 template<>
00159 struct signed_of_helper<char> {
00160   typedef signed char type; };
00161 
00162 template<>
00163 struct signed_of_helper<signed char> {
00164   typedef signed char type; };
00165 
00166 template<>
00167 struct signed_of_helper<short int> {
00168   typedef short int type; };
00169 
00170 template<>
00171 struct signed_of_helper<int> {
00172   typedef int type; };
00173 
00174 template<>
00175 struct signed_of_helper<long int> {
00176   typedef long int type; };
00177 
00178 template<>
00179 struct signed_of_helper<long long int> {
00180   typedef long long int type; };
00181 
00182 template<>
00183 struct signed_of_helper<unsigned char> {
00184   typedef signed char type; };
00185 
00186 template<>
00187 struct signed_of_helper<short unsigned int> {
00188   typedef short int type; };
00189 
00190 template<>
00191 struct signed_of_helper<unsigned int> {
00192   typedef int type; };
00193 
00194 template<>
00195 struct signed_of_helper<long unsigned int> {
00196   typedef long int type; };
00197 
00198 template<>
00199 struct signed_of_helper<long long unsigned int> {
00200   typedef long long int type; };
00201 
00203 template<typename C>
00204 struct unsigned_of_helper {
00205   typedef void type; };
00206 
00207 template<>
00208 struct unsigned_of_helper<char> {
00209   typedef unsigned char type; };
00210 
00211 template<>
00212 struct unsigned_of_helper<signed char> {
00213   typedef unsigned char type; };
00214 
00215 template<>
00216 struct unsigned_of_helper<short int> {
00217   typedef unsigned short int type; };
00218 
00219 template<>
00220 struct unsigned_of_helper<int> {
00221   typedef unsigned int type; };
00222 
00223 template<>
00224 struct unsigned_of_helper<long int> {
00225   typedef long unsigned int type; };
00226 
00227 template<>
00228 struct unsigned_of_helper<long long int> {
00229   typedef long long unsigned int type; };
00230 
00231 template<>
00232 struct unsigned_of_helper<unsigned char> {
00233   typedef unsigned char type; };
00234 
00235 template<>
00236 struct unsigned_of_helper<short unsigned int> {
00237   typedef short unsigned int type; };
00238 
00239 template<>
00240 struct unsigned_of_helper<unsigned int> {
00241   typedef unsigned int type; };
00242 
00243 template<>
00244 struct unsigned_of_helper<long unsigned int> {
00245   typedef long unsigned int type; };
00246 
00247 template<>
00248 struct unsigned_of_helper<long long unsigned int> {
00249   typedef long long unsigned int type; };
00250 
00252 template<typename C>
00253 struct is_signed_helper {
00254   static const bool value = false; };
00255 
00256 template<>
00257 struct is_signed_helper<signed char> {
00258   static const bool value = true; };
00259 
00260 template<>
00261 struct is_signed_helper<char> {
00262   static const bool value = true; };
00263 
00264 template<>
00265 struct is_signed_helper<short int> {
00266   static const bool value = true; };
00267 
00268 template<>
00269 struct is_signed_helper<int> {
00270   static const bool value = true; };
00271 
00272 template<>
00273 struct is_signed_helper<long int> {
00274   static const bool value = true; };
00275 
00276 template<>
00277 struct is_signed_helper<long long int> {
00278   static const bool value = true; };
00279 
00281 template<typename uC, nat s, uC x>
00282 struct int_bitsize_helper_rec {
00283   static const nat s2  = s >> 1;
00284   static const uC mask = (uC) MMX_SAFE_LEFT_SHIFT_INT(uC,-1,s2);
00285   static const uC size = (x & mask) ?
00286     (s2 + int_bitsize_helper_rec< uC, s2, (x >> s2)>::size) :
00287     int_bitsize_helper_rec< uC, s2, x >::size;
00288 };
00289   
00290 template<typename uC, uC x>
00291 struct int_bitsize_helper_rec <uC, 1, x> {
00292   static const uC size = (x == 0) ? 0 : 1; };
00293 
00294 template<typename uC, uC x>
00295 struct int_bitsize_helper_rec <uC, 0, x> {
00296   static const uC size = 0; };
00297 
00298 template<typename C, C p>
00299 struct int_bitsize_helper {
00300   typedef typename unsigned_of_helper<C>::type uC;
00301   static const uC up = (p < 0) ? -p : p;
00302 
00303   static const nat value =
00304     int_bitsize_helper_rec < uC, 8*sizeof(C), up >::size;
00305 };
00306 
00308 template<typename C> static inline nat
00309 bit_size (const C& p) {
00310   typedef typename unsigned_of_helper<C>::type uC;
00311   if (p == 0) return 0;
00312   uC up = abs (p);
00313   nat s = 0;
00314   nat k = 4 * sizeof(C);
00315   uC mask = ((uC) -1) << k;
00316   while (k != 0) {
00317     if (up & mask) { up >>= k; s += k; }
00318     k >>= 1;
00319     mask >>= k;
00320   }
00321   return up == 0 ? s : s + 1;
00322 }
00323 
00324 template<typename C>
00325 class long_int_mul_op {
00326   typedef typename unsigned_of_helper<C>::type uC;
00327   static const nat n2 = 4 * sizeof (uC);
00328   static const nat n  = 2 * n2;
00329   static const uC lo_mask = ((uC) -1) >> n2;
00330   static const uC carry = ((uC) 1) << n2;
00331   // Let B = (1 << n2),  R = B^2.
00332 
00333   static inline uC
00334   lo (const uC& s) { return s & lo_mask; }
00335     
00336   static inline uC
00337   hi (const uC& s) { return s >> n2; }
00338 
00339 public:
00341   static void
00342   op (uC& dest1, uC& dest0, const uC& a, const uC& b) {
00343     uC t0, t1;
00344     uC lo_a = lo (a);
00345     uC hi_a = hi (a);
00346     uC lo_b = lo (b);
00347     uC hi_b = hi (b);
00348     dest0 = lo_a * lo_b;
00349     dest1 = hi_a * hi_b;
00350     t0 = hi_a * lo_b;
00351     t1 = t0 + lo_a * hi_b;
00352     if (t1 < t0) dest1 += carry;
00353     t0 = lo (t1) << n2;
00354     dest0 += t0;
00355     dest1 += hi (t1);
00356     if (dest0 < t0) dest1++;
00357   }
00358 };
00359 
00360 template<typename C>
00361 struct long_int_rshift_op {
00362   typedef typename unsigned_of_helper<C>::type uC;
00363   static const nat n2 = 4 * sizeof (uC);
00364   static const nat n  = 2 * n2;
00365   // Let B = (1 << n2),  R = B^2.
00366 
00368   static inline void
00369   op (uC& dest1, uC& dest0, nat s) {
00370     if (s == 0) return;
00371     if (s < n) {
00372       dest0 = (dest0 >> s) | (dest1 << (n - s));
00373       dest1 = dest1 >> s;
00374     }
00375     else {
00376       dest0 = dest1 >> (s - n);
00377       dest1 = 0;
00378     }
00379   }
00380     
00383   static inline bool
00384   op_b (uC& dest1, uC& dest0, nat s) {
00385     bool ans;
00386     if (s == 0) return false;
00387     if (s < n) {
00388       ans = (dest0 & (((uC) -1) >> (n - s))) != 0;
00389       dest0 = (dest0 >> s) | (dest1 << (n - s));
00390       dest1 = dest1 >> s;
00391     }
00392     else {
00393       ans = dest0 != 0 || (dest1 & (((uC) -1) >> (2*n - s))) != 0;
00394       dest0 = dest1 >> (s - n);
00395       dest1 = 0;
00396     }
00397     return ans;
00398   }
00399 };
00400 
00401 template<typename C>
00402 struct long_int_lshift_op {
00403   typedef typename unsigned_of_helper<C>::type uC;
00404   static const nat n2 = 4 * sizeof (uC);
00405   static const nat n  = 2 * n2;
00406   // Let B = (1 << n2),  R = B^2.
00407 
00409   static inline void
00410   op (uC& dest1, uC& dest0, nat s) {
00411     if (s == 0) return;
00412     if (s < n) {
00413       dest1 = (dest1 << s) | (dest0 >> (n - s));
00414       dest0 = dest0 << s;
00415     }
00416     else {
00417       dest1 = dest0 << (s - n);
00418       dest0 = 0;
00419     }
00420   }
00421 };
00422 
00423 template<typename C>
00424 struct long_int_sub_op {
00425   typedef typename unsigned_of_helper<C>::type uC;
00426 
00428   static inline void
00429   op (uC& dest1, uC& dest0, const uC& s1, const uC& s0) {
00430     if (dest0 >= s0) {
00431       dest0 -= s0;
00432       dest1 -= s1;
00433     }
00434     else {
00435       dest0 -= s0;
00436       dest1 -= s1 + 1;
00437     }
00438   }
00439 };
00440 
00441 template<typename C>
00442 struct long_int_ge_op {
00443   typedef typename unsigned_of_helper<C>::type uC;
00444 
00446   static inline bool
00447   op (const uC& s1, const uC& s0, const uC& t1, const uC& t0) {
00448     if (s1 > t1) return true;
00449     if (s1 < t1) return false;
00450     return s0 >= t0;
00451   }
00452 };
00453 
00454 /******************************************************************************
00455 * Square root approximation
00456 ******************************************************************************/
00457 
00458 template<typename I>
00459 struct int_floor_sqrt_helper {
00460   static I floor_sqrt (const I& x) {
00461     ASSERT (x >= 0, "wrong sign");
00462     I y= x, r= 0, b= ((I) 1) << (8 * sizeof (I) - 2);
00463     while (b > y) b >>= 2;
00464     while (b != 0) {
00465         if (y >= r + b) { y -= r + b; r = (r >> 1) + b; }
00466         else r >>= 1;
00467         b >>= 2;
00468     }
00469     return r; }
00470 };
00471 
00472 #define INT_FLOOR_SQRT_DECLARE(I)                                       \
00473   inline I floor_sqrt (const I& x) {                                    \
00474     return int_floor_sqrt_helper<I>::floor_sqrt (x); }
00475 
00476 INT_FLOOR_SQRT_DECLARE(signed char)
00477 INT_FLOOR_SQRT_DECLARE(short int)
00478 INT_FLOOR_SQRT_DECLARE(int)
00479 INT_FLOOR_SQRT_DECLARE(long int)
00480 INT_FLOOR_SQRT_DECLARE(long long int)
00481 #undef INT_FLOOR_SQRT_DECLARE
00482 
00483 template<typename I>
00484 struct unsigned_int_floor_sqrt_helper {
00485   static I floor_sqrt (const I& x) {
00486     I y= x, r= 0, b= ((I) 1) << (8 * sizeof (I) - 2);
00487     while (b > y) b >>= 2;
00488     while (b != 0) {
00489         if (y >= r + b) { y -= r + b; r = (r >> 1) + b; }
00490         else r >>= 1;
00491         b >>= 2;
00492     }
00493     return r; }
00494 };
00495 
00496 #define UNSIGNED_INT_FLOOR_SQRT_DECLARE(I)                              \
00497   inline I floor_sqrt (const I& x) {                                    \
00498     return unsigned_int_floor_sqrt_helper<I>::floor_sqrt (x); }
00499 
00500 UNSIGNED_INT_FLOOR_SQRT_DECLARE(unsigned char)
00501 UNSIGNED_INT_FLOOR_SQRT_DECLARE(unsigned short int)
00502 UNSIGNED_INT_FLOOR_SQRT_DECLARE(unsigned int)
00503 UNSIGNED_INT_FLOOR_SQRT_DECLARE(unsigned long int)
00504 UNSIGNED_INT_FLOOR_SQRT_DECLARE(unsigned long long int)
00505 #undef UNSIGNED_INT_FLOOR_SQRT_DECLARE
00506 
00507 /******************************************************************************
00508 * Division
00509 ******************************************************************************/
00510 
00511 #define INT_SIGNED_DIV_DECLARE(I)                                       \
00512   inline bool divides (const I& n, const I& m) { return (m % n) == 0; } \
00513   inline I rem (const I& n, const I& m, I& q) {                         \
00514     I _q= 0;                                                            \
00515     if (m == 0) { q= 0; return n; }                                     \
00516     if (m < 0) {                                                        \
00517       I r= rem (n, -m, _q);                                             \
00518       q= -_q;                                                           \
00519       if (r == 0) return r;                                             \
00520       q--; return r + m;                                                \
00521     }                                                                   \
00522     if (n < 0) {                                                        \
00523       _q= (-n) / m; _q= -_q;                                            \
00524       I r= n - _q * m; q= _q;                                           \
00525       if (r == 0) return r;                                             \
00526       q--; return r + m;                                                \
00527     }                                                                   \
00528     _q= n / m; I r= n - _q * m; q= _q; return r; }                      \
00529   inline I quo (const I& n, const I& m) {                               \
00530     if (m == 0) return 0;                                               \
00531     if (m < 0) {                                                        \
00532       I q, r= rem (n, -m, q);                                           \
00533       q= -q;                                                            \
00534       if (r != 0) q--;                                                  \
00535       return q;                                                         \
00536     }                                                                   \
00537     if (n < 0) {                                                        \
00538       I q= n / m;                                                       \
00539       I r= n - q * m;                                                   \
00540       if (r != 0) q--;                                                  \
00541       return q;                                                         \
00542     }                                                                   \
00543     return n / m; }                                                     \
00544   inline I rem (const I& n, const I& m) {                               \
00545     if (m == 0) { return n; }                                           \
00546     if (m < 0) {                                                        \
00547       I r= rem (n, -m);                                                 \
00548       if (r == 0) return r;                                             \
00549       return r + m;                                                     \
00550     }                                                                   \
00551     if (n < 0) {                                                        \
00552       I r= n % m;                                                       \
00553       if (r == 0) return r;                                             \
00554       return r + m;                                                     \
00555     }                                                                   \
00556     return n % m; }
00557 INT_SIGNED_DIV_DECLARE(signed char)
00558 INT_SIGNED_DIV_DECLARE(short int)
00559 INT_SIGNED_DIV_DECLARE(int)
00560 INT_SIGNED_DIV_DECLARE(long int)
00561 INT_SIGNED_DIV_DECLARE(long long int)
00562 #undef INT_SIGNED_DIV_DECLARE
00563 
00564 #define INT_UNSIGNED_DIV_DECLARE(I)                                     \
00565   inline bool divides (const I& n, const I& m) { return (m % n) == 0; } \
00566   inline I rem (const I& n, const I& m) { return m == 0 ? n : n % m; }  \
00567   inline I quo (const I& n, const I& m) { return m == 0 ? 0 : n / m; }  \
00568   inline I rem (const I& n, const I& m, I& q) {                         \
00569     I _q= quo (n, m); I _r= n - _q * m; q= _q; return _r; }
00570 INT_UNSIGNED_DIV_DECLARE(unsigned char)
00571 INT_UNSIGNED_DIV_DECLARE(unsigned short int)
00572 INT_UNSIGNED_DIV_DECLARE(unsigned int)
00573 INT_UNSIGNED_DIV_DECLARE(unsigned long int)
00574 INT_UNSIGNED_DIV_DECLARE(unsigned long long int)
00575 #undef INT_UNSIGNED_DIV_DECLARE
00576 
00577 /******************************************************************************
00578 * Gcd
00579 ******************************************************************************/
00580 
00581 struct int_gcd_helper {
00582   template<typename I>
00583   static I gcd (const I& a, const I& b) {
00584     typedef typename unsigned_of_helper<I>::type U;
00585     I r0 = a, r1 = b, q, t;
00586     if ((r0 == 0) && (r1 != 0)) {
00587       q = (((U) r0) - ((U) r1)) / ((U) r1) + 1;
00588       t = r0 - q * r1;
00589       r0 = r1;
00590       r1 = t;
00591     }
00592     while (r1 != 0) {
00593       q = r0 / r1;
00594       t = r0 - q * r1;
00595       r0 = r1;
00596       r1 = t;
00597     }
00598     return r0;
00599   }
00600 
00601   template<typename I>
00602   static I gcd (const I& a, const I& b, I& co_a) {
00603     typedef typename unsigned_of_helper<I>::type U;
00604     I r0 = a, r1 = b, co_a0 = 1, co_a1 = 0, q, t;
00605     if ((r0 == 0) && (r1 != 0)) {
00606       q = (((U) r0) - ((U) r1)) / ((U) r1) + 1;
00607       t = r0 - q * r1;
00608       r0 = r1;
00609       r1 = t;
00610       t = co_a1;
00611       co_a1 = co_a0 - q * co_a1;
00612       co_a0 = t;
00613     }
00614     while (r1 != 0) {
00615       q = r0 / r1;
00616       t = r0 - q * r1;
00617       r0 = r1;
00618       r1 = t;
00619       t = co_a1;
00620       co_a1 = co_a0 - q * co_a1;
00621       co_a0 = t;
00622     }
00623     co_a = co_a0;
00624     return r0;
00625   }
00626   
00627   template<typename I>
00628   static I gcd (const I& a, const I& b, I& co_a, I& co_b) {
00629     typedef typename unsigned_of_helper<I>::type U;
00630     I r0 = a, r1 = b, co_a0 = 1, co_a1 = 0, co_b0 = 0, co_b1 = 1, q, t;
00631     if ((r0 == 0) && (r1 != 0)) {
00632       q = (((U) r0) - ((U) r1)) / ((U) r1) + 1;
00633       t = r0 - q * r1;
00634       r0 = r1;
00635       r1 = t;
00636       t = co_a1;
00637       co_a1 = co_a0 - q * co_a1;
00638       co_a0 = t;
00639       t = co_b1;
00640       co_b1 = co_b0 - q * co_b1;
00641       co_b0 = t;
00642     }
00643     while (r1 != 0) {
00644       q = r0 / r1;
00645       t = r0 - q * r1;
00646       r0 = r1;
00647       r1 = t;
00648       t = co_a1;
00649       co_a1 = co_a0 - q * co_a1;
00650       co_a0 = t;
00651       t = co_b1;
00652       co_b1 = co_b0 - q * co_b1;
00653       co_b0 = t;
00654     }
00655     co_a = co_a0;
00656     co_b = co_b0;
00657     return r0;
00658   }
00659 };
00660 
00661 struct unsigned_int_gcd_helper {
00662   template<typename I>
00663   static I gcd (const I& a, const I& b) {
00664     I r0=a, r1=b, q, t; 
00665     if ((r0 == 0) && (r1 != 0)) {
00666       q = (r0-r1) / r1 + 1;
00667       t = r0 - q * r1;
00668       r0 = r1;
00669       r1 = t;
00670     }
00671     while (r1 != 0) {
00672       q = r0 / r1;
00673       t = r0 - q * r1;
00674       r0 = r1;
00675       r1 = t;
00676     }
00677     return r0;
00678   }
00679 };
00680 
00681 #define INT_GCD_DECLARE(I)                                              \
00682   inline I gcd (const I a, const I b) {                                 \
00683     return int_gcd_helper::gcd (a, b); }                                \
00684                                                                         \
00685   inline I gcd (const I a, const I b, I& co_a) {                        \
00686     return int_gcd_helper::gcd (a, b, co_a); }                          \
00687                                                                         \
00688   inline I gcd (const I a, const I b, I& co_a, I& co_b) {               \
00689     return int_gcd_helper::gcd (a, b, co_a, co_b); }
00690   
00691 INT_GCD_DECLARE(signed char)
00692 INT_GCD_DECLARE(short int)
00693 INT_GCD_DECLARE(int)
00694 INT_GCD_DECLARE(long int)
00695 INT_GCD_DECLARE(long long int)
00696 #undef INT_GCD_DECLARE
00697 
00698 #define UINT_GCD_DECLARE(I)                                             \
00699   inline I gcd (const I a, const I b) {                                 \
00700     return unsigned_int_gcd_helper::gcd (a, b); }
00701 
00702 UINT_GCD_DECLARE(unsigned char)
00703 UINT_GCD_DECLARE(unsigned short int)
00704 UINT_GCD_DECLARE(unsigned int)
00705 UINT_GCD_DECLARE(unsigned long int)
00706 UINT_GCD_DECLARE(unsigned long long int)
00707 #undef UINT_GCD_DECLARE
00708 
00709 /******************************************************************************
00710 * Rational reconstruction
00711 ******************************************************************************/
00712 
00713 template<typename I>
00714 struct int_reconstruct_helper {
00715   static bool reconstruct (I& n, I& d, const I& u, const I& m,
00716                            const I& N, const I& D) {
00717   /* Return in r = n / d such that abs (n) <= N, 0 < d <= D, gcd (n, d) = 1,
00718      and n / d = u mod m. If such a rational exists then the function returns
00719      true and false otherwise. */
00720   ASSERT (u >= 0, "invalid input");
00721   ASSERT (m > u, "invalid input");
00722   ASSERT (N > 0, "invalid input");
00723   ASSERT (D > 0, "invalid input");
00724   I aux= N * D;
00725   ASSERT (aux / D == N && (aux << 1) / 2 == aux && (aux << 1) < m,
00726           "invalid input");
00727   I r0= m, t0= 0;
00728   I r1= u, t1= 1;
00729   I q;
00730   while (r1 > N) {
00731     I aux= rem (r0, r1, q);
00732     r0= r1; r1= aux;
00733     aux= t0;
00734     t0= t1; t1= aux - q * t1;
00735   }
00736   n= r1; d= t1;
00737   if (d < 0) { n= -n; d= -d; }
00738   if (d > D) return false;
00739   return abs (gcd (n, d)) == 1; }
00740 };
00741 
00742 #define INT_RECONSTRUCT_DECLARE(I)                                      \
00743   inline bool reconstruct (I& n, I& d, const I& u, const I& m,          \
00744                            const I& N, const I& D) {                    \
00745     return int_reconstruct_helper<I>::reconstruct (n, d, u, m, N, D); } \
00746   inline bool reconstruct (I& n, I& d, const I& u, const I& m) {        \
00747     I N= floor_sqrt (m >> 1);                                           \
00748     return int_reconstruct_helper<I>::reconstruct (n, d, u, m, N, N); }
00749 
00750 INT_RECONSTRUCT_DECLARE(signed char)
00751 INT_RECONSTRUCT_DECLARE(signed short int)
00752 INT_RECONSTRUCT_DECLARE(signed int)
00753 INT_RECONSTRUCT_DECLARE(signed long int)
00754 INT_RECONSTRUCT_DECLARE(signed long long int)
00755 #undef INT_RECONSTRUCT_DECLARE
00756 
00757 #define I typename signed_of_helper<U>::type
00758 template<typename U>
00759 struct unsigned_int_reconstruct_helper {
00760   static bool reconstruct (I& n, I& d, const U& u, const U& m,
00761                            const U& N, const U& D) {
00762   /* Return in r = n / d such that abs (n) <= N, 0 < d <= D, gcd (n, d) = 1,
00763      and n / d = u mod m. If such a rational exists then the function returns
00764      true and false otherwise. */
00765   ASSERT (m > u, "invalid input");
00766   ASSERT (N > 0, "invalid input");
00767   ASSERT (D > 0, "invalid input");
00768   U aux= N * D;
00769   ASSERT (aux / D == N && (aux << 1) / 2 == aux && (aux << 1) < m,
00770           "invalid input");
00771   U r0= m, t0= 0;
00772   U r1= u, t1= 1;
00773   U q;
00774   while (r1 > N) {
00775     U aux= rem (r0, r1, q);
00776     r0= r1; r1= aux;
00777     aux= t0;
00778     t0= t1; t1= aux - q * t1;
00779   }
00780   n= (I) r1; d= (I) t1;
00781   if (d < 0) { n= -n; d= -d; }
00782   if ((U) d > D) return false;
00783   return abs (gcd (n, d)) == 1; }
00784 };
00785 #undef I
00786 
00787 #define UNSIGNED_INT_RECONSTRUCT_DECLARE(U)                             \
00788   inline bool reconstruct (signed_of_helper<U>::type& n,                \
00789                            signed_of_helper<U>::type& d,                \
00790                            const U& u, const U& m,                      \
00791                            const U& N, const U& D) {                    \
00792     return unsigned_int_reconstruct_helper<U>                           \
00793       ::reconstruct (n, d, u, m, N, D); }                               \
00794   inline bool reconstruct (signed_of_helper<U>::type& n,                \
00795                            signed_of_helper<U>::type& d,                \
00796                            const U& u, const U& m) {                    \
00797     U N= floor_sqrt (m >> 1);                                           \
00798     return unsigned_int_reconstruct_helper<U>::                         \
00799       reconstruct (n, d, u, m, N, N); }
00800 
00801 UNSIGNED_INT_RECONSTRUCT_DECLARE(unsigned char)
00802 UNSIGNED_INT_RECONSTRUCT_DECLARE(unsigned short int)
00803 UNSIGNED_INT_RECONSTRUCT_DECLARE(unsigned int)
00804 UNSIGNED_INT_RECONSTRUCT_DECLARE(unsigned long int)
00805 UNSIGNED_INT_RECONSTRUCT_DECLARE(unsigned long long int)
00806 #undef UNSIGNED_INT_RECONSTRUCT_DECLARE
00807 
00808 /******************************************************************************
00809 * stdint helpers
00810 ******************************************************************************/
00811 
00812 #if BASIX_HAVE_STDINT_H
00813 template<int s>
00814 struct unsigned_stdint_of_size_at_least_helper {
00815 private:
00816   template<bool b, typename T, typename E>
00817   struct if_helper { typedef T ret; };
00818   template<typename T, typename E>
00819   struct if_helper<false,T,E> { typedef E ret; };
00820 public:
00821   typedef
00822     typename if_helper<s <= 8, uint8_t,
00823     typename if_helper<s <= 16, uint16_t,
00824     typename if_helper<s <= 32, uint32_t,
00825     typename if_helper<s <= 64, uint64_t,
00826       void >::ret >::ret >::ret >::ret type ; };
00827 
00828 template<int s>
00829 struct signed_stdint_of_size_at_least_helper {
00830 private:
00831   template<bool b, typename T, typename E>
00832   struct if_helper { typedef T ret; };
00833   template<typename T, typename E>
00834   struct if_helper<false,T,E> { typedef E ret; };
00835 public:
00836   typedef
00837     typename if_helper<s <= 8, int8_t,
00838     typename if_helper<s <= 16, int16_t,
00839     typename if_helper<s <= 32, int32_t,
00840     typename if_helper<s <= 64, int64_t,
00841       void >::ret >::ret >::ret >::ret type ; };
00842 
00843 template<typename C>
00844 struct stdint_of_helper {
00845 private:
00846   template<bool b, typename T, typename E>
00847   struct if_helper { typedef T ret; };
00848   template<typename T, typename E>
00849   struct if_helper<false,T,E> { typedef E ret; };
00850 public:
00851   typedef typename if_helper<is_signed_helper<C>::value,
00852     typename signed_stdint_of_size_at_least_helper<8*sizeof(C)>::type,
00853     typename unsigned_stdint_of_size_at_least_helper<8*sizeof(C)>::type>::ret
00854   type;
00855 };
00856 
00857 #endif // BASIX_HAVE_STDINT_H
00858 
00859 }// namespace mmx
00860 #endif //__MMX__INT__HPP