algebramix: include/algebramix/crt

00001 
00002 /******************************************************************************
00003 * MODULE     : crt_int.hpp
00004 * DESCRIPTION: CRT subroutines for hardware integers
00005 * COPYRIGHT  : (C) 2008  Joris van der Hoeven and 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_CRT_INT_HPP
00014 #define __MMX_CRT_INT_HPP
00015 #include <numerix/integer.hpp>
00016 #include <numerix/modular_int.hpp>
00017 #include <algebramix/crt_naive.hpp>
00018 
00019 namespace mmx {
00020 
00021 /******************************************************************************
00022 * Variant for hardware integers
00023 ******************************************************************************/
00024 
00025 DEFINE_VARIANT(crt_int,crt_signed<crt_naive>)
00026 DEFINE_VARIANT(crt_unsigned_int,crt_naive)
00027 
00028 #define DECLARE_HELPER(C)                       \
00029   STMPL struct crt_naive_variant_helper<C> {    \
00030     typedef crt_int CV; };
00031 DECLARE_HELPER (signed char)
00032 DECLARE_HELPER (short int)
00033 DECLARE_HELPER (int)
00034 DECLARE_HELPER (long int)
00035 DECLARE_HELPER (long long int)
00036 #undef DECLARE_HELPER
00037 
00038 #define DECLARE_HELPER(C)                       \
00039   STMPL struct crt_naive_variant_helper<C> {    \
00040     typedef crt_unsigned_int CV; };
00041 DECLARE_HELPER (unsigned char)
00042 DECLARE_HELPER (unsigned short int)
00043 DECLARE_HELPER (unsigned int)
00044 DECLARE_HELPER (unsigned long int)
00045 DECLARE_HELPER (unsigned long long int)
00046 #undef DECLARE_HELPER
00047 
00048 /******************************************************************************
00049 * Special transform taking care of casual overflows
00050 ******************************************************************************/
00051 
00052 template<typename V>
00053 struct implementation<crt_transform,V,crt_int> :
00054     public implementation<crt_project,V> {
00055   typedef implementation<crt_project,V> Crt;
00056 
00057   template<typename C, typename M, typename I> static inline void
00058   direct (I* c, const C& a, const M* p, nat n) {
00059     for (nat i= 0; i < n; i++)
00060       c[i]= Crt::mod (a, p[i]); }
00061 
00062   template<typename C, typename M, typename I, typename K> static inline void
00063   inverse (C& a, const I* c, const M* p, const K& P,
00064            const C* q, const I* m, nat n) {
00065     a= 0; I t;
00066     for (nat i= 0; i < n; i++) {
00067       mul_mod (t, m[i], c[i], p[i]);
00068       add_mod (a, t * q[i], P); } }
00069 };
00070 
00071 template<typename V>
00072 struct implementation<crt_transform,V,crt_unsigned_int> :
00073     public implementation<crt_transform,V,crt_int> {};
00074 
00075 /******************************************************************************
00076 * Sequence of prime numbers
00077 ******************************************************************************/
00078 
00079 #define TMPL template<typename M>
00080 
00081 struct prime_sequence_int {
00082 private:
00083   TMPL static inline bool
00084   is_prime (const vector<M>& v, nat n) {
00085     for (nat i= 0; true; i++) {
00086       nat p= (nat) *v[i];
00087       if (n % p == 0) return false;
00088       if (p * p >= n) break;
00089     }
00090     return true; }
00091 
00092 public:
00093   TMPL static bool
00094   extend (vector<M>& v, nat i) {
00095     if (N(v) == 0) v= vec (M(2), M(3));
00096     nat p= *(v[N(v) - 1]) + 2;
00097     while (N(v) < i) {
00098       if (is_prime (v, p)) v << M(p);
00099       p += 2;
00100     }
00101     return true; }
00102 };
00103 
00104 #define DECLARE_HELPER(I)                               \
00105 template<typename V>                                    \
00106  struct coprime_moduli_helper<modulus<I, V> > {         \
00107   typedef prime_sequence_int sequence;                  \
00108 };
00109 DECLARE_HELPER (signed char)
00110 DECLARE_HELPER (unsigned char)
00111 DECLARE_HELPER (short int)
00112 DECLARE_HELPER (unsigned short int)
00113 DECLARE_HELPER (int)
00114 DECLARE_HELPER (unsigned int)
00115 DECLARE_HELPER (long int)
00116 DECLARE_HELPER (unsigned long int)
00117 DECLARE_HELPER (long long int)
00118 DECLARE_HELPER (unsigned long long int)
00119 #undef DECLARE_HELPER
00120 
00121 /******************************************************************************
00122 * Covering
00123 ******************************************************************************/
00124 
00125 template<typename C, typename M, typename W>
00126 struct moduli_unsigned_integer_helper {
00127   // C is an unsigned int type
00128   template<typename V> static bool
00129   covering (vector<M, V>& v, nat s) {
00130     static coprime_moduli_sequence<M,W> seq;
00131     v= vector<M, V> ();
00132     C p (1), q; nat i= 0;
00133     while (bit_size (p) < s + 1) {
00134       if (seq [i] == 0) return false;
00135       v << seq [i];
00136       q = p * C (* seq [i]);
00137       if (quo (q, p) != C (* seq [i])) return false; // means overflow!
00138       p = q;
00139       i++;
00140     }
00141     return true; }
00142 };
00143 
00144 template<typename C, typename M, typename W>
00145 struct moduli_signed_integer_helper {
00146   // C is a signed int type
00147   template<typename V> static bool
00148   covering (vector<M, V>& v, nat s) {
00149     static coprime_moduli_sequence<M,W> seq;
00150     v= vector<M, V> ();
00151     C p (1), q; nat i= 0;
00152     while (bit_size (p) < s + 2) { // s + 2 is needed for sign decoding!
00153       if (seq [i] == 0) return false;
00154       v << seq [i];
00155       q = p * C (* (seq [i]));
00156       if (quo (q, p) != C (* seq [i])) return false; // means overflow!
00157       p = q;
00158       i++;
00159     }
00160     return true; }
00161 };
00162 
00163 #define DECLARE_HELPER(C)                                               \
00164   template<typename M, typename W>                                      \
00165   struct moduli_helper<C,M,W> : moduli_signed_integer_helper<C,M,W> {};
00166 DECLARE_HELPER (signed char)
00167 DECLARE_HELPER (short int)
00168 DECLARE_HELPER (int)
00169 DECLARE_HELPER (long int)
00170 DECLARE_HELPER (long long int)
00171 #undef DECLARE_HELPER
00172 
00173 #define DECLARE_HELPER(C)                                               \
00174   template<typename M, typename W>                                      \
00175   struct moduli_helper<C,M,W> : moduli_unsigned_integer_helper<C,M,W> {};
00176 DECLARE_HELPER (unsigned char)
00177 DECLARE_HELPER (unsigned short int)
00178 DECLARE_HELPER (unsigned int)
00179 DECLARE_HELPER (unsigned long int)
00180 DECLARE_HELPER (unsigned long long int)
00181 #undef DECLARE_HELPER
00182 
00183 /******************************************************************************
00184 * Naive primality test
00185 ******************************************************************************/
00186 
00187 bool is_prime (nat n);
00188 bool is_probable_prime (unsigned long int n);
00189 
00190 /******************************************************************************
00191 * Sequence of prime numbers of bit size s
00192 ******************************************************************************/
00193 
00194 template <nat s>
00195 struct fft_prime_sequence_int {
00196   // prime numbers of bit-size at most s,
00197   // starting with primes suitable to fft.
00198 private:
00199   static nat
00200   bit_reverse (nat i, nat n) {
00201     if (n <= 1) return i;
00202     nat n2= n >> 1;
00203     nat n1= n - n2;
00204     return (bit_reverse (i & ((1 << n2) - 1), n2) << n1) +
00205       bit_reverse (i >> n2, n1);
00206   }
00207 public:
00208   template<typename C, typename V> static bool
00209   extend (vector<modulus<C,V> >& v, nat n) {
00210     // construct prime moduli of s bits exactly,
00211     // starting with those which have a maximal number of
00212     // primitive 2^k-th roots of unity
00213     typedef modulus<C,V> M;
00214     const nat b= V::template maximum_size_helper<C>::value;
00215     ASSERT (s <= b, "bitsize overflow");
00216     ASSERT (b <= 8 * sizeof (nat), "bitsize overflow");
00217     if (s < 2 && n > 0) return false;
00218     const C m= 1 + MMX_SAFE_LEFT_SHIFT_INT (C, 1, (s-1));
00219     const C l= MMX_SAFE_LEFT_SHIFT_INT (C, 1, s);
00220     for (nat i= N(v); i < n; i++) {
00221       C j= i == 0 ? m : (bit_reverse (* v[i-1], s) + 2), k;
00222       for (; j <= l - 1 && j >= m; j += 2)
00223         if (is_prime (k= bit_reverse (j, s))) break;
00224       if (j > l - 1 || j < m) return false;
00225       v << M(k); }
00226     return true; }
00227 };
00228 
00229 /******************************************************************************
00230 * Sequence of probable prime numbers of bit size s
00231 ******************************************************************************/
00232 
00233 template <nat s>
00234 struct probable_prime_sequence_int {
00235   // probable prime numbers of bit-size at least s
00236 
00237   template<typename C, typename V> static bool
00238   extend (vector<modulus<C,V> >& v, nat n) {
00239     // construct probable prime moduli of s bits exactly.
00240     typedef modulus<C,V> M;
00241     const nat b= V::template maximum_size_helper<C>::value;
00242     ASSERT (s <= b, "bitsize overflow");
00243     if (s < 2 && n > 0) return false;
00244     integer a (N(v) == 0 ? (integer (1) << (s-1)) : integer (* v[N(v)-1]));
00245     integer p= probable_next_prime (a);
00246     if (bit_size (p) > b) return false;
00247     v << M(as <C> (p));
00248     return true; }
00249 };
00250 
00251 #undef TMPL
00252 } // namespace mmx
00253 #endif // __MMX_CRT_INT_HPP
include/algebramix/crt_int.hpp
