gso_interface.h

Go to the documentation of this file.

/* Copyright (C) 2005-2008 Damien Stehle.
   Copyright (C) 2007 David Cade.
   Copyright (C) 2011 Xavier Pujol.
   Copyright (C) 2019 Koen de Boer
 
   This file is part of fplll. fplll is free software: you
   can redistribute it and/or modify it under the terms of the GNU Lesser
   General Public License as published by the Free Software Foundation,
   either version 2.1 of the License, or (at your option) any later version.
 
   fplll is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
   GNU Lesser General Public License for more details.
 
   You should have received a copy of the GNU Lesser General Public License
   along with fplll. If not, see <http://www.gnu.org/licenses/>. */
 
#ifndef FPLLL_GSOInterface_H
#define FPLLL_GSOInterface_H
 
#include "nr/matrix.h"
 
FPLLL_BEGIN_NAMESPACE
 
enum MatGSOInterfaceFlags
{
  GSO_DEFAULT       = 0,
  GSO_INT_GRAM      = 1,
  GSO_ROW_EXPO      = 2,
  GSO_OP_FORCE_LONG = 4
};
 
template <class FT>
void adjust_radius_to_gh_bound(FT &max_dist, long max_dist_expo, int block_size, const FT &root_det,
                               double gh_factor);
 
