OKlibrary  0.2.1.6
GreenTao.hpp
Go to the documentation of this file.
00001 // Oliver Kullmann, 17.10.2009 (Swansea)
00002 /* Copyright 2009 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007 
00013 #ifndef GREENTAOHYPERGRAPH_mmnVre09h
00014 #define GREENTAOHYPERGRAPH_mmnVre09h
00015 
00016 #include <vector>
00017 #include <cassert>
00018 #include <stdexcept>
00019 #include <limits>
00020 #include <utility>
00021 
00022 #include <gmpxx.h>
00023 
00024 namespace OKlib {
00025   namespace Combinatorics {
00026     namespace Hypergraphs {
00027       namespace Generators {
00028 
00037         template <typename UInt = unsigned long>
00038         struct First_prime_numbers {
00039           typedef UInt uint_type;
00040           typedef std::vector<uint_type> range_type;
00041           range_type operator()(uint_type n) const {
00042             range_type primes;
00043             primes.reserve(n);
00044             const uint_type max = std::numeric_limits<uint_type>::max();
00045             for (mpz_class p = 2; n > 0; --n) {
00046               if (p > max)
00047                 throw std::overflow_error("ERROR[First_prime_numbers]: n too big.");
00048               primes.push_back(p.get_ui());
00049               mpz_nextprime(p.get_mpz_t(), p.get_mpz_t());
00050             }
00051             return primes;
00052           }
00053         };
00054         template <typename UInt>
00055         typename First_prime_numbers<UInt>::range_type first_prime_numbers(const UInt n) {
00056           return First_prime_numbers<UInt>()(n);
00057         }
00058 
00093         template <typename UInt = unsigned long>
00094         struct GreenTao {
00095 
00096           typedef UInt vertex_type;
00097           typedef std::vector<vertex_type> hyperedge_type;
00098           typedef std::vector<hyperedge_type> set_system_type;
00099           typedef typename set_system_type::size_type size_type;
00100 
00101           const vertex_type k; // length of arithmetic progressions
00102           const bool prime_k;
00103 
00104           static const unsigned int max_k = 19; // for k=20 the first
00105           // arithmetic progressions appears for n=22,009,064,470.
00106           static const unsigned int rounds_compositeness_test = 10;
00107 
00108           GreenTao(const vertex_type k, const vertex_type n)
00109             : k(k), prime_k(prime(k)), vertex_set_(first_prime_numbers(n)), n(n), max_prime((n==0)?0:vertex_set_.back()) {
00110             // C++0X: the following assert should become compile-time
00111             assert(std::numeric_limits<vertex_type>::max() >= 4294967295UL);
00112             assert(k <= max_k);
00113           }
00114           GreenTao(const vertex_type k, const vertex_type n, const hyperedge_type& V)
00115             : k(k), prime_k(prime(k)), vertex_set_(V), n(n), max_prime((n==0)?0:vertex_set_.back()) {
00116             // C++0X: the following assert should become compile-time
00117             assert(std::numeric_limits<vertex_type>::max() >= 4294967295UL);
00118             assert(k <= max_k);
00119             assert(vertex_set_.size() == n);
00120           }
00121 
00122           size_type nver() const { return n; }
00123           size_type max_index() const { return max_prime; }
00124           size_type nhyp() const; // XXX
00125 
00126           const hyperedge_type& vertex_set() const { return vertex_set_; }
00127 
00128           set_system_type hyperedge_set() const {
00129             set_system_type result;
00130             if (k==0) { result.push_back(hyperedge_type()); return result; }
00131             if (n==0) return result;
00132             assert(not vertex_set_.empty());
00133             typedef typename hyperedge_type::const_iterator v_it_type;
00134             const v_it_type vbegin(vertex_set_.begin());
00135             const v_it_type vend(vertex_set_.end());
00136             if (k==1) {
00137               for (v_it_type i = vbegin; i != vend; ++i)
00138                 result.push_back(hyperedge_type(1,*i));
00139               return result;
00140             }
00141             if (k==2) {
00142               for (v_it_type j = vbegin; j != vend; ++j)
00143                 for (v_it_type i = vbegin; i != j; ++i) {
00144                   // C++0X: result.push_back(hyperedge_type {*i,*j})
00145                   hyperedge_type H;
00146                   H.push_back(*i); H.push_back(*j);
00147                   result.push_back(H);
00148                 }
00149               return result;
00150             }
00151             if (n <= k) return result;
00152             assert(n >= 4);
00153             assert(k >= 3);
00154             assert(n > k);
00155             std::vector<bool> primes_table(vertex_set_.back()+1,false);
00156               for (v_it_type i = vbegin; i != vend; ++i)
00157                 primes_table[*i] = true;
00158             vertex_type prd_primes = 1;
00159             for (v_it_type i = vbegin; *i <= k; ++i)
00160               prd_primes *= *i;
00161             // Now the main loop:
00162             for (v_it_type pi = vbegin; pi != vend; ++pi) {
00163               const vertex_type p = *pi;
00164               if (p <= k) continue;
00165               set_system_type F; // collects all progressions ending in p
00166               if (p > prd_primes) {
00167                 const vertex_type max_slope = (p-k)/(k-1);
00168                 for (vertex_type d = prd_primes; d <= max_slope; d+=prd_primes) {
00169                   bool all_primes = true;
00170                   for (vertex_type i = 1, e = p-d; i < k; ++i, e -= d)
00171                     if (not primes_table[e]) {all_primes = false; break;}
00172                   if (all_primes) {
00173                     // TODO: use progression-iterator here
00174                     hyperedge_type H;
00175                     H.reserve(k);
00176                     for (vertex_type i = 0, e = p - (k-1)*d; i < k; ++i, e += d)
00177                       H.push_back(e);
00178                     F.push_back(H);
00179                   }
00180                 }
00181               }
00182               if (prime_k) { // the progression starting with k
00183                 const vertex_type q = (p-k)/(k-1);
00184                 if (q >= 2 and (k-1)*q == p-k) {
00185                   bool all_primes = true;
00186                   for (vertex_type e = k+q; e <= p; e+= q)
00187                     if (not primes_table[e]) {all_primes = false; break;}
00188                   if (all_primes) {
00189                     // TODO: use progression-iterator here
00190                     hyperedge_type H;
00191                     H.reserve(k);
00192                     for (vertex_type e = k; e <= p; e += q)
00193                       H.push_back(e);
00194                     F.push_back(H);
00195                   }
00196                 }
00197               }
00198               result.insert(result.end(), F.rbegin(), F.rend());
00199             }
00200             return result;
00201           }
00202 
00203           bool prime(const vertex_type e) const {
00204             const unsigned int res = mpz_probab_prime_p(mpz_class(e).get_mpz_t(), rounds_compositeness_test);
00205             if (res == 2) return true;
00206             if (res == 0) return false;
00207             throw std::runtime_error("ERROR[GreenTao]: Uncertain about primality.");
00208           }
00209 
00210         private :
00211 
00212           const hyperedge_type vertex_set_;
00213           const vertex_type n; // number of vertices
00214           const vertex_type max_prime; // maximal value of a vertex
00215 
00216         };
00217 
00218 
00235         template <class Vertices, class Hyperedges>
00236         struct Sizes_strata_indmon {
00237           typedef Vertices vertex_container_type;
00238           typedef Hyperedges hyperedge_container_type;
00239           typedef typename vertex_container_type::value_type vertex_type;
00240           typedef typename vertex_container_type::size_type v_size_type;
00241           typedef typename hyperedge_container_type::size_type h_size_type;
00242 
00243           typedef std::pair<v_size_type, h_size_type> pair_type;
00244           typedef std::vector<pair_type> result_type;
00245           
00246           result_type operator()(const vertex_container_type& V, const hyperedge_container_type& G) {
00247             result_type result;
00248             if (G.empty()) return result;
00249             typedef typename vertex_container_type::const_iterator v_it_t;
00250             v_it_t v_it = V.begin();
00251             vertex_type v = *v_it;
00252             v_size_type i = 1; // i is the index of v
00253             typedef typename hyperedge_container_type::const_iterator h_it_t;
00254             h_it_t h_it = G.begin();
00255             vertex_type m = h_it -> back();
00256             while (v != m) { ++v_it; ++i; v = *v_it; }
00257             h_size_type c = 1;
00258             ++h_it;
00259             const h_it_t h_end = G.end();
00260             for(; h_it != h_end; ++h_it) {
00261               const vertex_type nm = h_it -> back();
00262               if (nm == m) ++c;
00263               else {
00264                 result.push_back(pair_type(i,c));
00265                 m = nm; c = 1;
00266                 while (v != m) { ++v_it; ++i; v = *v_it; }
00267               }
00268             }
00269             result.push_back(pair_type(i,c));
00270             return result;
00271           }
00272         };
00273 
00274         template <class Vertices, class Hyperedges>
00275         inline typename Sizes_strata_indmon<Vertices, Hyperedges>::result_type
00276           sizes_strata_indmon(const Vertices& V, const Hyperedges& H) {
00277           return Sizes_strata_indmon<Vertices, Hyperedges>()(V,H);
00278         }
00279 
00296         template <typename V, typename Num>
00297         struct Accumulate_l {
00298           typedef V label_type;
00299           typedef Num number_type;
00300           typedef std::pair<label_type, number_type> pair_type;
00301           typedef std::vector<pair_type> vector_type;
00302           void operator()(vector_type& vector) {
00303             number_type sum(0);
00304             typedef typename vector_type::iterator it_t;
00305             const it_t end = vector.end();
00306             for (it_t i = vector.begin(); i != end; ++i) {
00307               sum += i -> second;
00308               i -> second = sum;
00309             }
00310           }
00311         };
00312 
00313         template <typename V, typename Num>
00314         inline void accumulate_l(
00315           std::vector<std::pair<V, Num> >& vector) {
00316           Accumulate_l<V,Num>()(vector);
00317         }
00318         
00319       }
00320     }
00321   }
00322 }
00323 
00324 #endif