From 7cde96e48f2f526b9fbb9aa128f7a80be7329f5d Mon Sep 17 00:00:00 2001 From: Wojciech Politarczyk Date: Mon, 28 Nov 2016 09:26:41 +0100 Subject: [PATCH] Verification of perioidcity congruences is finished Verification of periodicity congruences is finished and works ok. --- .gitignore | 1 + Makefile | 4 +- algebra/Z.h | 19 +++ algebra/multivariate_laurentpoly.h | 3 - kk.cpp | 111 ++++---------- lib/map.h | 1 - lib/map_wrapper.h | 2 +- periodicity.cpp | 223 +++++++++++++++++++++++------ periodicity.h | 111 ++++++++++---- 9 files changed, 312 insertions(+), 163 deletions(-) diff --git a/.gitignore b/.gitignore index 62f5ee3..7805d31 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ #* *~ +.* *.o GPATH GRTAGS diff --git a/Makefile b/Makefile index bdee858..d79926f 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ CXX = g++ -std=c++11 INCLUDES = -I. -I/opt/local/include # OPTFLAGS = -g -OPTFLAGS = -O2 -g -pg +OPTFLAGS = -O2 -g # OPTFLAGS = -O2 -g -DNDEBUG LDFLAGS = -L/opt/local/lib @@ -120,6 +120,6 @@ realclean: clean $(LIB_OBJS): $(LIB_HEADERS) $(ALGEBRA_OBJS): $(ALGEBRA_HEADERS) $(LIB_HEADERS) $(KNOTKIT_OBJS) main.o mpimain.o kk.o: $(KNOTKIT_HEADERS) $(ALGEBRA_HEADERS) $(LIB_HEADERS) $(PERIODICITY_HEADERS) -$(PERIODICITY_OBJS) : $(PERIODICITY_HEADERS) +$(PERIODICITY_OBJS) : $(PERIODICITY_HEADERS) $(ALGEBRA_HEADERS) $(LIB_HEADERS) mpimain.o mpi_aux.o: mpi_aux.h diff --git a/algebra/Z.h b/algebra/Z.h index db61310..660d997 100644 --- a/algebra/Z.h +++ b/algebra/Z.h @@ -14,6 +14,7 @@ class Z private: std::shared_ptr impl; +#ifdef DEBUG_Z void write_state() const { std::cout << "I store the following value " << *this << "\n"; std::cout << "Number of objects pointing to the same value " << impl.use_count() << "\n"; @@ -22,6 +23,7 @@ class Z /* std::cout << "Size of std::shared_ptr " << sizeof(impl) << "\n"; */ /* std::cout << "Size of mpz_class " << sizeof(*impl) << "\n"; */ } +#endif public: Z() : impl(new mpz_class) { @@ -52,6 +54,12 @@ class Z #ifdef DEBUG_Z std::cout << "Z(const Z& z)" << "\n"; write_state(); +#endif + } + Z(copy, const Z& z) : impl(z.impl) { +#ifdef DEBUG_Z + std::cout << "Z(COPY, const Z& z)\n"; + write_state(); #endif } Z(Z&& z) : impl(std::move(z.impl)) { @@ -115,6 +123,17 @@ class Z return *impl.get() < *z.impl.get(); } + bool operator <= (const Z& z) const { + return *this < z || *this == z; + } + + bool operator > (const Z& z) const { + return !(*this <= z); + } + bool operator >= (const Z& z) const { + return !(*this < z); + } + bool is_unit () const { return *this == 1 || *this == -1; diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h index 5fe65f6..6c9d548 100644 --- a/algebra/multivariate_laurentpoly.h +++ b/algebra/multivariate_laurentpoly.h @@ -548,7 +548,6 @@ multivariate_laurentpoly invert_variable(const multivariate_laurentpoly& p return result; } -<<<<<<< HEAD template multivariate_laurentpoly multivariate_laurentpoly::evaluate(T val, unsigned index) const { @@ -569,6 +568,4 @@ multivariate_laurentpoly::evaluate(T val, unsigned index) const { return res; } -======= ->>>>>>> 297e02dc05d79c4546d04aeb720680a9c231284e #endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H diff --git a/kk.cpp b/kk.cpp index 5e9c341..4ac28ca 100644 --- a/kk.cpp +++ b/kk.cpp @@ -1,10 +1,6 @@ - #include -<<<<<<< HEAD #include #include -======= ->>>>>>> 297e02dc05d79c4546d04aeb720680a9c231284e #include #include @@ -322,6 +318,17 @@ int compute_s_inv(knot_diagram& kd) { } void check_periodicity(std::string out_file) { + if(periodicity_test == "all") { + Kh_periodicity_checker Kh_pc(kd); + Przytycki_periodicity_checker P_pc(Kh_pc.get_KhP().evaluate(-1, eval_index)); + for(auto& p : primes_list) { + std::cout << "Przytycki's criterion: " + << P_pc(p) << std::endl + << "Kh criterion: " + << Kh_pc(p) << std::endl; + } + } + else { if(period == 2 || period == 3) { std::cout << "Sorry, the criteria don't work for period " << period << "...\n"; @@ -332,73 +339,26 @@ void check_periodicity(std::string out_file) { std::cout << "For now you can only check periodicity for primes up to 31..." << "\n"; exit(EXIT_FAILURE); } - std::ofstream out(out_file); + std::ofstream out(out_file, std::ios_base::app); - if(periodicity_test == "all") { - std::cout << "I will perform both test..." << "\n"; - } - else if(periodicity_test == "Przytycki") { - switch(period) { - case 5: { - Przytycki_periodicity_checker<5> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - case 7: { - Przytycki_periodicity_checker<7> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - case 11: { - Przytycki_periodicity_checker<11> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - case 13: { - Przytycki_periodicity_checker<13> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - case 17: { - Przytycki_periodicity_checker<17> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - case 19: { - Przytycki_periodicity_checker<19> pcc; - if(out_file.size() != 0) - out << pcc(compute_jones(kd)) << std::endl; - else - std::cout << pcc(compute_jones(kd)) << std::endl; - break; - } - } + if(periodicity_test == "Przytycki") { + Przytycki_periodicity_checker P_pc(compute_jones(kd)); + if(out_file.size() != 0) + out << P_pc(period) << std::endl; + else + std::cout << P_pc(period) << std::endl; } else if(periodicity_test == "Kh") { - Kh_periodicity_checker pc(kd); + Kh_periodicity_checker Kh_pc(kd); if(out_file.size() != 0) - out << pc(period) << std::endl; + out << Kh_pc(period) << std::endl; else - std::cout << pc(period) << std::endl; + std::cout << Kh_pc(period) << std::endl; } else { - std::cout << "Sorry, I don't recognize this option..." << "\n"; - exit(EXIT_FAILURE); + std::cout << "Sorry, I don't recognize this option..." << "\n"; + exit(EXIT_FAILURE); + } } } @@ -726,7 +686,7 @@ main (int argc, char **argv) else if(!strcmp (argv[i], "-p")) { i++; if(i == argc) { - fprintf (stderr, "error: missing argument to option `-o'\n"); + fprintf (stderr, "error: missing argument to option `-p'\n"); exit (EXIT_FAILURE); } period = std::stoi(argv[i]); @@ -734,7 +694,7 @@ main (int argc, char **argv) else if (!strcmp(argv[i], "-t")) { i++; if(i == argc) { - fprintf (stderr, "error: missing argument to option `-o'\n"); + fprintf (stderr, "error: missing argument to option `-t'\n"); exit (EXIT_FAILURE); } periodicity_test = argv[i]; @@ -844,25 +804,6 @@ main (int argc, char **argv) << ") of " << knot << " = " << std::endl << khp << std::endl; } - else if(!strcmp(invariant, "periodicity_congruence")) { - if(!strcmp(field, "Z2")) { - // first we check whether period is a prime - if(period == 2 || period == 3) { - std::cerr << "The criterion does not work for period = " << period << "\n"; - exit(EXIT_FAILURE); - } - auto result = find(primes.begin(), primes.end(), period); - if(result == primes.end()) { - std::cerr << "For now it is possible to check periodicity for primes up to 31" << "\n"; - exit(EXIT_FAILURE); - } - check_periodicity_criterion(kd,period); - } - else { - std::cerr << "error: for now this function is only defined for Z2 coefficients..." << "\n"; - exit(EXIT_FAILURE); - } - } else { if (!strcmp (field, "Z2")) diff --git a/lib/map.h b/lib/map.h index 859829a..ddc1c9f 100644 --- a/lib/map.h +++ b/lib/map.h @@ -1,4 +1,3 @@ - /* wrapper for std::map */ template diff --git a/lib/map_wrapper.h b/lib/map_wrapper.h index e263412..a8d51cc 100644 --- a/lib/map_wrapper.h +++ b/lib/map_wrapper.h @@ -64,7 +64,7 @@ class map_wrapper /* returns the pair associated to the largest key */ pair tail () const { - typename M::const_reverse_iterator i = impl->t.rbegin (); + typename M::const_reverse_iterator i = impl->t.rbegin(); assert (i != impl->t.rend ()); return pair (i->first, i->second); } diff --git a/periodicity.cpp b/periodicity.cpp index 6d98527..085324c 100644 --- a/periodicity.cpp +++ b/periodicity.cpp @@ -1,6 +1,150 @@ #include #include +bool Przytycki_periodicity_checker::check(int period) const { + switch(period) { + case 5: { + periodic_congruence_checker> pcc(5); + return pcc(jones_pol); + } + case 7: { + periodic_congruence_checker> pcc(7); + return pcc(jones_pol); + } + case 11: { + periodic_congruence_checker> pcc(11); + return pcc(jones_pol); + } + case 13: { + periodic_congruence_checker> pcc(13); + return pcc(jones_pol); + } + case 17: { + periodic_congruence_checker> pcc(17); + return pcc(jones_pol); + } + case 19: { + periodic_congruence_checker> pcc(19); + return pcc(jones_pol); + } + } +} + +std::string Przytycki_periodicity_checker::operator () (int period) const { + std::ostringstream res; + res << knot << ": period = " << period << ": " + << (check(period) ? "Maybe" : "No"); + return res.str(); +} + +template +polynomial_iterator::polynomial_iterator(const multivariate_laurentpoly& pol, + start_pos sp) { + for(typename map::const_iter i = pol.coeffs; i; ++i) { + monomials.push_back(i.key()); + bounds.push_back(i.val()); + current_pos.push_back(Z(0)); + } + + if(sp == start_pos::begin) { + level = 0; + } + else { + level = bounds.size(); + } +#ifndef NDEBUG + check_current_pos(); +#endif +} + +#ifndef NDEBUG +template +void polynomial_iterator::check_current_pos() { + assert(bounds.size() == monomials.size()); + assert(bounds.size() == current_pos.size()); + assert(level <= current_pos.size()); + for(unsigned i = 0; i < current_pos.size(); i++) { + if(i < level) { + assert((current_pos[i] <= bounds[i]) && + Z(0) <= (current_pos[i])); + } + else if(i == level) { + if(level > 0) + assert(current_pos[i] <= bounds[i] && Z(0) < current_pos[i]); + else + assert((current_pos[i] <= bounds[i]) && Z(0) <= (current_pos[i])); + } + else + assert(current_pos[i] == Z(0)); + } +} +#endif // NDEBUG + +template +polynomial_iterator& +polynomial_iterator::operator ++ () { +#ifndef NDEBUG + check_current_pos(); +#endif + if(level == monomials.size()) + return *this; + + unsigned i = 0; + while(i <= level) { +#ifndef NDEBUG + check_current_pos(); +#endif + if(current_pos[i] < bounds[i]) { + current_pos[i] += 1; + break; + } + else { + if(i == level) { + if(level < monomials.size() - 1) { + current_pos[i] = 0; + current_pos[i+1] += 1; + level++; + break; + } + else { + level++; + break; + } + } + else { + current_pos[i] = 0; + } + } + i++; + } + return *this; +} + +template +multivariate_laurentpoly polynomial_iterator::operator *() const { + polynomial res; + for(unsigned i = 0; i <= level; i++) { + res += polynomial(current_pos[i], monomials[i]); + } + return res; +} + +template +std::string polynomial_iterator::write_self() const { + std::ostringstream res; + res << "level = " << level << std::endl + << "monomials:" << std::endl; + for(auto& mon : monomials) + res << mon << std::endl; + res << "bounds: " << std::endl; + for(auto& b : bounds) + res << b << std::endl; + res << "current_pos: " << std::endl; + for(auto& pos : current_pos) + res << pos << std::endl; + return res.str(); +} + void Kh_periodicity_checker::compute_knot_polynomials(knot_diagram& kd) { unsigned m = kd.num_components (); @@ -57,8 +201,14 @@ void Kh_periodicity_checker::compute_quot() { quot += polynomial(m1.second, m1.first); polynomial p = polynomial(m1.second, m1.first) * mul; diff -= p; - m1 = diff.tail(); + if(diff != 0) + m1 = diff.tail(); + else break; } + if(diff != 0) + m = diff.head(); + else + break; } quot += polynomial(m.second, m.first); polynomial p = polynomial(m.second, m.first) * mul; @@ -91,33 +241,6 @@ Kh_periodicity_checker::compute_quotient_and_remainder(const polynomial& quot, return std::make_pair(quotient, remainder); } -std::list> -Kh_periodicity_checker::generate_candidates(const polynomial &q) const { - std::list result; - Z size = 0; - map::const_iter i = q.coeffs; - - for(int j = 0; Z(j) <= i.val(); j++) { - result.push_back(polynomial(Z(j), i.key())); - size += 1; - } - - i++; - - for( ; i; i++) { - for(int j = 1; Z(j) <= i.val(); j++) { - std::list temp_list; - for_each(result.begin(), result.end(), - [&result, &temp_list, j, i](const polynomial& p){ - temp_list.push_back(p + polynomial(Z(j), i.key())); - }); - result.splice(result.end(), temp_list); - size += 1; - } - } - return result; -} - bool Kh_periodicity_checker::check(const polynomial& q, const polynomial& r, int period) const { @@ -126,26 +249,34 @@ bool Kh_periodicity_checker::check(const polynomial& q, if(q == 0) { return pcc(t.evaluate(-1,1)); } - else { - // generate all polynomials - if(verbose) { - std::cout << "Generating all candidates..." << std::endl; + if(verbose) + std::cout << "Checking congruences..."; + polynomial_iterator pi(q); + polynomial_iterator pi_end(q, polynomial_iterator::start_pos::end); + Z count = pi.get_count(); + if(verbose) + std::cout << count << " candidates..." << std::endl; + Z step = 1; + if(count >= 32) + step = count / 32; + Z c = 0; + while(pi != pi_end) { + //std::cout << "pi: " << std::endl << pi; + polynomial temp = t + polynomial(period - 1) * mul * (*pi); + if(pcc(temp.evaluate(-1,1))) { + std::cout << "Candidates:" << std::endl + << "EKhP_1 = " << temp << std::endl + << "EKhP_" << (period - 1) << " = " + << (polynomial(period - 1) * mul *(r - *pi)) + << std::endl; + return true; } - std::list candidates = generate_candidates(q); - // Z i = 0; - // for(auto& l : candidates) { - // i += 1; - // std::cout << i << ": " << l << std::endl; - // } - // and check each one of them - if(verbose) - std::cout << "Checking congruences..." << std::endl; - polynomial m = mul; - return any_of(candidates.begin(), candidates.end(), - [&pcc, &t, &m, &period](const polynomial& p){ - return pcc((t + Z(period - 1) * m * p).evaluate(-1,1)); - }); + ++pi; + c += 1; + if(verbose && c % step == 0) + std::cout << c << "/" << count << "..." << std::endl; } + return false; } std::string Kh_periodicity_checker::operator () (int period) const { diff --git a/periodicity.h b/periodicity.h index 616b5bb..6c83e57 100644 --- a/periodicity.h +++ b/periodicity.h @@ -5,9 +5,6 @@ #include #include #include -#include -#include -#include extern bool verbose; extern const char* knot; @@ -16,15 +13,12 @@ extern std::string periodicity_test; const std::vector primes_list = {5, 7, 11, 13, 17, 19}; -enum class Test_type { - Przytycki_criterion, - BKP_criterion, - all -}; +const unsigned eval_index = 1; +const unsigned invert_index = 2; template class periodic_congruence_checker { -public: +protected: using polynomial = multivariate_laurentpoly; using monomial = multivariate_laurent_monomial; @@ -32,19 +26,20 @@ public: unsigned index; polynomial prepare_polynomial(const polynomial& pol) const { - return (pol - invert_variable(pol, index)); + polynomial inv = invert_variable(pol, index); + return pol - inv; } bool reduce(const polynomial& pol) const; public: periodic_congruence_checker(int pprime = 5, - unsigned ind = 2) : + unsigned ind = invert_index) : prime(pprime), index(ind) {} - ~periodic_congruence_checker() {}; + virtual ~periodic_congruence_checker() {}; bool operator() (const polynomial& pol) const { return reduce(prepare_polynomial(pol)); @@ -61,32 +56,90 @@ bool periodic_congruence_checker::reduce(const multivariate_laurentpoly& p monomial mon = monomial(VARIABLE, index, c); res += polynomial(i.val(), mon); } - if(verbose) - std::cout << "reduced = " << res << "\n"; + // if(verbose) + // std::cout << "reduced = " << res << "\n"; return res == 0; } -template class Przytycki_periodicity_checker { - using polynomial = multivariate_laurentpoly>; + using polynomial = multivariate_laurentpoly; using monomial = multivariate_laurent_monomial; - periodic_congruence_checker> cong_checker; + polynomial jones_pol; + + bool check(int period) const; public: - Przytycki_periodicity_checker() : cong_checker(p) {} + Przytycki_periodicity_checker(polynomial j) : jones_pol(j) {} ~Przytycki_periodicity_checker() {} - std::string operator() (const polynomial& pol) const { - std::ostringstream out; - out << knot << ": period = " << p - << ": " - << (cong_checker(pol) ? "Maybe" : "No"); - return out.str(); + std::string operator() (int period) const; +}; + +template +class polynomial_iterator { + using polynomial = multivariate_laurentpoly; + using monomial = multivariate_laurent_monomial; + + std::vector monomials; + std::vector bounds; + std::vector current_pos; + unsigned level; + + void check_current_pos(); + +public: + enum class start_pos { begin, end }; + + polynomial_iterator(const polynomial& init, + start_pos sp = start_pos::begin); + polynomial_iterator(const polynomial_iterator& pi) =default; + polynomial_iterator(polynomial_iterator&& pi) =default; + + ~polynomial_iterator() {} + + polynomial_iterator& operator = (const polynomial_iterator& pi) =default; + polynomial_iterator& operator = (polynomial_iterator&& pi) =default; + + polynomial_iterator& operator ++(); + + bool operator == (const polynomial_iterator& pi) const { + if(level == monomials.size() || pi.level == pi.monomials.size()) { + return level == pi.level && + monomials == pi.monomials && + bounds == pi.bounds; + } + else { + return level == pi.level && + bounds == pi.bounds && + monomials == pi.monomials && + current_pos == pi.current_pos; + } + } + bool operator != (const polynomial_iterator& pi) const { + return !(*this == pi); + } + polynomial operator*() const; + + T get_count() const { + Z res = 1; + for(auto& v : bounds) + res *= (v + 1); + return res; + } + + std::string write_self() const; + friend inline std::ostream& operator << (std::ostream& os, const polynomial_iterator& pi) { + return os << pi.write_self(); } }; +template +std::ostream& operator << (std::ostream& os, const polynomial_iterator& pi) { + return os << *pi; +} + class Kh_periodicity_checker { using polynomial = multivariate_laurentpoly; using monomial = multivariate_laurent_monomial; @@ -101,7 +154,7 @@ class Kh_periodicity_checker { void compute_quot(); std::pair compute_quotient_and_remainder(const polynomial& p, int period) const; - std::list generate_candidates(const polynomial& q) const; + // std::list generate_candidates(const polynomial& q) const; bool check(const polynomial& q, const polynomial& r, int period) const; public: @@ -119,6 +172,14 @@ class Kh_periodicity_checker { ~Kh_periodicity_checker() {} std::string operator () (int period) const; + + polynomial get_KhP() const { + return khp; + } + + polynomial get_LeeP() const { + return leep; + } }; #endif // _KNOTKIT_PERIODICITY_H