Merge pull request #3992 from pleroy/LLL

An implementation of the LLL algorithm

numerics/lattices.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include "numerics/concepts.hpp"
+
+namespace principia {
+namespace numerics {
+namespace _lattices {
+namespace internal {
+
+using namespace principia::numerics::_concepts;
+
+template<typename Matrix>
+  requires two_dimensional<Matrix>
+Matrix LenstraLenstraLovász(Matrix const& L);
+
+}  // namespace internal
+
+using internal::LenstraLenstraLovász;
+
+}  // namespace _lattices
+}  // namespace numerics
+}  // namespace principia
+
+#include "numerics/lattices_body.hpp"

numerics/lattices_body.hpp

@@ -0,0 +1,93 @@
+#pragma once
+
+#include "numerics/lattices.hpp"
+
+#include <algorithm>
+
+#include "numerics/fixed_arrays.hpp"
+#include "numerics/matrix_computations.hpp"
+#include "numerics/matrix_views.hpp"
+#include "numerics/unbounded_arrays.hpp"
+#include "quantities/elementary_functions.hpp"
+
+namespace principia {
+namespace numerics {
+namespace _lattices {
+namespace internal {
+
+using namespace principia::numerics::_fixed_arrays;
+using namespace principia::numerics::_matrix_computations;
+using namespace principia::numerics::_matrix_views;
+using namespace principia::numerics::_unbounded_arrays;
+using namespace principia::quantities::_elementary_functions;
+
+template<typename Matrix>
+struct LenstraLenstraLovászGenerator;
+
+template<typename Scalar>
+struct LenstraLenstraLovászGenerator<UnboundedMatrix<Scalar>> {
+  using Vector = UnboundedVector<Scalar>;
+};
+
+template<typename Scalar, int rows, int columns>
+struct LenstraLenstraLovászGenerator<
+    FixedMatrix<Scalar, rows, columns>> {
+  using Vector = FixedVector<Scalar, rows>;
+};
+
+// This implements [HPS], theorem 7.71, figure 7.8.  Note that figures 7.9 and
+// 7.10 are supposedly more efficient, but they are significantly more
+// complicated.  If performance is an issue, we should look into recent
+// improvements of LLL.
+template<typename Matrix>
+  requires two_dimensional<Matrix>
+Matrix LenstraLenstraLovász(Matrix const& L) {
+  using G = LenstraLenstraLovászGenerator<Matrix>;
+  auto const n = L.columns();
+  auto const m = L.rows();
+  auto v = L;
+  for (int k = 1; k < n;) {
+    auto qr = UnitriangularGramSchmidt(v);
+    auto vₖ = ColumnView{.matrix = v,
+                        .first_row = 0,
+                        .last_row = m - 1,
+                        .column = k};
+    for (int j = k - 1; j >= 0; --j) {
+      auto const μₖⱼ = qr.R(j, k);
+      auto vⱼ = ColumnView{.matrix = v,
+                           .first_row = 0,
+                           .last_row = m - 1,
+                           .column = j};
+      auto const round_μₖⱼ = std::round(μₖⱼ);
+      if (round_μₖⱼ != 0) {
+        vₖ -= round_μₖⱼ * typename G::Vector(vⱼ);
+        qr = UnitriangularGramSchmidt(v);
+      }
+    }
+    auto const μₖₖ₋₁ = qr.R(k - 1, k);
+    auto v𐌟ₖ = ColumnView{.matrix = qr.Q,
+                         .first_row = 0,
+                         .last_row = m - 1,
+                         .column = k};
+    auto v𐌟ₖ₋₁ = ColumnView{.matrix = qr.Q,
+                           .first_row = 0,
+                           .last_row = m - 1,
+                           .column = k - 1};
+    if (v𐌟ₖ.Norm²() >= (0.75 - Pow<2>(μₖₖ₋₁)) * v𐌟ₖ₋₁.Norm²()) {
+      ++k;
+    } else {
+      auto vₖ₋₁ = ColumnView{.matrix = v,
+                            .first_row = 0,
+                            .last_row = m - 1,
+                            .column = k - 1};
+      SwapColumns(vₖ₋₁, vₖ);
+      k = std::max(k - 1, 1);
+    }
+  }
+  return v;
+}
+
+}  // namespace internal
+}  // namespace _lattices
+}  // namespace numerics
+}  // namespace principia

numerics/lattices_test.cpp

