rewriting polynomial.h

parent aa6475f0
......@@ -28,7 +28,7 @@ void launchF4(istream &in, ostream &out) {
F4<T>(sys, G);
for (typename std::vector<Polynomial <T>>::const_iterator it = G.begin(); \
it != G.end(); ++it) {
it->print(out);
it->Print(out);
out << std::endl;
}
}
......
......@@ -7,7 +7,10 @@ fields/gfp_05.inl fields/gfp_1.inl fields/gfp_2.inl fields/gf2_10.tpp \
fields/gf2_30.tpp fields/gfp_05.tpp fields/gfp_1.tpp fields/gfp_2.tpp \
annotator.h F4_launcher.h F4_launcher.tpp
tinygb_CXXFLAGS = ${SP_CXXFLAGS} -I../include -I../include/mpfq -O3 -funroll-loops -std=c++17 -DNDEBUG \
# tinygb_CXXFLAGS = ${SP_CXXFLAGS} -I../include -I../include/mpfq -O3 -funroll-loops -std=c++17 -DNDEBUG \
# -march=native -mtune=native
tinygb_CXXFLAGS = ${SP_CXXFLAGS} -I../include -I../include/mpfq -O3 -funroll-loops -std=c++17 -g \
-march=native -mtune=native
tinygb_LDADD = ../lib/mpfq-1.2/libmpfq_p_0_5.a ../lib/mpfq-1.2/libmpfq_pm_1.a \
......
......@@ -26,7 +26,7 @@ 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()))
if (i != j && GB[i].LeadingMonomial().IsDivisibleBy(GB[j].LeadingMonomial()))
return false;
return true;
}
......@@ -36,7 +36,7 @@ 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())) {
if (G[i].get_LM().IsDivisibleBy(G[G.size()-1].get_LM())) {
G.erase(G.begin() + i);
--i;
flag = true;
......@@ -53,9 +53,9 @@ void gebmol_gb_insert(list<critpair<T> > &vCrit, vector<polymat<T> > &G,
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_S1().second.get_LM().LCM(f.get_LM()))
&&
it->LCM().isDivisibleBy(it->get_S2().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--;
}
......@@ -73,14 +73,14 @@ void gebmol_gb_insert(list<critpair<T> > &vCrit, vector<polymat<T> > &G,
bool flag = false;
Monomial tmp_lcm = cp1.LCM();
for (auto& it : std::as_const(S0)) {
if (tmp_lcm.isDivisibleBy(it.LCM())) {
if (tmp_lcm.IsDivisibleBy(it.LCM())) {
flag = true;
break;
}
}
if (!flag) {
for (auto& it : std::as_const(S1)) {
if (tmp_lcm.isDivisibleBy(it.LCM())) {
if (tmp_lcm.IsDivisibleBy(it.LCM())) {
flag = true;
break;
}
......@@ -88,7 +88,7 @@ void gebmol_gb_insert(list<critpair<T> > &vCrit, vector<polymat<T> > &G,
}
if (!flag) {
for (auto& it : std::as_const(S2)) {
if (tmp_lcm.isDivisibleBy(it.LCM())) {
if (tmp_lcm.IsDivisibleBy(it.LCM())) {
flag = true;
break;
}
......@@ -116,28 +116,28 @@ unsigned F4_critpair_selection(list<critpair<T> > &vCrit_step,
#ifdef BILIN
// TODO(pj): make const
// TODO(pj): improve
stepdegree = vCrit.front().degree_bilin();
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();
if (it->Degree_bilin() < stepdegree)
stepdegree = it->Degree_bilin();
#else
stepdegree = vCrit.front().degree();
stepdegree = vCrit.front().Degree();
for (auto& it : std::as_const(vCrit))
if (it.degree() < stepdegree)
stepdegree = it.degree();
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) {
if (it->Degree_bilin() == stepdegree) {
vCrit_step.push_back(*it);
it = vCrit.erase(it);
it--;
}
#else
if (it->degree() == stepdegree) {
if (it->Degree() == stepdegree) {
vCrit_step.push_back(*it);
it = vCrit.erase(it);
it--;
......@@ -152,7 +152,7 @@ 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();
it.DivideByLeadingCoefficient();
vector<polymat<T> > res;
matrixF4<T> *M = new matrixF4<T>(_P);
......@@ -188,7 +188,7 @@ void F4_reductors_selection(const list<critpair<T> > &vCrit_step,
for (auto it = std::as_const(globred).begin();
it != std::as_const(globred).end();
++it) {
if (m.isDivisibleBy(it->get_LM())) {
if (m.IsDivisibleBy(it->get_LM())) {
if (!flag) {
flag = true;
itred = it;
......
......@@ -40,7 +40,7 @@ class critpair {
}
#ifdef BILIN
unsigned degree_bilin() {
unsigned Degree_bilin() {
if (degreebilin == 0) {
set<Monomial> setmon;
S1.second.merge_sets(setmon, S1.first);
......@@ -48,13 +48,13 @@ class critpair {
degreebilin = 0;
for (set<Monomial>::const_iterator it = setmon.begin();
it != setmon.end(); ++it)
if (it->degree_bilin() > degreebilin)
degreebilin = it->degree_bilin();
if (it->Degree_bilin() > degreebilin)
degreebilin = it->Degree_bilin();
}
return degreebilin;
}
#else
unsigned degree() const {return S1.second.degree()+S1.first.degree();}
unsigned Degree() const {return S1.second.Degree()+S1.first.Degree();}
#endif
Polynomial<T> eval() const;
......@@ -77,13 +77,13 @@ critpair<T>::critpair(const polymat<T> &p1, const polymat<T> &p2) {
#ifdef BILIN
template <class T>
bool critpair<T>::operator<(critpair<T> &cp2) {
if (degree_bilin() < cp2.degree_bilin()) return true;
if (Degree_bilin() < cp2.Degree_bilin()) return true;
return false;
}
#else
template <class T>
bool critpair<T>::operator<(critpair<T> &cp2) {
if (degree() < cp2.degree()) return true;
if (Degree() < cp2.Degree()) return true;
return false;
}
#endif
......@@ -108,7 +108,7 @@ Polynomial<T> critpair<T>::eval() const {
S2.second.to_polynomial(R);
R = R*S2.first;
R -= tmp;
R.normalize();
R.Normalize();
return R;
}
} // namespace tinygb
......
......@@ -34,10 +34,10 @@ class Monomial {
Monomial();
~Monomial();
inline expo_int degree() const;
inline expo_int Degree() const;
#ifdef BILIN
expo_int degree_bilin() const;
expo_int Degree_bilin() const;
#endif
Monomial& operator=(const Monomial& m2);
......@@ -52,7 +52,7 @@ class Monomial {
Monomial LCM(const Monomial&) const;
std::pair<Monomial, Monomial> cofact(const Monomial&) const;
inline bool isDivisibleBy(const Monomial&) const;
inline bool IsDivisibleBy(const Monomial&) const;
void print(std::ostream&) const;
......@@ -105,7 +105,7 @@ exp = new expo_int[NBVAR];
return *this;
}
inline expo_int Monomial::degree() const {
inline expo_int Monomial::Degree() const {
expo_int R = 0;
for (unsigned i = 0; i < NBVAR; ++i)
R += exp[i];
......@@ -116,7 +116,7 @@ inline expo_int Monomial::degree() const {
// defines a "bilinear degree" used for critical pair selection for bilinear
// systems
expo_int Monomial::degree_bilin() const {
expo_int Monomial::Degree_bilin() const {
expo_int Rx = 0, Ry = 0;
for (unsigned i = 0; i < Nx; ++i)
Rx += exp[i];
......@@ -139,7 +139,7 @@ void Monomial::operator*=(const Monomial &m2) {
}
Monomial Monomial::operator/(const Monomial &m2) const {
assert(isDivisibleBy(m2));
assert(IsDivisibleBy(m2));
Monomial R;
for (unsigned i = 0; i < NBVAR; ++i)
R.exp[i] = exp[i] - m2.exp[i];
......@@ -189,8 +189,8 @@ inline bool Monomial::operator<(const Monomial &m2) const {
# else
inline bool Monomial::operator<(const Monomial &m2) const {
expo_int d1 = degree();
expo_int d2 = m2.degree();
expo_int d1 = Degree();
expo_int d2 = m2.Degree();
if (d1 < d2) return true;
if (d1 > d2) return false;
int i = NBVAR-1;
......@@ -255,7 +255,7 @@ std::pair<Monomial, Monomial> Monomial::cofact(const Monomial &m2) const {
return std::pair<Monomial, Monomial>(R1, R2);
}
inline bool Monomial::isDivisibleBy(const Monomial &m2) const {
inline bool Monomial::IsDivisibleBy(const Monomial &m2) const {
for (unsigned i = 0; i < NBVAR; ++i)
if (exp[i] < m2.exp[i])
return false;
......
......@@ -37,7 +37,7 @@ vector<Polynomial<T> > Parser<T>::parse() const {
p.insert(m, e);
in_stream >> s;
}
p.normalize();
p.Normalize();
res.push_back(p);
}
return res;
......
......@@ -34,7 +34,7 @@ class polymat {
polymat(const matrixF4<T> *_M, const unsigned _row) : M(_M), row(_row) {
Polynomial<T> P;
to_polynomial(P);
LM = P.LM();
LM = P.LeadingMonomial();
}
polymat(const polymat<T> &p) : M(p.M), row(p.row) {
......@@ -48,11 +48,11 @@ class polymat {
void to_polynomial(Polynomial<T>&) const;
#ifdef BILIN
unsigned degree_bilin() const;
unsigned Degree_bilin() const;
#endif
unsigned degree() const {
return LM.degree();
unsigned Degree() const {
return LM.Degree();
} // TODO(pj): improve
bool istopred(const vector<polymat<T> >&) const;
......@@ -72,20 +72,20 @@ class polymat {
#ifdef BILIN
template <class T>
unsigned polymat<T>::degree_bilin() const {
unsigned polymat<T>::Degree_bilin() const {
unsigned res = 0;
if (row < M->Asize) {
if (M->v1[row].degree_bilin() > res)
res = M->v1[row].degree_bilin();
if (M->v1[row].Degree_bilin() > res)
res = M->v1[row].Degree_bilin();
for (unsigned j = 0; j < M->B->get_csize(); ++j)
if (!(*(M->B))(row, j).is_zero())
if (M->v2[j].degree_bilin() > res)
res = M->v2[j].degree_bilin();
if (M->v2[j].Degree_bilin() > res)
res = M->v2[j].Degree_bilin();
} else {
for (unsigned j = 0; j < M->D->get_csize(); ++j)
if (!(*(M->D))(row - M->Asize, j).is_zero())
if (M->v2[j].degree_bilin() > res)
res = M->v2[j].degree_bilin();
if (M->v2[j].Degree_bilin() > res)
res = M->v2[j].Degree_bilin();
}
return res;
}
......@@ -102,8 +102,8 @@ void polymat<T>::to_polynomial(Polynomial<T> &P) const {
if (!(*(M->get_B()))(row, j).is_zero())
P.insert(M->get_v2()[j], (*(M->get_B()))(row, j));
}
P.normalize();
P.set_monic();
P.Normalize();
P.DivideByLeadingCoefficient();
} else {
assert(row < M->get_rsize());
P.insert(M->get_pivots_LM()[row-M->get_Asize()], 1);
......@@ -112,8 +112,8 @@ void polymat<T>::to_polynomial(Polynomial<T> &P) const {
if (!(*(M->get_D()))(row-M->get_Asize(), j).is_zero())
P.insert(M->get_v2()[j], (*(M->get_D()))(row-M->get_Asize(), j));
}
P.normalize();
P.set_monic();
P.Normalize();
P.DivideByLeadingCoefficient();
}
}
......@@ -179,7 +179,7 @@ void polymat<T>::merge_sets(set<Monomial> &prio, const Monomial &m) const {
template<class T>
bool polymat<T>::istopred(const vector<polymat<T> > &G) const {
for (unsigned i = 0; i < G.size(); ++i)
if (LM.isDivisibleBy(G[i].get_LM()))
if (LM.IsDivisibleBy(G[i].get_LM()))
return true;
return false;
}
......
......@@ -20,16 +20,19 @@ class Polynomial {
typedef std::pair<Monomial, T> Term;
Polynomial() {
terms_ = std::list<std::pair<Monomial, T> >();
terms_ = std::list<Term>();
}
explicit Polynomial(const std::list<Term>& terms_)
: terms_(terms_) {}
const std::list<std::pair<Monomial, T> >& terms() const {
const std::list<Term>& terms() const {
return terms_;
}
void normalize();
// Put the list of terms in standard form. Depends on the monomial ordering.
// The list is ordered from smallest to largest
// TODO(pj): leading monomial should be first
void Normalize();
void SetToZero() {
terms_.clear();
......@@ -39,51 +42,58 @@ class Polynomial {
terms_.push_back(std::pair<Monomial, T>(m, e));
}
const Monomial& LM() const {
assert(!terms_.back().second.IsZero());
// Assumes that the polynomial's list of terms is normalized
const Monomial& LeadingMonomial() const {
assert(!terms_.back().second.is_zero());
return terms_.back().first;
}
const T& LC() const {
const T& LeadingCoefficient() const {
return terms_.back().second;
}
void set_monic();
void DivideByLeadingCoefficient();
unsigned degree() const {
// TODO(pj): extend for non-degree orderings
return LM().degree();
unsigned Degree() const {
// TODO(pj): Extend for non-degree orderings.
return LeadingMonomial().Degree();
}
#ifdef BILIN
unsigned degree_bilin() const;
unsigned Degree_bilin() const;
#endif
bool topred(const Polynomial &g);
bool istopred(const Polynomial &g) const;
T coeff(const Monomial& m) const;
// Reduces leading monomial by the polynomial given.
// Returns true if leading monomial is divisible by the leading monomial of
// the argument. When this function returns false, then nothing has been
// modified.
// Assumes that terms_ is normalized.
bool TopReductionBy(const Polynomial &);
// Returns true if the leading monomial is divisible by the leading monomial
// of the argument.
// Assumes that terms_ is normalized.
bool IsTopReducibleBy(const Polynomial &) const;
Polynomial<T> operator*(const Monomial&) const;
void operator+=(const Polynomial<T>&);
void operator-=(const Polynomial<T>&);
void operator*=(const T&);
void operator/=(const T&);
void print(std::ostream& out) const;
void Print(std::ostream& out) const;
inline bool operator<(const Polynomial<T>& P) const;
friend std::ostream& operator<< <T>(std::ostream&, const Polynomial<T>&);
private:
// Adds coefficients of the two leading terms if they share the same
// monomials. Useful for adding polynomials.
void CombineLeadingMonomials();
private:
// Adds coefficients of consecutive terms if they share the same monomials.
// Useful for adding polynomials. Do not use if terms_ is not sorted.
void CombineMonomials();
bool IsEmpty() const {
return terms_.empty();
}
// Removes terms with coefficient 0
void RemoveZeroes();
void RemoveZeroTerms();
std::list<Term> terms_;
......@@ -91,153 +101,115 @@ class Polynomial {
#ifdef BILIN
template<class T>
unsigned Polynomial<T>::degree_bilin() const {
unsigned Polynomial<T>::Degree_bilin() const {
unsigned R = 0;
for (typename std::list<std::pair<Monomial, T> >::const_iterator it = \
terms_.begin(); it != terms_.end(); ++it)
if (it->first.degree_bilin() > R)
R = it->first.degree_bilin();
for (auto& it = std::as_const(terms_))
if (it->first.Degree_bilin() > R)
R = it->first.Degree_bilin();
return R;
}
#endif
template<class T>
void Polynomial<T>::RemoveZeroes() {
void Polynomial<T>::RemoveZeroTerms() {
for (auto it = terms_.begin(); it != terms_.end(); ++it) {
if (it->second == T(0)) {
while (it != terms_.end() && it->second == T(0))
it = terms_.erase(it);
continue;
}
}
}
template<class T>
void Polynomial<T>::CombineLeadingMonomials() {
void Polynomial<T>::CombineMonomials() {
typename std::list<std::pair<Monomial, T> >::iterator it2;
for (auto it = terms_.begin(); it != terms_.end();) {
it2 = it;
++it;
if (it == terms_.end()) {
for (auto it1 = terms_.begin(); it1 != terms_.end();) {
it2 = it1;
++it1;
if (it1 == terms_.end()) {
break;
} else if (it->first == it2->first) {
it->second += it2->second;
} else if (it1->first == it2->first) {
it1->second += it2->second;
terms_.erase(it2);
}
}
}
template<class T>
T Polynomial<T>::coeff(const Monomial& m) const {
for (typename std::list<std::pair<Monomial, T> >::const_iterator it = \
terms_.begin(); it != terms_.end(); ++it)
if (it->first == m)
return it->second;
return T(0);
}
template<class T>
void Polynomial<T>::set_monic() {
RemoveZeroes();
CombineLeadingMonomials();
T leading_coeff = terms_.back().second;
leading_coeff.inv();
// TODO(pj): test if 1
void Polynomial<T>::DivideByLeadingCoefficient() {
T leading_coefficient = terms_.back().second;
assert(leading_coefficient != T(0));
leading_coefficient.inv();
if (leading_coefficient == T(1)) return;
for (auto& it : terms_)
it.second *= leading_coeff;
it.second *= leading_coefficient;
}
template<class T>
void Polynomial<T>::normalize() {
void Polynomial<T>::Normalize() {
terms_.sort();
RemoveZeroes();
CombineLeadingMonomials();
CombineMonomials();
RemoveZeroTerms();
}
template<class T>
bool Polynomial<T>::topred(const Polynomial &g) {
if (g.IsEmpty() || IsEmpty())
return false;
if (!LM().isDivisibleBy(g.LM()))
return false;
Polynomial reductor = g*(LM()/g.LM());
reductor.set_monic();
reductor *= LC();
bool Polynomial<T>::TopReductionBy(const Polynomial &g) {
assert(!g.IsEmpty());
if (IsEmpty()) return false;
if (!LeadingMonomial().IsDivisibleBy(g.LeadingMonomial())) return false;
Polynomial reductor = g*(LeadingMonomial()/g.LeadingMonomial());
reductor.DivideByLeadingCoefficient();
reductor *= LeadingCoefficient();
operator-=(reductor);
return true;
}
template<class T>
bool Polynomial<T>::istopred(const Polynomial<T> &g) const {
if (g.IsEmpty() || IsEmpty()) return false;
return LM().isDivisibleBy(g.LM());
bool Polynomial<T>::IsTopReducibleBy(const Polynomial<T> &g) const {
assert(!g.IsEmpty());
if (IsEmpty()) return false;
return LeadingMonomial().IsDivisibleBy(g.LeadingMonomial());
}
template<class T>
Polynomial<T> Polynomial<T>::operator*(const Monomial &m) const {
Polynomial R;
R.terms_.insert(R.terms_.begin(), this->terms_.begin(), this->terms_.end());
for (typename std::list<std::pair<Monomial, T> >::iterator it = \
R.terms_.begin(); it != R.terms_.end(); ++it)
(it->first) *= m;
for (auto& it : R.terms_)
it.first *= m;
return R;
}
template<class T>
void Polynomial<T>::operator+=(const Polynomial<T> &p2) {
terms_.merge(p2.terms_);
RemoveZeroes();
CombineLeadingMonomials();
void Polynomial<T>::operator+=(const Polynomial<T> &P) {
std::list<Term> tmp_terms = P.terms_;
terms_.merge(tmp_terms);
CombineMonomials();
RemoveZeroTerms();
}
template<class T>
void Polynomial<T>::operator*=(const T &e) {
for (typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin(); \
it != terms_.end(); ++it)
it->second *= e;
for (auto& it : terms_)
it.second *= e;
}
template<class T>
void Polynomial<T>::operator/=(const T& e) {
for (typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin(); \
it != terms_.end(); ++it)
it->second /= e;
for (auto& it : terms_)
it.second /= e;
}
template <class T>
void Polynomial<T>::operator-=(const Polynomial<T>& p2) {
typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin();
typename std::list<std::pair<Monomial, T> >::const_iterator it2 = \
p2.terms_.begin();
typename std::list<std::pair<Monomial, T> >::const_iterator it2end = \
p2.terms_.end();
while (it2 != p2.terms_.end()) {
if (it == terms_.end()) {
terms_.insert(it, it2, it2end);
while (it != terms_.end()) {
it->second.neg(); ++it;
}