template <class ZT, class FT> class MatGSOInterface
{
public:
  virtual ~MatGSOInterface();
 
  MatGSOInterface(Matrix<ZT> &arg_u, Matrix<ZT> &arg_uinv_t, int flags)
      : enable_int_gram(flags & GSO_INT_GRAM), enable_row_expo(flags & GSO_ROW_EXPO),
        enable_transform(arg_u.get_rows() > 0), enable_inverse_transform(arg_uinv_t.get_rows() > 0),
        row_op_force_long(flags & GSO_OP_FORCE_LONG), u(arg_u), u_inv_t(arg_uinv_t),
        n_known_rows(0), n_source_rows(0), n_known_cols(0), cols_locked(false), alloc_dim(0),
        gptr(nullptr)
  {
#ifdef DEBUG
    row_op_first = row_op_last = -1;
#endif
  }
 
  int d;
 
  // Matrix<ZT> &b;
 
  virtual ZT &sqnorm_coordinates(ZT &sqnorm, vector<ZT> coordinates) = 0;
 
  virtual long get_max_exp_of_b() = 0;
 
  virtual bool b_row_is_zero(int i) = 0;
 
  virtual int get_cols_of_b() const = 0;
 
  virtual int get_rows_of_b() const = 0;
 
  virtual void negate_row_of_b(int i) = 0;
 
  vector<long> row_expo;
 
  inline void row_op_begin(int first, int last);
 
  void row_op_end(int first, int last);
 
  virtual FT &get_gram(FT &f, int i, int j) = 0;
 
  virtual ZT &get_int_gram(ZT &z, int i, int j) = 0;
 
  const Matrix<FT> &get_mu_matrix() const { return mu; }
 
  const Matrix<FT> &get_r_matrix() const { return r; }
 
  const Matrix<ZT> &get_g_matrix() const
  {
    if (gptr == nullptr)
    {
      throw std::runtime_error("Error: gptr == nullpointer.");
    }
    return *gptr;
  }
 
  inline const FT &get_mu_exp(int i, int j, long &expo) const;
  inline const FT &get_mu_exp(int i, int j) const;
 
  inline FT &get_mu(FT &f, int i, int j);
 
  ZT get_max_gram();
 
  FT get_max_bstar();
 
  inline const FT &get_r_exp(int i, int j, long &expo) const;
  inline const FT &get_r_exp(int i, int j) const;
 
  inline FT &get_r(FT &f, int i, int j);
 
  long get_max_mu_exp(int i, int n_columns);
 
  bool update_gso_row(int i, int last_j);
 
  inline bool update_gso_row(int i);
 
  inline bool update_gso();
 
  inline void discover_all_rows();
 
  void set_r(int i, int j, FT &f);
 
  virtual void move_row(int old_r, int new_r) = 0;
 
  virtual inline void row_addmul(int i, int j, const FT &x);
 
  virtual void row_addmul_we(int i, int j, const FT &x, long expo_add) = 0;
 
  // b[i] += b[j] / b[i] -= b[j] (i > j)
  virtual void row_add(int i, int j) = 0;
  virtual void row_sub(int i, int j) = 0;
 
  void lock_cols();
  void unlock_cols();
 
  inline void create_row();
  virtual void create_rows(int n_new_rows) = 0;
 
  inline void remove_last_row();
  virtual void remove_last_rows(int n_removed_rows) = 0;
 
  void apply_transform(const Matrix<FT> &transform, int src_base, int target_base);
 
  void apply_transform(const Matrix<FT> &transform, int src_base)
  {
    apply_transform(transform, src_base, src_base);
  }
 
  inline void dump_mu_d(double *mu, int offset = 0, int block_size = -1);
  inline void dump_mu_d(vector<double> mu, int offset = 0, int block_size = -1);
 
  inline void dump_r_d(double *r, int offset = 0, int block_size = -1);
  inline void dump_r_d(vector<double> &r, int offset = 0, int block_size = -1);
 
  double get_current_slope(int start_row, int stop_row);
 
  FT get_root_det(int start_row, int end_row);
 
  FT get_log_det(int start_row, int end_row);
 
  FT get_slide_potential(int start_row, int end_row, int block_size);
 
  inline void print_mu_r_g(ostream &os);
 
  void virtual babai(vector<ZT> &v, int start = 0, int dimension = -1);
 
  void virtual babai(vector<ZT> &w, const vector<FT> &v, int start = 0, int dimension = 1);
 
  const bool enable_int_gram;
 
  const bool enable_row_expo;
 
  const bool enable_transform;
 
  const bool enable_inverse_transform;
 
  const bool row_op_force_long;
 
protected:
  virtual void size_increased() = 0;
 
  /* Increases known rows and invalidates the last
   * gram row (gf) when enable_int_gram is false.
   * When enable_int_gram is true, it computes
   * the new inner products for the last gram
   * row (g).
   */
  virtual void discover_row() = 0;
 
  // Marks mu(i, j) and r(i, j) as invalid for j >= new_valid_cols
  void invalidate_gso_row(int i, int new_valid_cols = 0);
  /*  Upates the i-th row of bf. It does not invalidate anything, so the caller
   * must take into account that it might change row_expo.
   */
  virtual void update_bf(int i) = 0;
  /* Marks g(i, j) for all j <= i (but NOT for j > i) */
  virtual void invalidate_gram_row(int i) = 0;
 
  // b[i] <- b[i] + x * b[j] (i > j)
  virtual void row_addmul_si(int i, int j, long x) = 0;
  // b[i] <- b[i] + (2^expo * x) * b[j] (i > j)
 
  // TODO check if these must be scratched here!
  // virtual void row_addmul_si_2exp(int i, int j, long x, long expo) = 0;
  // virtual void row_addmul_2exp(int i, int j, const ZT &x, long expo) = 0;
 
public:
  // made public for lll.cpp and bkz.cpp
  void symmetrize_g();
 
  // Made public for bkz.cpp (dbkz_postprocessing)
  // b[i] <-> b[j] (i < j)
  virtual void row_swap(int i, int j) = 0;
 
protected:
  inline ZT &sym_g(int i, int j);
 
  /* Floating-point representation of the basis. It is used when
     enable_int_gram=false. */
  Matrix<FT> bf;
 
  Matrix<ZT> &u;        // Transform
  Matrix<ZT> &u_inv_t;  // Transposed inverse transform
 
  // init_row_size[i] = (last non-zero column in the i-th row of b) + 1
  vector<int> init_row_size;
 
  // bf[i], g[i], gf[i], mu[i] and r[i] are invalid for i >= n_known_rows
  int n_known_rows;
  int n_source_rows;  // Known rows before the beginning of early reduction
  int n_known_cols;
  bool cols_locked;
  int alloc_dim;
 
  Matrix<FT> mu;
 
  Matrix<FT> r;
 
public:
  /* Gram matrix (dot products of basis vectors, lower triangular matrix)
   g(i, j) is valid if 0 <= i < n_known_rows and j <= i */
  Matrix<ZT> *gptr;
  //  Matrix<ZT> g;
 
protected:
  Matrix<FT> gf;
 
  /* Number of valid columns of the i-th row of mu and r.
     Valid only for 0 <= i < n_known_rows */
  vector<int> gso_valid_cols;
 
  /* Used by update_gso_row (+ update_gso), get_max_mu_exp and row_addmul_we. */
  FT ftmp1, ftmp2;
  /* Used by row_add, row_sub, row_addmul_si_2exp, row_addmul_2exp and
     indirectly by row_addmul. */
  ZT ztmp1;
  /* Used by row_addmul. */
  ZT ztmp2;
  /* Used by update_bf. */
  vector<long> tmp_col_expo;
 
#ifdef DEBUG
  /* Used only in debug mode. */
  int row_op_first, row_op_last;
  bool in_row_op_range(int i) { return i >= row_op_first && i < row_op_last; }
#endif
};
 
template <class ZT, class FT> inline MatGSOInterface<ZT, FT>::~MatGSOInterface()
{
  // delete gptr;
}
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::symmetrize_g()
{
  if (gptr == nullptr)
  {
    throw std::runtime_error("Error: gptr is equal to the nullpointer.");
  }
  Matrix<ZT> &gr = *gptr;
  for (int i = 0; i < d; i++)
  {
    for (int j = 0; j < d; j++)
    {
      gr(i, j) = sym_g(i, j);
    }
  }
}
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::print_mu_r_g(ostream &os)
{
  os << "mu = " << endl;
  mu.print(os);
  os << endl << "r = " << endl;
  r.print(os);
  os << endl;
  if (gptr != nullptr)
  {
    os << "g = " << endl;
    symmetrize_g();
    gptr->print(os);
    os << endl << endl;
  }
}
 