@@ -0,0 +1,40 @@
+#include "numerics/lattices.hpp"
+
+#include "gtest/gtest.h"
+#include "numerics/fixed_arrays.hpp"
+#include "testing_utilities/almost_equals.hpp"
+
+namespace principia {
+namespace numerics {
+namespace _lattices {
+
+using namespace principia::numerics::_fixed_arrays;
+using namespace principia::numerics::_lattices;
+using namespace principia::testing_utilities::_almost_equals;
+
+class LatticesTest : public ::testing::Test {};
+
+// [HPS14], example 7.75.
+TEST_F(LatticesTest, Example_7_75) {
+  FixedMatrix<double, 6, 6> l({19, 15, 43, 20,  0, 48,
+                                2, 42, 15, 44, 48, 33,
+                               32, 11,  0, 44, 35, 32,
+                               46,  0, 24,  0, 16,  9,
+                                3,  3,  4, 18, 31,  1,
+                               33, 24, 16, 15, 31, 29});
+
+  auto const reduced = LenstraLenstraLovász(l);
+  EXPECT_THAT(reduced,
+              AlmostEquals(
+                  (FixedMatrix<double, 6, 6>({  7, -20,  5,  -6, -10,   7,
+                                              -12,   4,  2,  -7, -24,   4,
+                                               -8,  -9, 33, -20,  21,  -9,
+                                                4,  16,  0, -21, -15, -11,
+                                               19,  13, 15,   8,  -6,   1,
+                                                9,  16, -9, -12, -11,  31})),
+                  0));
+}
+
+}  // namespace _lattices
+}  // namespace numerics
+}  // namespace principia

numerics/matrix_views.hpp

