algebramix: src/kronecker_int.cpp Source File

00001 
00002 /******************************************************************************
00003 * MODULE     : kronecker_int.cpp
00004 * DESCRIPTION: Breaking big integers into hardware parts and recomposing them
00005 * COPYRIGHT  : (C) 2009  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 #include <numerix/integer.hpp>
00014 namespace mmx {
00015 
00016 /******************************************************************************
00017 * Naive encoding and decoding of integer polynomials using Kronecker's method
00018 ******************************************************************************/
00019 
00020 #if !defined(__GNU_MP__)
00021 template<typename I> static inline void
00022 encode_kronecker_int (integer& dest, const I* src, nat n, xnat bits) {
00023   if (n == 0) dest= 0;
00024   else if (n == 1) dest= src[0];
00025   else {
00026     nat h= n>>1;
00027     integer aux;
00028     encode_kronecker_int (aux , src+h, n-h, bits);
00029     encode_kronecker_int (dest, src  , h  , bits);
00030     dest += (aux << (h*bits));
00031   }
00032 }
00033 
00034 #define IMPLEMENTATION_HELPER(I)                                        \
00035   void                                                                  \
00036   encode_kronecker (integer& dest, const I* src, nat n, xnat bits) {    \
00037     encode_kronecker_int (dest, src, n, bits); }
00038 IMPLEMENTATION_HELPER (unsigned char)
00039 IMPLEMENTATION_HELPER (unsigned short)
00040 IMPLEMENTATION_HELPER (unsigned int)
00041 IMPLEMENTATION_HELPER (unsigned long int)
00042 IMPLEMENTATION_HELPER (unsigned long long int)
00043 IMPLEMENTATION_HELPER (signed char)
00044 IMPLEMENTATION_HELPER (short int)
00045 IMPLEMENTATION_HELPER (int)
00046 IMPLEMENTATION_HELPER (long int)
00047 IMPLEMENTATION_HELPER (long long int)
00048 #undef IMPLEMENTATION_HELPER
00049 
00050 template<typename I> static inline void
00051 decode_kronecker_int (I* dest, nat n, xnat bits, const integer& src) {
00052   static const integer mask= (integer (1) << (8 * sizeof (I))) - 1;
00053   if (n == 0);
00054   else if (n == 1) 
00055     dest[0]= src > 0 ? as<I> (src & mask) : - as<I> ((-src) & mask);
00056   else if (src > 0) {
00057     nat h= n>>1;
00058     integer aux= src >> (h*bits);
00059     if (src[h*bits-1]) aux += 1;
00060     decode_kronecker (dest+h, n-h, bits, aux);
00061     aux= src - (aux << (h*bits));
00062     decode_kronecker (dest, h, bits, aux);
00063   }
00064   else {
00065     integer bis= -src;
00066     nat h= n>>1;
00067     integer aux= bis >> (h*bits);
00068     if (bis[h*bits-1]) aux += 1;
00069     decode_kronecker (dest+h, n-h, bits, -aux);
00070     aux= bis - (aux << (h*bits));
00071     decode_kronecker (dest, h, bits, -aux);
00072   }
00073 }
00074 
00075 template<typename I> static inline void
00076 decode_kronecker_uint (I* dest, nat n, xnat bits, const integer& src) {
00077   if (n == 0);
00078   else if (n == 1) dest[0]= as<I> (src);
00079   else {
00080     nat h= n>>1;
00081     integer aux= src >> (h*bits);
00082     decode_kronecker_int (dest+h, n-h, bits, aux);
00083     aux= src - (aux << (h*bits));
00084     decode_kronecker_int (dest, h, bits, aux);
00085   }
00086 }
00087 
00088 #define IMPLEMENTATION_HELPER(I)                                        \
00089   void                                                                  \
00090   decode_kronecker (I* dest, nat n, xnat bits, const integer& src) {    \
00091     decode_kronecker_uint (dest, n, bits, src); }
00092 IMPLEMENTATION_HELPER (unsigned char)
00093 IMPLEMENTATION_HELPER (unsigned short int)
00094 IMPLEMENTATION_HELPER (unsigned int)
00095 IMPLEMENTATION_HELPER (unsigned long int)
00096 IMPLEMENTATION_HELPER (unsigned long long int)
00097 #undef IMPLEMENTATION_HELPER
00098 
00099 #define IMPLEMENTATION_HELPER(I)                                        \
00100   void                                                                  \
00101   decode_kronecker (I* dest, nat n, xnat bits, const integer& src) {    \
00102     decode_kronecker_int (dest, n, bits, src); }
00103 IMPLEMENTATION_HELPER (signed char)
00104 IMPLEMENTATION_HELPER (short int)
00105 IMPLEMENTATION_HELPER (int)
00106 IMPLEMENTATION_HELPER (long int)
00107 IMPLEMENTATION_HELPER (long long int)
00108 #undef IMPLEMENTATION_HELPER
00109 
00110 #else
00111 /******************************************************************************
00112 * Specific implementation for GMP
00113 ******************************************************************************/
00114 
00115 // from polynomial_gmp.cpp
00116 #define DECLARE_HELPER(I)                                               \
00117   void mpz_encode_kronecker (mpz_t dest, const I* src,                  \
00118                              unsigned long n, unsigned long bits);      \
00119   void mpz_decode_kronecker (I* dest, unsigned long n,                  \
00120                              unsigned long bits, const mpz_t src);
00121 DECLARE_HELPER(unsigned char)
00122 DECLARE_HELPER(unsigned short)
00123 DECLARE_HELPER(unsigned int)
00124 DECLARE_HELPER(unsigned long int)
00125 DECLARE_HELPER(unsigned long long int)
00126 #undef DECLARE_HELPER
00127 
00128 #define IMPLEMENTATION_HELPER(I)                                        \
00129   void encode_kronecker (integer& dest,                                 \
00130                          const I* src, nat n, xnat bits) {              \
00131     typedef unsigned_of_helper<I>::type U;                              \
00132     mpz_encode_kronecker (*dest, (const U*) (const I*) src,             \
00133                           (unsigned long) n,                            \
00134                           (unsigned long) bits); }                      \
00135   void decode_kronecker (I* dest, nat n, xnat bits,                     \
00136                          const integer& src) {                          \
00137     typedef unsigned_of_helper<I>::type U;                              \
00138     mpz_decode_kronecker ((U*) (void*) dest, (unsigned long) n,         \
00139                           (unsigned long) bits, *src); }
00140 IMPLEMENTATION_HELPER(signed char)
00141 IMPLEMENTATION_HELPER(unsigned char)
00142 IMPLEMENTATION_HELPER(short int)
00143 IMPLEMENTATION_HELPER(unsigned short)
00144 IMPLEMENTATION_HELPER(int)
00145 IMPLEMENTATION_HELPER(unsigned int)
00146 IMPLEMENTATION_HELPER(long int)
00147 IMPLEMENTATION_HELPER(unsigned long int)
00148 IMPLEMENTATION_HELPER(long long int)
00149 IMPLEMENTATION_HELPER(unsigned long long int)
00150 #undef IMPLEMENTATION_HELPER
00151 
00152 #endif // __GNU_MP__
00153 
00154 /******************************************************************************
00155 * Multiply two integer polynomials stored in segments using Kronecker's method
00156 ******************************************************************************/
00157 
00158 template<typename I> static inline void
00159 mul_kronecker_int (I* dest,
00160                    const I* src1, nat n1,
00161                    const I* src2, nat n2)
00162 {
00163   if (n1 == 0 && n2 == 0) return;
00164   for (nat i= 0; i < n1 + n2 - 1; i++) dest[i]= 0;
00165   while (n1 > 0 && src1[n1-1] == 0) n1--;
00166   while (n2 > 0 && src2[n2-1] == 0) n2--;
00167   if (n1 == 0 || n2 == 0) return;
00168 
00169   xnat bits= 16 * sizeof (I) + bit_size (min (n1, n2));
00170   integer aux1, aux2;
00171   encode_kronecker (aux1, src1, n1, bits);
00172   encode_kronecker (aux2, src2, n2, bits);
00173   integer aux= aux1 * aux2;
00174   decode_kronecker (dest, n1+n2-1, bits, aux);
00175 }
00176 
00177 #define IMPLEMENTATION_HELPER(I)                                        \
00178   void                                                                  \
00179   mul_kronecker (I* dest, const I* src1, nat n1, const I* src2, nat n2) { \
00180     mul_kronecker_int (dest, src1, n1, src2, n2); }
00181 IMPLEMENTATION_HELPER (signed char)
00182 IMPLEMENTATION_HELPER (unsigned char)
00183 IMPLEMENTATION_HELPER (signed short int)
00184 IMPLEMENTATION_HELPER (unsigned short int)
00185 IMPLEMENTATION_HELPER (int)
00186 IMPLEMENTATION_HELPER (unsigned int)
00187 IMPLEMENTATION_HELPER (long int)
00188 IMPLEMENTATION_HELPER (unsigned long int)
00189 IMPLEMENTATION_HELPER (long long int)
00190 IMPLEMENTATION_HELPER (unsigned long long int)
00191 #undef IMPLEMENTATION_HELPER
00192 
00193 template<typename I> static inline void
00194 square_kronecker_int (I* dest, const I* src, nat n) {
00195   if (n == 0) return;
00196   for (nat i= 0; i < 2 * n - 1; i++) dest[i]= 0;
00197   while (n > 0 && src[n-1] == 0) n--;
00198   if (n == 0) return;
00199   if (n == 1) { dest[0]= square (src[0]); return; }
00200 
00201   xnat bits= 16 * sizeof (I) + bit_size (n);
00202   integer aux1;
00203   encode_kronecker (aux1, src, n, bits);
00204   integer aux= square (aux1);
00205   decode_kronecker (dest, 2*n - 1, bits, aux);
00206 }
00207 
00208 #define IMPLEMENTATION_HELPER(I)                                        \
00209   void                                                                  \
00210   square_kronecker (I* dest, const I* src, nat n) {                     \
00211     square_kronecker_int (dest, src, n); }
00212 IMPLEMENTATION_HELPER (signed char)
00213 IMPLEMENTATION_HELPER (unsigned char)
00214 IMPLEMENTATION_HELPER (signed short int)
00215 IMPLEMENTATION_HELPER (unsigned short int)
00216 IMPLEMENTATION_HELPER (int)
00217 IMPLEMENTATION_HELPER (unsigned int)
00218 IMPLEMENTATION_HELPER (long int)
00219 IMPLEMENTATION_HELPER (unsigned long int)
00220 IMPLEMENTATION_HELPER (long long int)
00221 IMPLEMENTATION_HELPER (unsigned long long int)
00222 #undef IMPLEMENTATION_HELPER
00223 
00224 } // namespace mmx