template <class ZT, class FT> inline ZT &MatGSOInterface<ZT, FT>::sym_g(int i, int j)
{
  if (gptr == nullptr)
  {
    throw std::runtime_error("Error: gptr is equal to the nullpointer.");
  }
  Matrix<ZT> &gr = *gptr;
  return (i >= j) ? gr(i, j) : gr(j, i);
}
 
template <class ZT, class FT>
inline const FT &MatGSOInterface<ZT, FT>::get_mu_exp(int i, int j, long &expo) const
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < i && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  if (enable_row_expo)
    expo = row_expo[i] - row_expo[j];
  else
    expo = 0;
  return mu(i, j);
}
 
template <class ZT, class FT>
inline const FT &MatGSOInterface<ZT, FT>::get_mu_exp(int i, int j) const
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < i && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  return mu(i, j);
}
 
template <class ZT, class FT> inline FT &MatGSOInterface<ZT, FT>::get_mu(FT &f, int i, int j)
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < i && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  f = mu(i, j);
  if (enable_row_expo)
    f.mul_2si(f, row_expo[i] - row_expo[j]);
  return f;
}
 
template <class ZT, class FT>
inline const FT &MatGSOInterface<ZT, FT>::get_r_exp(int i, int j, long &expo) const
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  if (enable_row_expo)
    expo = row_expo[i] + row_expo[j];
  else
    expo = 0;
  return r(i, j);
}
 
template <class ZT, class FT>
inline const FT &MatGSOInterface<ZT, FT>::get_r_exp(int i, int j) const
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  return r(i, j);
}
 
template <class ZT, class FT> inline FT &MatGSOInterface<ZT, FT>::get_r(FT &f, int i, int j)
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && j >= 0 && j < gso_valid_cols[i] &&
                    !in_row_op_range(i));
  f = r(i, j);
  if (enable_row_expo)
    f.mul_2si(f, row_expo[i] + row_expo[j]);
  return f;
}
 
template <class ZT, class FT> inline bool MatGSOInterface<ZT, FT>::update_gso_row(int i)
{
  return update_gso_row(i, i);
}
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::set_r(int i, int j, FT &f)
{
  FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && gso_valid_cols[i] >= j && j >= 0 &&
                    j < n_source_rows);
  r(i, j) = f;
  if (gso_valid_cols[i] == j)
    gso_valid_cols[i]++;
}
 
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::row_addmul(int i, int j, const FT &x)
{
  row_addmul_we(i, j, x, 0);
}
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::create_row() { create_rows(1); }
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::remove_last_row()
{
  remove_last_rows(1);
}
 
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::discover_all_rows()
{
  while (n_known_rows < d)
    discover_row();
}
 
template <class ZT, class FT> inline bool MatGSOInterface<ZT, FT>::update_gso()
{
  for (int i = 0; i < d; i++)
  {
    if (!update_gso_row(i))
      return false;
  }
  return true;
}
 
#ifdef DEBUG
template <class ZT, class FT> inline void MatGSOInterface<ZT, FT>::row_op_begin(int first, int last)
{
  FPLLL_DEBUG_CHECK(row_op_first == -1);
  row_op_first = first;
  row_op_last  = last;
}
#else
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::row_op_begin(int /*first*/, int /*last*/)
{
}
#endif
 
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::dump_mu_d(double *mu, int offset, int block_size)
{
  FT e;
  if (block_size <= 0)
  {
    block_size = get_rows_of_b();
  }
 
  for (int i = 0; i < block_size; ++i)
  {
    for (int j = 0; j < block_size; ++j)
    {
      get_mu(e, offset + i, offset + j);
      mu[i * block_size + j] = e.get_d();
    }
  }
}
 
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::dump_mu_d(vector<double> mu, int offset, int block_size)
{
  FT e;
  if (block_size <= 0)
  {
    block_size = get_rows_of_b();
  }
 
  mu.reserve(mu.size() + block_size * block_size);
  for (int i = 0; i < block_size; ++i)
  {
    for (int j = 0; j < block_size; ++j)
    {
      get_mu(e, offset + i, offset + j);
      mu.push_back(e.get_d());
    }
  }
}
 
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::dump_r_d(double *r, int offset, int block_size)
{
  FT e;
  if (block_size <= 0)
  {
    block_size = get_rows_of_b();
  }
 
  for (int i = 0; i < block_size; ++i)
  {
    get_r(e, offset + i, offset + i);
    r[i] = e.get_d();
  }
}
 
template <class ZT, class FT>
inline void MatGSOInterface<ZT, FT>::dump_r_d(vector<double> &r, int offset, int block_size)
{
  FT e;
  if (block_size <= 0)
  {
    block_size = get_rows_of_b();
  }
 
  r.reserve(r.size() + block_size * block_size);
  for (int i = 0; i < block_size; ++i)
  {
    get_r(e, offset + i, offset + i);
    r.push_back(e.get_d());
  }
}
 
FPLLL_END_NAMESPACE
 
#endif