diff --git a/src_basic/Basic_functions.h b/src_basic/Basic_functions.h
index 455c574..4f09b53 100644
--- a/src_basic/Basic_functions.h
+++ b/src_basic/Basic_functions.h
@@ -68,6 +68,13 @@ template <typename T> T T_abs(T const &eVal) {
   return fVal;
 }
 
+template <typename T> T T_sign(T const &eVal) {
+  T fVal(1);
+  if (eVal < 0)
+    fVal = T(-1);
+  return fVal;
+}
+
 // use std::min for generic types (int, long, float, ...)
 template <typename T> T T_min(T const &eVal1, T const &eVal2) {
   if (eVal1 > eVal2)
diff --git a/src_matrix/MAT_Matrix.h b/src_matrix/MAT_Matrix.h
index a2ab7a9..d9c478a 100644
--- a/src_matrix/MAT_Matrix.h
+++ b/src_matrix/MAT_Matrix.h
@@ -2766,6 +2766,22 @@ template <typename T> MyVector<T> SignCanonicalizeVector(const MyVector<T> &V) {
   throw TerminalException{1};
 }
 
+template<typename T>
+MyVector<T> GetDiagonal(MyMatrix<T> const& M) {
+  int n_col = M.cols();
+  int n_row = M.rows();
+  if (n_row != n_col) {
+    std::cerr << "The matrix is rectangular but not square, so diagonal is not defined\n";
+    throw TerminalException{1};
+  }
+  MyMatrix<T> V(n_row);
+  for (int i=0; i<n_row; i++) {
+    V(i) = M(i, i);
+  }
+  return V;
+}
+
+
 // clang-format off
 #endif  // SRC_MATRIX_MAT_MATRIX_H_
 // clang-format on
diff --git a/src_number/TestQuadraticResidue.cpp b/src_number/TestQuadraticResidue.cpp
index 78960a0..141d672 100644
--- a/src_number/TestQuadraticResidue.cpp
+++ b/src_number/TestQuadraticResidue.cpp
@@ -17,7 +17,6 @@ void process(int TheMod) {
   for (int u=0; u<TheMod; u++) {
     T u_T(u);
     bool test = is_quadratic_residue(u_T, TheMod_T);
-    std::cerr << "u=" << u << " test=" << test << "\n";
     if (test) {
       n_quad += 1;
     }
@@ -36,6 +35,7 @@ void process(int TheMod) {
 int main(int argc, char *argv[]) {
   try {
     int TheMod = 23;
+    std::cerr << "TheMod=" << TheMod << "\n";
     process<mpz_class>(TheMod);
     process<SafeInt64>(TheMod);
     //    process<boost::multiprecision::cpp_int>(TheMod);
diff --git a/src_number/factorizations.h b/src_number/factorizations.h
index 7f3cba1..5fef6fb 100644
--- a/src_number/factorizations.h
+++ b/src_number/factorizations.h
@@ -95,6 +95,18 @@ template <typename T> std::vector<T> FactorsInt(T const &N) {
   }
 }
 
+template <typename T> std::map<T, size_t> FactorsIntMap(T const &N) {
+  static_assert(is_implementation_of_Z<T>::value, "Requires T to be a Z ring");
+  std::vector<T> vect = FactorsInt(N);
+  std::map<T,size_t> map;
+  for (auto & eV : vect) {
+    map[eV] += 1;
+  }
+  return map;
+}
+
+
+
 template <typename T>
 std::vector<T> GetAllFactors(std::map<T, int> const &eMap) {
   static_assert(is_implementation_of_Z<T>::value, "Requires T to be a Z ring");
diff --git a/src_number/quadratic_residue.h b/src_number/quadratic_residue.h
index fb4b0d7..2b53060 100644
--- a/src_number/quadratic_residue.h
+++ b/src_number/quadratic_residue.h
@@ -9,6 +9,7 @@
 #include <vector>
 // clang-format on
 
+// This is an exhaustive search that works even if m is not prime.
 template <typename T>
 std::optional<T> find_quadratic_residue(T const& a, T const& m_in) {
   static_assert(is_implementation_of_Z<T>::value, "Requires T to be a Z ring");
@@ -26,11 +27,8 @@ std::optional<T> find_quadratic_residue(T const& a, T const& m_in) {
   T x(0);
   T TwoXpOne(1);
   T xSqr(0);
-  //  std::cerr << "m=" << m << " a=" << a << " a_mod=" << a_mod << "\n";
   while (x != upper) {
-    //    std::cerr << "  x=" << x << " xSqr=" << xSqr << "\n";
     if (xSqr == a_mod) {
-      //      std::cerr << "Returning x=" << x << "\n";
       return x;
     }
     xSqr += TwoXpOne;