lll.cpp

Go to the documentation of this file.
00001 /*****************************************************************************
00002     
00003     lll.cpp -- LLL basis reduction.
00004 
00005     This file is a part of the Arageli library.
00006 
00007     Copyright (C) Nikolai Yu. Zolotykh, 1999--2006
00008 
00009 *****************************************************************************/
00010 
00011 #include "config.hpp"
00012 
00013 #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
00014     defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_LLL)
00015 
00016 
00017 #include <cmath>
00018 
00019 #include "vector.hpp"
00020 #include "matrix.hpp"
00021 #include "lll.hpp"
00022 
00023 
00024 namespace Arageli
00025 {
00026 
00027 
00028 #define _ARAGELI_LLL_RED(k, l)                           \
00029 {                                                        \
00030   T v_1_2(1, 2);                                         \
00031   if(v_1_2 < std::abs(Mu(k, l)))                         \
00032   {                                                      \
00033     T q = std::floor(Mu(k, l) + v_1_2);                  \
00034     B.assign_col(k, B.copy_col(k) - B.copy_col(l) * q);  \
00035     H.assign_col(k, H.copy_col(k) - H.copy_col(l) * q);  \
00036     Mu(k, l) = Mu(k, l) - q;                             \
00037     for (index i = 0; i < l; i++)                        \
00038       Mu(k, i) = Mu(k, i) - q * Mu(l, i);                \
00039   }                                                      \
00040 }
00041 
00042 
00043 #define _ARAGELI_LLL_SWAP(k)                                        \
00044 {                                                                   \
00045   B.swap_cols(k, k - 1);                                            \
00046   H.swap_cols(k, k - 1);                                            \
00047   if (k > 1)                                                        \
00048     for (index j = 0; j < k - 1; j++)                               \
00049     {                                                               \
00050       T tmp = Mu(k, j);                                             \
00051       Mu(k, j) = Mu(k - 1, j);                                      \
00052       Mu(k - 1, j) = tmp;                                           \
00053     }                                                               \
00054   T mu = Mu(k, k - 1);                                              \
00055   T b2 = B2[k] + mu * mu * B2[k - 1];                               \
00056   Mu(k, k - 1) = mu * B2[k - 1] / b2;                               \
00057   vector<T> b = Bst.copy_col(k - 1);                                \
00058   Bst.assign_col(k - 1, Bst.copy_col(k) + b * mu);                  \
00059   Bst.assign_col                                                    \
00060     (k, Bst.copy_col(k) * (-Mu(k, k - 1)) + b * (B2[k]/b2));        \
00061   B2[k] = B2[k - 1] * B2[k] / b2;                                   \
00062   B2[k - 1] = b2;                                                   \
00063                                                                     \
00064   for (index i = k + 1; i <= k_max; i++)                            \
00065   {                                                                 \
00066     T t = Mu(i, k);                                                 \
00067     Mu(i, k) = Mu(i, k - 1) - mu * t;                               \
00068     Mu(i, k - 1) = t + Mu(k, k - 1) * Mu(i, k);                     \
00069   }                                                                 \
00070 }
00071 
00072 
00073 template <typename B_type, typename H_type>
00074 bool lll_reduction (B_type& B, H_type& H)
00075 {
00076   typedef typename B_type::difference_type index;
00077   typedef typename B_type::value_type T;
00078     
00079   index m = B.nrows();
00080   index n = B.ncols();
00081   matrix<T, false> Bst(m, n, fromsize);
00082   matrix<T, false> Mu(n);
00083   vector<T, false> B2(n);
00084 
00085   Bst.assign_col(0, B.copy_col(0));
00086   B2[0] = dotprod(B.copy_col(0), B.copy_col(0));
00087   H = H_type(n, eye);
00088   //output_aligned(std::cout << "\nH =\n", H);
00089 
00090   index k_max = 0;
00091 
00092   for(index k = 1; k < n; k++)
00093   {
00094 
00095     // Incremental Gram--Schmidt
00096 
00097     if(k > k_max)
00098     {
00099       k_max = k;
00100       Bst.assign_col(k, B.copy_col(k));
00101       
00102       for(index j = 0; j < k; j++)
00103       {
00104         Mu(k, j) = dotprod(B.copy_col(k), Bst.copy_col(j)) / B2[j];
00105         Bst.assign_col(k, Bst.copy_col(k) - Bst.copy_col(j) * Mu(k, j));
00106         B2[k] = dotprod(Bst.copy_col(k), Bst.copy_col(k));
00107       }
00108       
00109       if(is_null(B2[k]))
00110         return false; // B(i) did not form the basis
00111     }
00112 
00113     // Test LLL-condition
00114 
00115     T v_3_4 = T(3, 4);
00116     
00117     for(;;)
00118     {
00119       //_ARAGELI_LLL_RED(k, k - 1)
00120         {                                                        
00121             T v_1_2(1, 2);                                         
00122             if(v_1_2 < std::abs(Mu(k, k-1)))                         
00123             {                                                      
00124                 T q = std::floor(Mu(k, k-1) + v_1_2);                  
00125                 B.assign_col(k, B.copy_col(k) - B.copy_col(k-1) * q);  
00126                 H.assign_col(k, H.copy_col(k) - H.copy_col(k-1) * q);  
00127                 Mu(k, k-1) = Mu(k, k-1) - q;                             
00128                 for (index i = 0; i < k-1; i++)                        
00129                 Mu(k, i) = Mu(k, i) - q * Mu(k-1, i);                
00130             }                                                      
00131         }
00132 
00133 
00134       if(B2[k] >= (v_3_4 - Mu(k, k - 1) * Mu(k, k - 1)) * B2[k - 1])
00135         break;
00136 
00137       //_ARAGELI_LLL_SWAP(k)
00138       {                                                                   
00139         B.swap_cols(k, k - 1);                                            
00140         H.swap_cols(k, k - 1);                                            
00141         if (k > 1)                                                        
00142             for (index j = 0; j < k - 1; j++)                               
00143             {                                                               
00144             T tmp = Mu(k, j);                                             
00145             Mu(k, j) = Mu(k - 1, j);                                      
00146             Mu(k - 1, j) = tmp;                                           
00147             }                                                               
00148         T mu = Mu(k, k - 1);                                              
00149         T b2 = B2[k] + mu * mu * B2[k - 1];                               
00150         Mu(k, k - 1) = mu * B2[k - 1] / b2;                               
00151         vector<T> b = Bst.copy_col(k - 1);                                
00152         Bst.assign_col(k - 1, Bst.copy_col(k) + b * mu);                  
00153         Bst.assign_col                                                    
00154             (k, Bst.copy_col(k) * (-Mu(k, k - 1)) + b * (B2[k]/b2));        
00155         B2[k] = B2[k - 1] * B2[k] / b2;                                   
00156         B2[k - 1] = b2;                                                   
00157                                                                             
00158         for (index i = k + 1; i <= k_max; i++)                            
00159         {                                                                 
00160             T t = Mu(i, k);                                                 
00161             Mu(i, k) = Mu(i, k - 1) - mu * t;                               
00162             Mu(i, k - 1) = t + Mu(k, k - 1) * Mu(i, k);                     
00163         }                                                                 
00164       }
00165 
00166 
00167       k = (2 > k) ? 1 : k - 1; // MAX(1, k - 1);
00168     }
00169 
00170     for(index l = k - 1; l > 0; l--)
00171       //_ARAGELI_LLL_RED(k, l - 1)
00172     {                                                        
00173         T v_1_2(1, 2);                                         
00174         if(v_1_2 < std::abs(Mu(k, l-1)))                         
00175         {                                                      
00176             T q = std::floor(Mu(k, l-1) + v_1_2);                  
00177             B.assign_col(k, B.copy_col(k) - B.copy_col(l-1) * q);  
00178             H.assign_col(k, H.copy_col(k) - H.copy_col(l-1) * q);  
00179             Mu(k, l-1) = Mu(k, l-1) - q;                             
00180             for (index i = 0; i < l-1; i++)                        
00181             Mu(k, i) = Mu(k, i) - q * Mu(l-1, i);                
00182         }                                                      
00183     }
00184 
00185 
00186   }
00187 
00188   return true;
00189 }
00190 
00191 #undef _ARAGELI_LLL_RED
00192 #undef _ARAGELI_LLL_SWAP
00193 
00194 }
00195 
00196 
00197 #endif // #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE

Generated on Thu Aug 31 17:38:07 2006 for Arageli by  doxygen 1.4.7