cleaned tpp files

src/F4_launcher.h

 #ifndef F4_LAUNCHER_H_
 #define F4_LAUNCHER_H_
 
+#include <string>
+#include <vector>
 #include "./algos.h"
 #include "./polynomial.h"
 #include "./parser.h"
@@ -16,8 +18,57 @@ void field_init();
 template <class T>
 void field_clear();
 
-}  // namespace tinygb
+template <class T>
+void launchF4(istream &in, ostream &out) {
+  log(LOG_INFO) << "begin parsing" << endl;
+  Parser<T> pp(in);
+  vector<polynomial<T> > sys = pp.parse();
+  log(LOG_INFO) << "parsing done" << endl;
+  vector<polynomial<T> > G;
+  F4<T>(sys, G);
+  for (typename std::vector<polynomial <T>>::const_iterator it = G.begin(); \
+      it != G.end(); ++it) {
+    it->print(out);
+    out << std::endl;
+  }
+}
+
+void initF4(istream &in, ostream &out) {
+  log(LOG_INFO) << "Initialization" << endl;
+  string s;
+  in >> s;
+#ifndef NBVAR
+  NBVAR = (unsigned)atoi(s.c_str());
+#endif
+  in >> s;
+  // set cardinality of finite field
+  // at most UINT_MAX at the moment
+  // TODO(pj): multiprecision
+  mpz_set_str(cardfield, s.c_str(), 10);
+  log(LOG_INFO) << NBVAR << " variables" << endl;
 
-#include "F4_launcher.tpp"
+  mpz_t max64b;
+  mpz_init(max64b);
+  mpz_set_str(max64b, "18446744073709551615", 10);
+  out << s << std::endl;
+  out << NBVAR << std::endl;
+  if (mpz_cmp_ui(cardfield, 65521) <= 0) {
+    mpfq_p_0_5_field_init(GFp_05::k);
+    mpfq_p_0_5_field_specify(GFp_05::k, MPFQ_PRIME_MPZ, cardfield);
+    BLASFLAG = true;
+    launchF4<GFp_05>(in, out);
+    mpfq_p_0_5_field_clear(GFp_1::k);  // clear field
+  } else if (mpz_cmp(cardfield, max64b) <= 0) {
+    mpfq_pm_1_field_init(GFp_1::k);
+    mpfq_pm_1_field_specify(GFp_1::k, MPFQ_PRIME_MPZ, cardfield);
+    launchF4<GFp_1>(in, out);
+    mpfq_pm_1_field_clear(GFp_1::k);
+  } else {
+    log(LOG_ERROR) << "field not supported" << endl;
+  }
+  mpz_clear(max64b);
+}
+
+}  // namespace tinygb
 
 #endif  // F4_LAUNCHER_H_

src/F4_launcher.tpp

-#include <string>
-#include <vector>
-
-namespace tinygb {
-
-void initF4(std::istream& in, std::ostream& out);
-
-template <class T>
-void field_init();
-
-template <class T>
-void field_clear();
-
-template <class T>
-void launchF4(istream &in, ostream &out) {
-  log(LOG_INFO) << "begin parsing" << endl;
-  Parser<T> pp(in);
-  vector<polynomial<T> > sys = pp.parse();
-  log(LOG_INFO) << "parsing done" << endl;
-  vector<polynomial<T> > G;
-  F4<T>(sys, G);
-  for (typename std::vector<polynomial <T>>::const_iterator it = G.begin(); \
-      it != G.end(); ++it) {
-    it->print(out);
-    out << std::endl;
-  }
-}
-
-void initF4(istream &in, ostream &out) {
-  log(LOG_INFO) << "Initialization" << endl;
-  string s;
-  in >> s;
-#ifndef NBVAR
-  NBVAR = (unsigned)atoi(s.c_str());
-#endif
-  in >> s;
-  // set cardinality of finite field
-  // at most UINT_MAX at the moment
-  // TODO(pj): multiprecision
-  mpz_set_str(cardfield, s.c_str(), 10);
-  log(LOG_INFO) << NBVAR << " variables" << endl;
-
-  mpz_t max64b;
-  mpz_init(max64b);
-  mpz_set_str(max64b, "18446744073709551615", 10);
-  out << s << std::endl;
-  out << NBVAR << std::endl;
-  if (mpz_cmp_ui(cardfield, 65521) <= 0) {
-    mpfq_p_0_5_field_init(GFp_05::k);
-    mpfq_p_0_5_field_specify(GFp_05::k, MPFQ_PRIME_MPZ, cardfield);
-    BLASFLAG = true;
-    launchF4<GFp_05>(in, out);
-    mpfq_p_0_5_field_clear(GFp_1::k);  // clear field
-  } else if (mpz_cmp(cardfield, max64b) <= 0) {
-    mpfq_pm_1_field_init(GFp_1::k);
-    mpfq_pm_1_field_specify(GFp_1::k, MPFQ_PRIME_MPZ, cardfield);
-    launchF4<GFp_1>(in, out);
-    mpfq_pm_1_field_clear(GFp_1::k);
-  } else {
-    log(LOG_ERROR) << "field not supported" << endl;
-  }
-  mpz_clear(max64b);
-}
-
-}  // namespace tinygb