@@ -34,9 +34,6 @@ struct BlockView {
   template<typename T>
     requires two_dimensional<T> && same_elements_as<T, Matrix>
   BlockView& operator=(T const& right);
-  template<typename T>
-    requires two_dimensional<T> && same_elements_as<T, Matrix>
-  BlockView& operator=(T&& right);
 
   template<typename T>
     requires two_dimensional<T> && same_elements_as<T, Matrix>
@@ -69,9 +66,7 @@ struct ColumnView {
   template<typename T>
     requires one_dimensional<T> && same_elements_as<T, Matrix>
   ColumnView& operator=(T const& right);
-  template<typename T>
-    requires one_dimensional<T> && same_elements_as<T, Matrix>
-  ColumnView& operator=(T&& right);
+  ColumnView& operator=(ColumnView const& right);
 
   template<typename T>
     requires one_dimensional<T> && same_elements_as<T, Matrix>
@@ -88,6 +83,12 @@ struct ColumnView {
   constexpr int size() const;
 };
 
+// TODO(phl): This should probably be just |swap|.  The semantics of BlockView
+// and ColumnView do imply an implicit dereferencing, so swapping should work
+// the same.
+template<typename Matrix>
+void SwapColumns(ColumnView<Matrix>& m1, ColumnView<Matrix>& m2);
+
 template<typename LMatrix, typename RMatrix>
 Product<typename LMatrix::Scalar, typename RMatrix::Scalar> operator*(
     TransposedView<ColumnView<LMatrix>> const& left,
@@ -117,6 +118,7 @@ std::ostream& operator<<(std::ostream& out,
 
 using internal::BlockView;
 using internal::ColumnView;
+using internal::SwapColumns;
 
 }  // namespace _matrix_views
 }  // namespace numerics

numerics/matrix_views_body.hpp

@@ -49,21 +49,6 @@ BlockView<Matrix>& BlockView<Matrix>::operator=(T const& right) {
   return *this;
 }
 
-template<typename Matrix>
-  requires two_dimensional<Matrix>
-template<typename T>
-  requires two_dimensional<T> && same_elements_as<T, Matrix>
-BlockView<Matrix>& BlockView<Matrix>::operator=(T&& right) {
-  DCHECK_EQ(rows(), right.rows());
-  DCHECK_EQ(columns(), right.columns());
-  for (int i = 0; i < rows(); ++i) {
-    for (int j = 0; j < columns(); ++j) {
-      matrix(first_row + i, first_column + j) = std::move(right(i, j));
-    }
-  }
-  return *this;
-}
-
 template<typename Matrix>
   requires two_dimensional<Matrix>
 template<typename T>
@@ -158,12 +143,10 @@ ColumnView<Matrix>& ColumnView<Matrix>::operator=(T const& right) {
 
 template<typename Matrix>
   requires two_dimensional<Matrix>
-template<typename T>
-  requires one_dimensional<T> && same_elements_as<T, Matrix>
-ColumnView<Matrix>& ColumnView<Matrix>::operator=(T&& right) {
+ColumnView<Matrix>& ColumnView<Matrix>::operator=(ColumnView const& right) {
   DCHECK_EQ(size(), right.size());
   for (int i = 0; i < size(); ++i) {
-    matrix(first_row + i, column) = std::move(right[i]);
+    matrix(first_row + i, column) = right[i];
   }
   return *this;
 }
@@ -232,6 +215,15 @@ constexpr auto ColumnView<Matrix>::size() const -> int {
   return last_row - first_row + 1;
 }
 
+template<typename Matrix>
+void SwapColumns(ColumnView<Matrix>& m1, ColumnView<Matrix>& m2) {
+  DCHECK_EQ(m1.first_row, m2.first_row);
+  DCHECK_EQ(m1.last_row, m2.last_row);
+  for (int i = 0; i < m1.size(); ++i) {
+    std::swap(m1[i], m2[i]);
+  }
+}
+
 template<typename LMatrix, typename RMatrix>
 Product<typename LMatrix::Scalar, typename RMatrix::Scalar> operator*(
     TransposedView<ColumnView<LMatrix>> const& left,

numerics/matrix_views_test.cpp

@@ -142,6 +142,22 @@ TEST_F(MatrixViewsTest, ColumnView_Indexing) {
   EXPECT_EQ(-9, cum34[1]);
 }
 
+TEST_F(MatrixViewsTest, ColumnView_Assignment) {
+  ColumnView<FixedMatrix<double, 3, 4>> cfm34{.matrix = fm34_,
+                                              .first_row = 1,
+                                              .last_row = 2,
+                                              .column = 2};
+  ColumnView<UnboundedMatrix<double>> cum34{.matrix = um34_,
+                                            .first_row = 1,
+                                            .last_row = 2,
+                                            .column = 3};
+  cum34 = cfm34;
+  EXPECT_EQ(9, cum34[0]);
+  EXPECT_EQ(-2, cum34[1]);
+
+  cfm34 = cfm34;
+}
+
 TEST_F(MatrixViewsTest, ColumnView_Addition) {
   ColumnView<FixedMatrix<double, 3, 4>> cfm34{.matrix = fm34_,
                                               .first_row = 1,

numerics/numerics.vcxproj

@@ -24,6 +24,8 @@
     <ClInclude Include="gradient_descent_body.hpp" />
     <ClInclude Include="hermite2.hpp" />
     <ClInclude Include="hermite2_body.hpp" />
+    <ClInclude Include="lattices.hpp" />
+    <ClInclude Include="lattices_body.hpp" />
     <ClInclude Include="matrix_computations.hpp" />
     <ClInclude Include="matrix_computations_body.hpp" />
     <ClInclude Include="combinatorics.hpp" />
@@ -119,6 +121,7 @@
     <ClCompile Include="gradient_descent_test.cpp" />
     <ClCompile Include="hermite2_test.cpp" />
     <ClCompile Include="hermite3_test.cpp" />
+    <ClCompile Include="lattices_test.cpp" />
     <ClCompile Include="legendre_test.cpp" />
     <ClCompile Include="matrix_computations_test.cpp" />
     <ClCompile Include="matrix_views_test.cpp" />

numerics/numerics.vcxproj.filters

@@ -113,9 +113,6 @@
     <ClInclude Include="apodization.hpp">
       <Filter>Header Files</Filter>
     </ClInclude>
-    <ClInclude Include="apodization_body.hpp">
-      <Filter>Source Files</Filter>
-    </ClInclude>
     <ClInclude Include="fast_fourier_transform.hpp">
       <Filter>Header Files</Filter>
     </ClInclude>
@@ -266,6 +263,12 @@
     <ClInclude Include="matrix_views_body.hpp">
       <Filter>Source Files</Filter>
     </ClInclude>
+    <ClInclude Include="lattices.hpp">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="lattices_body.hpp">
+      <Filter>Source Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="fixed_arrays_test.cpp">
@@ -388,6 +391,12 @@
     <ClCompile Include="matrix_views_test.cpp">
       <Filter>Test Files</Filter>
     </ClCompile>
+    <ClCompile Include="apodization.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="lattices_test.cpp">
+      <Filter>Test Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <Text Include="xgscd.proto.txt">