lll.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    lll.cpp -- LLL basis reduction.

    This file is a part of the Arageli library.

    Copyright (C) Nikolai Yu. Zolotykh, 1999--2006

*****************************************************************************/

#include "config.hpp"

#if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
    defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_LLL)


#include <cmath>

#include "vector.hpp"
#include "matrix.hpp"
#include "lll.hpp"


namespace Arageli
{


#define _ARAGELI_LLL_RED(k, l)                           \
{                                                        \
  T v_1_2(1, 2);                                         \
  if(v_1_2 < std::abs(Mu(k, l)))                         \
  {                                                      \
    T q = std::floor(Mu(k, l) + v_1_2);                  \
    B.assign_col(k, B.copy_col(k) - B.copy_col(l) * q);  \
    H.assign_col(k, H.copy_col(k) - H.copy_col(l) * q);  \
    Mu(k, l) = Mu(k, l) - q;                             \
    for (index i = 0; i < l; i++)                        \
      Mu(k, i) = Mu(k, i) - q * Mu(l, i);                \
  }                                                      \
}


#define _ARAGELI_LLL_SWAP(k)                                        \
{                                                                   \
  B.swap_cols(k, k - 1);                                            \
  H.swap_cols(k, k - 1);                                            \
  if (k > 1)                                                        \
    for (index j = 0; j < k - 1; j++)                               \
    {                                                               \
      T tmp = Mu(k, j);                                             \
      Mu(k, j) = Mu(k - 1, j);                                      \
      Mu(k - 1, j) = tmp;                                           \
    }                                                               \
  T mu = Mu(k, k - 1);                                              \
  T b2 = B2[k] + mu * mu * B2[k - 1];                               \
  Mu(k, k - 1) = mu * B2[k - 1] / b2;                               \
  vector<T> b = Bst.copy_col(k - 1);                                \
  Bst.assign_col(k - 1, Bst.copy_col(k) + b * mu);                  \
  Bst.assign_col                                                    \
    (k, Bst.copy_col(k) * (-Mu(k, k - 1)) + b * (B2[k]/b2));        \
  B2[k] = B2[k - 1] * B2[k] / b2;                                   \
  B2[k - 1] = b2;                                                   \
                                                                    \
  for (index i = k + 1; i <= k_max; i++)                            \
  {                                                                 \
    T t = Mu(i, k);                                                 \
    Mu(i, k) = Mu(i, k - 1) - mu * t;                               \
    Mu(i, k - 1) = t + Mu(k, k - 1) * Mu(i, k);                     \
  }                                                                 \
}


template <typename B_type, typename H_type>
bool lll_reduction (B_type& B, H_type& H)
{
  typedef typename B_type::difference_type index;
  typedef typename B_type::value_type T;
    
  index m = B.nrows();
  index n = B.ncols();
  matrix<T, false> Bst(m, n, fromsize);
  matrix<T, false> Mu(n);
  vector<T, false> B2(n);

  Bst.assign_col(0, B.copy_col(0));
  B2[0] = dotprod(B.copy_col(0), B.copy_col(0));
  H = H_type(n, eye);
  //output_aligned(std::cout << "\nH =\n", H);

  index k_max = 0;

  for(index k = 1; k < n; k++)
  {

    // Incremental Gram--Schmidt

    if(k > k_max)
    {
      k_max = k;
      Bst.assign_col(k, B.copy_col(k));
      
      for(index j = 0; j < k; j++)
      {
        Mu(k, j) = dotprod(B.copy_col(k), Bst.copy_col(j)) / B2[j];
        Bst.assign_col(k, Bst.copy_col(k) - Bst.copy_col(j) * Mu(k, j));
        B2[k] = dotprod(Bst.copy_col(k), Bst.copy_col(k));
      }
      
      if(is_null(B2[k]))
        return false; // B(i) did not form the basis
    }

    // Test LLL-condition

    T v_3_4 = T(3, 4);
    
    for(;;)
    {
      //_ARAGELI_LLL_RED(k, k - 1)
        {                                                        
            T v_1_2(1, 2);                                         
            if(v_1_2 < std::abs(Mu(k, k-1)))                         
            {                                                      
                T q = std::floor(Mu(k, k-1) + v_1_2);                  
                B.assign_col(k, B.copy_col(k) - B.copy_col(k-1) * q);  
                H.assign_col(k, H.copy_col(k) - H.copy_col(k-1) * q);  
                Mu(k, k-1) = Mu(k, k-1) - q;                             
                for (index i = 0; i < k-1; i++)                        
                Mu(k, i) = Mu(k, i) - q * Mu(k-1, i);                
            }                                                      
        }


      if(B2[k] >= (v_3_4 - Mu(k, k - 1) * Mu(k, k - 1)) * B2[k - 1])
        break;

      //_ARAGELI_LLL_SWAP(k)
      {                                                                   
        B.swap_cols(k, k - 1);                                            
        H.swap_cols(k, k - 1);                                            
        if (k > 1)                                                        
            for (index j = 0; j < k - 1; j++)                               
            {                                                               
            T tmp = Mu(k, j);                                             
            Mu(k, j) = Mu(k - 1, j);                                      
            Mu(k - 1, j) = tmp;                                           
            }                                                               
        T mu = Mu(k, k - 1);                                              
        T b2 = B2[k] + mu * mu * B2[k - 1];                               
        Mu(k, k - 1) = mu * B2[k - 1] / b2;                               
        vector<T> b = Bst.copy_col(k - 1);                                
        Bst.assign_col(k - 1, Bst.copy_col(k) + b * mu);                  
        Bst.assign_col                                                    
            (k, Bst.copy_col(k) * (-Mu(k, k - 1)) + b * (B2[k]/b2));        
        B2[k] = B2[k - 1] * B2[k] / b2;                                   
        B2[k - 1] = b2;                                                   
                                                                            
        for (index i = k + 1; i <= k_max; i++)                            
        {                                                                 
            T t = Mu(i, k);                                                 
            Mu(i, k) = Mu(i, k - 1) - mu * t;                               
            Mu(i, k - 1) = t + Mu(k, k - 1) * Mu(i, k);                     
        }                                                                 
      }


      k = (2 > k) ? 1 : k - 1; // MAX(1, k - 1);
    }

    for(index l = k - 1; l > 0; l--)
      //_ARAGELI_LLL_RED(k, l - 1)
    {                                                        
        T v_1_2(1, 2);                                         
        if(v_1_2 < std::abs(Mu(k, l-1)))                         
        {                                                      
            T q = std::floor(Mu(k, l-1) + v_1_2);                  
            B.assign_col(k, B.copy_col(k) - B.copy_col(l-1) * q);  
            H.assign_col(k, H.copy_col(k) - H.copy_col(l-1) * q);  
            Mu(k, l-1) = Mu(k, l-1) - q;                             
            for (index i = 0; i < l-1; i++)                        
            Mu(k, i) = Mu(k, i) - q * Mu(l-1, i);                
        }                                                      
    }


  }

  return true;
}

#undef _ARAGELI_LLL_RED
#undef _ARAGELI_LLL_SWAP

}


#endif // #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE

---