src/algos.h

 #ifndef ALGOS_H_
 #define ALGOS_H_
 
+#include <algorithm>
+#include <list>
+#include <map>
+#include <set>
+#include <utility>
 #include <vector>
+#include "./critpair.h"
+#include "./linalg.h"
 #include "./polynomial.h"
+#include "./timer.h"
 
 namespace tinygb {
 
@@ -13,7 +21,321 @@ void F4(const std::vector<polynomial<T> > &sys,
 template <class T>
 bool gb_is_minimal(const std::vector<polynomial<T> >&);
 
+// test whether a GB is minimal
+template <class T>
+bool gb_is_minimal(const vector<polynomial<T> > &GB) {
+  for (unsigned i = 0; i < GB.size(); ++i)
+    for (unsigned j = 0; j < GB.size(); ++j)
+      if (i != j && GB[i].LM().isDivisibleBy(GB[j].LM()))
+        return false;
+  return true;
+}
+
+// computes a minimal GB from a GB
+template <class T>
+bool GB_minimization(vector<polymat<T> > &G) {
+  bool flag = false;
+  for (unsigned i = 0; i < G.size() - 1; ++i) {
+    if (G[i].get_LM().isDivisibleBy(G[G.size()-1].get_LM())) {
+      G.erase(G.begin() + i);
+      --i;
+      flag = true;
+    }
+  }
+  return flag;
+}
+
+// Gebauer/Moller installation of Buchberger's criteria (see description
+// Vanessa Vitse's thesis)
+template <class T>
+void gebmol_gb_insert(list<critpair<T> > &vCrit, vector<polymat<T> > &G,
+    const polymat<T> &f) {
+  for (auto it = vCrit.begin(); it != vCrit.end(); ++it) {
+    if (it->LCM() != it->get_S1().second.get_LM().LCM(f.get_LM()) &&
+        it->LCM() != it->get_S2().second.get_LM().LCM(f.get_LM()) &&
+        it->LCM().isDivisibleBy(it->get_S1().second.get_LM().LCM(f.get_LM()))
+        &&
+        it->LCM().isDivisibleBy(it->get_S2().second.get_LM().LCM(f.get_LM()))) {
+      it = vCrit.erase(it);
+      it--;
+    }
+  }
+  vector<critpair<T> > S0, S1, S2;
+  for (unsigned i = 0; i < G.size(); ++i) {
+    if (G[i].get_LM().LCM(f.get_LM()) == G[i].get_LM() * f.get_LM())
+      S0.push_back(critpair<T>(f, G[i]));
+    else
+      S1.push_back(critpair<T>(f, G[i]));
+  }
+  while (!S1.empty()) {
+    critpair<T> cp1 = S1[S1.size() - 1];
+    S1.erase(S1.end()-1);
+    bool flag = false;
+    monomial tmp_lcm = cp1.LCM();
+    for (auto& it : std::as_const(S0)) {
+      if (tmp_lcm.isDivisibleBy(it.LCM())) {
+        flag = true;
+        break;
+      }
+    }
+    if (!flag) {
+      for (auto& it : std::as_const(S1)) {
+        if (tmp_lcm.isDivisibleBy(it.LCM())) {
+          flag = true;
+          break;
+        }
+      }
+    }
+    if (!flag) {
+      for (auto& it : std::as_const(S2)) {
+        if (tmp_lcm.isDivisibleBy(it.LCM())) {
+          flag = true;
+          break;
+        }
+      }
+    }
+    if (!flag) {
+      S2.push_back(cp1);
+      vCrit.push_back(cp1);
+    }
+  }
+  if (f.istopred(G))
+    return;
+  G.push_back(f);
+  GB_minimization(G);
+}
+
+// selects critical pairs to reduce in a F4 step from the current list of
+// critical pairs
+template <class T>
+unsigned F4_critpair_selection(list<critpair<T> > &vCrit_step,
+    list<critpair<T> > &vCrit) {
+  unsigned stepdegree;
+  vCrit_step.clear();
+
+#ifdef BILIN
+  // TODO(pj): make const
+  // TODO(pj): improve
+  stepdegree = vCrit.front().degree_bilin();
+  for (typename list<critpair<T> >::iterator it = vCrit.begin();
+      it != vCrit.end(); ++it)
+    if (it->degree_bilin() < stepdegree)
+      stepdegree = it->degree_bilin();
+#else
+  stepdegree = vCrit.front().degree();
+  for (auto& it : std::as_const(vCrit))
+    if (it.degree() < stepdegree)
+      stepdegree = it.degree();
+#endif
+
+  for (typename list<critpair<T> >::iterator it = vCrit.begin();
+      it != vCrit.end(); ++it)
+#ifdef BILIN
+    if (it->degree_bilin() == stepdegree) {
+      vCrit_step.push_back(*it);
+      it = vCrit.erase(it);
+      it--;
+    }
+#else
+  if (it->degree() == stepdegree) {
+    vCrit_step.push_back(*it);
+    it = vCrit.erase(it);
+    it--;
+  }
+#endif
+  return stepdegree;
+}
+
+// interreduce a list of polynomials (linear algebra)
+template <class T>
+void interreduce(list<polymat<T> > &Ginit, vector<matrixF4<T>* > &vMat,
+    list<polynomial<T> > &_P) {
+  log(LOG_INFO) << "interreducing --- ";
+  for (auto& it : _P)
+    it.set_monic();
+
+  vector<polymat<T> > res;
+  matrixF4<T> *M = new matrixF4<T>(_P);
+  vMat.push_back(M);
+  M->row_echelon();
+  polymat<T>::construct_polymats(res, M);
+  Ginit.clear();
+  for (unsigned i = 0; i < res.size(); ++i)
+    Ginit.push_back(res[i]);
+  log(LOG_INFO) << "DONE" << endl;
+}
+
+// selects reductors to reduce a list of critical pairs
+template <class T>
+void F4_reductors_selection(const list<critpair<T> > &vCrit_step,
+    const list<polymat<T> > &globred,
+    vector<pair<monomial, polymat<T> > > &reductors) {
+  reductors.clear();
+  polymat<T> reductor;
+  typename list<polymat<T> >::const_iterator itred;
+  polynomial<T> tmp;
+  set<monomial> lM;
+  for (auto& it : std::as_const(vCrit_step)) {
+    tmp = it.eval();
+    for (auto& it2 : std::as_const(tmp.get_vp()))
+      lM.insert(it2.first);
+  }
+
+  bool flag;
+  while (!lM.empty()) {
+    monomial m = *(lM.rbegin());
+    flag = false;
+    for (auto it  = std::as_const(globred).begin();
+              it != std::as_const(globred).end();
+              ++it) {
+      if (m.isDivisibleBy(it->get_LM())) {
+        if (!flag) {
+          flag = true;
+          itred = it;
+        } else if (it->nb_non_zero() < itred->nb_non_zero()) {
+          itred = it;
+        }
+      }
+    }
+    if (flag) {
+      reductor = polymat<T> (*itred);
+      reductors.push_back(pair<monomial, polymat<T> >(m/reductor.get_LM(),
+            reductor));
+      reductor.merge_sets(lM, m/reductor.get_LM());
+    }
+    lM.erase(--lM.end());
+  }
+}
+
+// toplevel F4
+template <class T>
+void F4(const vector<polynomial<T> > &sys, vector<polynomial<T> > &res) {
+  log(LOG_INFO) << "beginning F4: " << endl;
+
+  for (auto it = sys.begin(); it != sys.end(); ++it)
+    log(LOG_INFO) << *it << std::endl;
+
+  F4timers::init();
+  F4timers::timeTotal.start();
+
+  vector<polymat<T> > G;
+  vector<polymat<T> > new_pols, new_pols2;
+  list<critpair<T> > vCrit;
+  list<critpair<T> > vCrit_step;
+
+  vector<matrixF4<T>*> vMat;
+
+  list<polymat<T> > global_reductors;
+  vector<pair<monomial, polymat<T> > > reductors;
+
+  list<polynomial<T> > _G;
+
+  for (unsigned i = 0; i < sys.size(); ++i)
+    _G.push_back(sys[i]);
+
+  list<polymat<T> > Ginit;
+  interreduce(Ginit, vMat, _G);
+  _G.clear();
+
+  for (auto& it : std::as_const(Ginit)) {
+    gebmol_gb_insert(vCrit, G, it);
+    global_reductors.push_back(it);
+  }
+  vCrit.sort();
+  unsigned countGB = sys.size();
+  unsigned step_degree;
+
+  log(LOG_INFO) << "initialized" << endl;
+
+  while (!vCrit.empty()) {
+    log(LOG_INFO) << "------------- new step ------------ " << endl;
+    log(LOG_INFO) << "total number critpairs: " << vCrit.size() << endl;
+
+    F4timers::timepretreatment.start();
+
+    step_degree = F4_critpair_selection(vCrit_step, vCrit);
+    if (vCrit_step.empty())
+      continue;
+
+    log(LOG_INFO) << "step degree: " << step_degree << endl;
+
+    F4timers::timepretreatment.stop();
+    F4timers::timeSelectReductors.start();
+    F4_reductors_selection(vCrit_step, global_reductors, reductors);
+    F4timers::timeSelectReductors.stop();
+
+    list<pair<monomial, polymat<T> > > lmat;
+    for (auto& it : as_const(reductors)) {
+      lmat.push_back(it);
+    }
+
+    log(LOG_INFO) << "step critpairs: " << vCrit_step.size() << endl;
+    log(LOG_INFO) << "step reductors: " << reductors.size() << endl;
+    reductors.clear();
+    reductors.shrink_to_fit();
+
+    matrixF4<T> *M = new matrixF4<T>(lmat, vCrit_step);
+    vCrit_step.clear();
+    vMat.push_back(M);
+    log(LOG_INFO) << "constructed matrix of size " << M->get_rsize() << "x";
+    log(LOG_INFO) << M->get_csize() << "; " << "Asize: ";
+    log(LOG_INFO) << M->get_Asize() << "; " << endl;
+    log(LOG_INFO) << "Computing row echelon" << endl;
+    M->row_echelon();
+
+    log(LOG_INFO) << "constructed matrix of size " << M->get_rsize() << "x";
+    log(LOG_INFO) << M->get_csize() << "; " << "Asize: ";
+    log(LOG_INFO) << M->get_Asize() << "; " << endl;
+    log(LOG_INFO) << endl;
+    log(LOG_INFO) << "reconstructing polynomials: " << endl;
+    polymat<T>::construct_polymats(new_pols, M);
+    new_pols2.clear();
+    for (unsigned _i = 0 ; _i < new_pols.size(); ++_i) {
+      global_reductors.push_back(new_pols[_i]);
+      if (_i >= M->get_Asize())
+        new_pols2.push_back(new_pols[_i]);
+    }
+
+    new_pols.clear();
+    new_pols.shrink_to_fit();
+    log(LOG_INFO) << "new polynomials: " << new_pols2.size() << endl;
+
+    F4timers::timeposttreatment.start();
+    log(LOG_INFO) << "updating critical pairs" << endl;
+    for (unsigned ii = 0; ii < new_pols2.size(); ++ii) {
+      countGB++;
+      gebmol_gb_insert(vCrit, G, new_pols2[ii]);
+    }
+    vCrit.sort();  // TODO(pj): improve efficiency
+    F4timers::timeposttreatment.stop();
+    new_pols2.clear();
+    new_pols2.shrink_to_fit();
+  }
+  log(LOG_INFO) << "size GB: " << G.size() << endl;
+
+  F4timers::timeTotal.stop();
+
+  log(LOG_INFO) << F4timers::timeTotal << endl;
+  log(LOG_INFO) << F4timers::timeLinalg << endl;
+  log(LOG_INFO) << F4timers::timepretreatment << endl;
+  log(LOG_INFO) << F4timers::timeSelectReductors << endl;
+  log(LOG_INFO) << F4timers::timeposttreatment << endl;
+  log(LOG_INFO) << F4timers::timeMatrixConstruct << endl;
+  log(LOG_INFO) << F4timers::timeMatrixCompression << endl;
+
+  log(LOG_INFO) << "------------ cleaning reductors -----------" << endl;
+  global_reductors.clear();
+  polynomial<T> tmp_P;
+  res.clear();
+  for (unsigned i = 0; i < G.size(); ++i) {
+    G[i].to_polynomial(tmp_P);
+    res.push_back(tmp_P);
+  }
+  std::sort(res.begin(), res.end());
+  std::reverse(res.begin(), res.end());
+  log(LOG_DEBUG) << "Is minimal: " << gb_is_minimal<T>(res) << endl;
+}
 }  // namespace tinygb
 
-#include "algos.tpp"
+
 #endif  // ALGOS_H_

src/algos.tpp

-#include <algorithm>
-#include <list>
-#include <map>
-#include <set>
-#include <utility>
-#include <vector>
-#include "./critpair.h"
-#include "./linalg.h"
-#include "./timer.h"
-
-namespace tinygb {
-
-// test whether a GB is minimal
-template <class T>
-bool gb_is_minimal(const vector<polynomial<T> > &GB) {
-  for (unsigned i = 0; i < GB.size(); ++i)
-    for (unsigned j = 0; j < GB.size(); ++j)
-      if (i != j && GB[i].LM().isDivisibleBy(GB[j].LM()))
-        return false;
-  return true;
-}
-
-// computes a minimal GB from a GB
-template <class T>
-bool GB_minimization(vector<polymat<T> > &G) {
-  bool flag = false;
-  for (unsigned i = 0; i < G.size() - 1; ++i) {
-    if (G[i].get_LM().isDivisibleBy(G[G.size()-1].get_LM())) {
-      G.erase(G.begin() + i);
-      --i;
-      flag = true;
-    }
-  }
-  return flag;
-}
-
-// Gebauer/Moller installation of Buchberger's criteria (see description
-// Vanessa Vitse's thesis)
-template <class T>
-void gebmol_gb_insert(list<critpair<T> > &vCrit, vector<polymat<T> > &G,
-    const polymat<T> &f) {
-  for (auto it = vCrit.begin(); it != vCrit.end(); ++it) {
-    if (it->LCM() != it->get_S1().second.get_LM().LCM(f.get_LM()) &&
-        it->LCM() != it->get_S2().second.get_LM().LCM(f.get_LM()) &&
-        it->LCM().isDivisibleBy(it->get_S1().second.get_LM().LCM(f.get_LM()))
-        &&
-        it->LCM().isDivisibleBy(it->get_S2().second.get_LM().LCM(f.get_LM()))) {
-      it = vCrit.erase(it);
-      it--;
-    }
-  }
-  vector<critpair<T> > S0, S1, S2;
-  for (unsigned i = 0; i < G.size(); ++i) {
-    if (G[i].get_LM().LCM(f.get_LM()) == G[i].get_LM() * f.get_LM())
-      S0.push_back(critpair<T>(f, G[i]));
-    else
-      S1.push_back(critpair<T>(f, G[i]));
-  }
-  while (!S1.empty()) {
-    critpair<T> cp1 = S1[S1.size() - 1];
-    S1.erase(S1.end()-1);
-    bool flag = false;
-    monomial tmp_lcm = cp1.LCM();
-    for (auto& it : std::as_const(S0)) {
-      if (tmp_lcm.isDivisibleBy(it.LCM())) {
-        flag = true;
-        break;
-      }
-    }
-    if (!flag) {
-      for (auto& it : std::as_const(S1)) {
-        if (tmp_lcm.isDivisibleBy(it.LCM())) {
-          flag = true;
-          break;
-        }
-      }
-    }
-    if (!flag) {
-      for (auto& it : std::as_const(S2)) {
-        if (tmp_lcm.isDivisibleBy(it.LCM())) {
-          flag = true;
-          break;
-        }
-      }
-    }
-    if (!flag) {
-      S2.push_back(cp1);
-      vCrit.push_back(cp1);
-    }
-  }
-  if (f.istopred(G))
-    return;
-  G.push_back(f);
-  GB_minimization(G);
-}
-
-// selects critical pairs to reduce in a F4 step from the current list of
-// critical pairs
-template <class T>
-unsigned F4_critpair_selection(list<critpair<T> > &vCrit_step,
-    list<critpair<T> > &vCrit) {
-  unsigned stepdegree;
-  vCrit_step.clear();
-
-#ifdef BILIN
-  // TODO(pj): make const
-  // TODO(pj): improve
-  stepdegree = vCrit.front().degree_bilin();
-  for (typename list<critpair<T> >::iterator it = vCrit.begin();
-      it != vCrit.end(); ++it)
-    if (it->degree_bilin() < stepdegree)
-      stepdegree = it->degree_bilin();
-#else
-  stepdegree = vCrit.front().degree();
-  for (auto& it : std::as_const(vCrit))
-    if (it.degree() < stepdegree)
-      stepdegree = it.degree();
-#endif
-
-  for (typename list<critpair<T> >::iterator it = vCrit.begin();
-      it != vCrit.end(); ++it)
-#ifdef BILIN
-    if (it->degree_bilin() == stepdegree) {
-      vCrit_step.push_back(*it);
-      it = vCrit.erase(it);
-      it--;
-    }
-#else
-  if (it->degree() == stepdegree) {
-    vCrit_step.push_back(*it);
-    it = vCrit.erase(it);
-    it--;
-  }
-#endif