fork download
copy 
  1. // Classification of total antisymmetric relations
  2. // See http colon slash slash ruediger-plantiko.blogspot.ch/2016/07/jenseits-von-schere-stein-und-papier.html
  3. // And http colon slash slash ruediger-plantiko.net/relations/ for a JS version
  4.  
  5. #include <functional>
  6. #include <iostream>
  7. #include <vector>
  8. #include <cstdint>
  9. #include <math.h>
  10. #include <bitset>
  11. #include <algorithm>
  12. #include <unordered_map>
  13.  
  14. using namespace std;
  15.  
  16. // ----------------------------------------------------------------------------
  17. // class Bitset
  18. // Homebrewn class for a sequence of bits, can do increments (++ overloaded), 
  19. // bit get, set and unset, and comparisons (== overloaded)
  20. // ----------------------------------------------------------------------------
  21. class Bitset {
  22.   vector<uint32_t> data;
  23.   size_t length;
  24.   size_t maxIndex, maxPlaces;
  25.   uint32_t maxMask;
  26. public:
  27.   Bitset(size_t length, uint32_t word0 = 0 ) : length(length)  {
  28.     this->maxIndex = (size_t) (length/32);
  29.     this->maxPlaces = length % 32 == 0 ? 32 : ( length % 32 );
  30.     this->maxMask = this->maxPlaces < 32 ? ( 1 << this->maxPlaces ) - 1 : 0xffffffff;
  31.     this->data.resize( this->maxIndex + 1, 0 );
  32.     this->data[0] = word0;
  33.     for (size_t i=1;i<=this->maxIndex;i++) this->data[i] = 0;
  34.   }
  35.  
  36.   Bitset(const Bitset& bs) {
  37.     this->length = bs.length;
  38.     this->maxIndex = bs.maxIndex;
  39.     this->maxPlaces = bs.maxPlaces;
  40.     this->maxMask = bs.maxMask;
  41.     for (size_t i=0;i<=this->maxIndex;i++) this->data.push_back(bs.data[i]);
  42.   }
  43.  
  44.   void toStream(ostream &str) const {
  45.     uint32_t highestWord = this->data[maxIndex];
  46.     for (uint32_t bit = 1<<(this->maxPlaces-1);bit;bit >>=1) {
  47.       str << ((bit & highestWord) ? '1' : '0');
  48.     }
  49.     for (size_t i=this->maxIndex-1;i>=0;i--) str << bitset<32>(this->data[i]); 
  50.   }  
  51.  
  52.   bool equals(const Bitset& bf) const {
  53.     if (this->length != bf.length) return false;
  54.     for (size_t i=0;i<this->maxIndex; i++) 
  55.       if (this->data[i]!=bf.data[i]) return false;
  56.     return (((
  57.       this->data[this->maxIndex]^bf.data[this->maxIndex])
  58.       &this->maxMask)==0);
  59.   }
  60.  
  61.   bool isZero() {
  62.     for (size_t i=0;i<this->maxIndex;i++) if (this->data[i]!=0) return false;
  63.     return (this->data[this->maxIndex]&this->maxMask)==0;
  64.   }
  65.  
  66.   bool get(const size_t pos) const {
  67.     return (this->data[(size_t)(pos/32)] & (1 << (pos % 32) )) && true;
  68.   }
  69.  
  70.   friend bool operator==(const Bitset& lhs, const Bitset& rhs){ 
  71.     return lhs.equals(rhs); 
  72.   }
  73.   friend bool operator!=(const Bitset& lhs, const Bitset& rhs){ 
  74.     return !(lhs == rhs); 
  75.   }
  76.   void clear() {
  77.     for (size_t i=0;i<=this->maxIndex;i++) this->data[i] = 0;
  78.   }
  79.   void set(const size_t pos) {
  80.     this->data[(int)(pos/32)] |= (1 << (pos % 32));
  81.   }
  82.   void unset(const size_t pos) {
  83.     this->data[(int)(pos/32)] &= ~(1 << (pos % 32));
  84.   }
  85.   Bitset& operator++() {
  86.     for (size_t i=0;i<this->maxIndex;i++) {
  87.       if (this->data[i] != 0xffffffff) {
  88. 	this->data[i]++;
  89.         return *this;
  90.       }
  91.       else {
  92. 	this->data[i] = 0;      
  93.       }    
  94.     }
  95.     if (this->data[this->maxIndex] != this->maxMask) {
  96.       this->data[this->maxIndex]++;
  97.     }
  98.     else {
  99.       this->data[this->maxIndex] = 0;
  100.     }
  101.     return *this;
  102.   }
  103.   Bitset operator++(int) {
  104.     Bitset tmp(*this);
  105.     operator++();
  106.     return tmp;
  107.   }
  108.  
  109. };
  110.  
  111. ostream & operator<<(ostream & Str, Bitset const & b) { 
  112.   b.toStream(Str);
  113.   return Str;
  114. }
  115.  
  116. // ----------------------------------------------------------------------------
  117. // Auxiliary class IndexComputer
  118. // Enumerate all the places of the upper triangle of a N by N matrix
  119. // For each N, there is precisely one instance of IndexComputer
  120. // ----------------------------------------------------------------------------
  121. class IndexComputer {
  122.   typedef struct { int i,j; } pair;
  123.   int N;
  124.   int C,C2;
  125.   double C1,D;
  126.   IndexComputer(int N): N(N) {
  127.     C = 2*N-1;
  128.     C1 = N - 0.499999;
  129.     C2 = C - 2;
  130.     D = C*C/4;
  131.   }
  132.   static vector<IndexComputer*> inst;
  133. public:
  134.   static const IndexComputer* const get(int N) {
  135.     if (inst.size()<N+1){
  136.       for (int i=inst.size();i<N+1;i++)
  137.         inst.push_back( new IndexComputer(i));
  138.     }
  139.     return inst[N];
  140.   } 
  141.   const pair getPair(const int index) const {
  142.     int k = (int)(C1-sqrt(D-2.*index));
  143.     return {k,index+(k-C2)*k/2+1};
  144.   } 
  145.   const int getIndex(const int i,const int j) const {
  146.     return i*N-(i+1)*(i+2)/2+j;
  147.   }
  148. };
  149.  
  150. vector<IndexComputer*> IndexComputer::inst = {};
  151.  
  152. // Forward declaration of function permute()
  153. // Needed in DominanceRelation::EquivalenceCheck
  154. void permute(vector<int> numbers,function<bool(const vector<int>&)> f);
  155.  
  156. // ----------------------------------------------------------------------------
  157. // DominanceRelation is the main class of this program
  158. // An instance represents a total antisymmetric relation
  159. // The key part of the implementation is the check of 
  160. // two instances against each other for equivalence 
  161. // (i.e. is there a permutation of indices transforming the
  162. // first into the second DominanceRelation )
  163. // ----------------------------------------------------------------------------
  164. class DominanceRelation {
  165.  
  166.   const int N, NN;
  167.   const IndexComputer *ic;
  168.   vector<vector<int>> cardinality;
  169.   vector<vector<int>> cycles;
  170. public:
  171.   Bitset* bits;
  172.   DominanceRelation(const int N,Bitset* bits):N(N),NN(N*(N-1)/2),bits(bits) {
  173.     ic = IndexComputer::get(N);
  174.   }
  175.   DominanceRelation(const int N,const function<bool(int,int)> rel):N(N),NN(N*(N-1)/2) {
  176.     ic = IndexComputer::get(N);
  177.     bits = new Bitset(NN);
  178.     for (int i=0;i<NN;i++) {
  179.       auto pair = ic->getPair(i);
  180.       if (rel(pair.i,pair.j)) bits->set(i);
  181.     }
  182.   }
  183.   void bitsChanged() {
  184.     cardinality.resize(0);
  185.     cycles.resize(0);
  186.   }
  187.   const bool isEquivalent( DominanceRelation& r) {
  188.     return EquivalenceCheck(*this,r).areEquivalent();
  189.   }
  190.   const string getCardinalitySignature() {
  191.     const vector<vector<int>>* c = getCardinality();
  192.     vector<int> sizes(N);
  193.     for (int i=0;i<N;i++) sizes[i] = (c->at(i)).size();
  194.     return getSignature( &sizes );
  195.   }
  196.   const vector<vector<int>>* const getCardinality() {
  197.     if (cardinality.size()==0) {
  198.       cardinality.resize(N);
  199.       vector<vector<int>> greaterThan(N);
  200.       for (int i=0;i<NN;i++) {
  201.         auto p = ic->getPair(i);
  202.         if (bits->get(i)) greaterThan[p.i].push_back(p.j);
  203.         else greaterThan[p.j].push_back(p.i); 
  204.       }
  205.       for (int i=0;i<N;i++) {
  206.         cardinality[greaterThan[i].size()].push_back(i);
  207.       }
  208.     } 
  209.     return &cardinality;    
  210.   }
  211.  
  212.   const bool holds(const int i, const int j) const {
  213.     return i == j ? false : 
  214.       i > j ? !bits->get(ic->getIndex(j,i)) : bits->get(ic->getIndex(i,j));
  215.   }
  216.  
  217.   const string getCyclesSignature() {
  218.     const vector<vector<int>>* c = getCycles();
  219.     vector<int> sizes(N+1,0);
  220.     for (int i=0;i<c->size();i++)
  221.       sizes[(c->at(i)).size()]++;
  222.     return getSignature( &sizes );
  223.   }   
  224.  
  225.   const vector<vector<int>>* getCycles() {
  226.     if (cycles.size()==0) {
  227.       vector<vector<int>> _cycles; 
  228.       for (int i=0;i<N;i++) _getCycles(vector<int>(0),i,&_cycles);
  229.       for (int i=0;i<_cycles.size();i++) {
  230. // Ignore cyclic permutations
  231.         if (min_element(_cycles[i].begin(),_cycles[i].end())==_cycles[i].begin()) {
  232.           cycles.push_back(_cycles[i]); 
  233.         }  
  234.       }
  235.     }
  236.     return &cycles;
  237.   }
  238.  
  239.   void toStream(ostream & str) const {
  240.     bool doSpace = false;
  241.     str << '{';
  242.     for (int i=0;i<N;i++)
  243.       for (int j=i+1;j<N;j++) 
  244.         if (holds(i,j)) {
  245.           if (doSpace) str << ' ';
  246.           str << i << ':' << j;
  247.           doSpace = true;
  248.         }
  249.     str << '}'; 
  250.   }
  251.  
  252.  
  253. private: 
  254.   class EquivalenceCheck {
  255.     int N,NN;
  256.     DominanceRelation & r,& s;
  257.     vector<int> target;
  258.     vector<int>* targetslice;
  259.     vector<vector<int>>& rCard,& sCard;
  260.     vector<vector<int>> possible = { vector<int>(0) };
  261.     vector<vector<int>> new_possible;
  262.     int plen,hlen;
  263.   public:
  264.     EquivalenceCheck(DominanceRelation& s,DominanceRelation& r)
  265.         :r(r),s(s),N(s.N),NN(s.NN),sCard(s.cardinality),rCard(r.cardinality) { }
  266.     bool areEquivalent(vector<int>* perm = 0) {
  267.  
  268.       if (s.N!=r.N) return false;
  269.  
  270.       if (sCard.size()==0) s.getCardinality();
  271.       if (rCard.size()==0) r.getCardinality();
  272.  
  273.       for (int i=0;i<N;i++) 
  274. 	    if (sCard[i].size()!=rCard[i].size())
  275. 	      return false;
  276.  
  277.       buildTargetFromCardinality();
  278.  
  279.       function<bool(const vector<int>&)> findNewPossible = [this](const vector<int>& p){
  280.         // Check the new digits against themselves
  281.          if (!this->equalOnSubgroup(p)) return false;
  282.          // Seems to work - combine with formerly proven heads
  283.          this->combineNewWithOld(p);
  284.          return false;
  285.       };
  286.  
  287.       for (int i=0;i<N;i++) {
  288.         targetslice = & rCard[i];
  289.         plen = targetslice->size( );
  290.         if (plen == 0) continue;
  291.         hlen = possible[0].size(); 
  292.         new_possible.resize(0);
  293.         permute( sCard[i], findNewPossible);
  294.         if (new_possible.size()==0) return false;
  295.         possible = new_possible; 
  296.       }
  297.  
  298.       // We made it! 
  299.       if (perm) { // if the permutation is requested
  300.         perm->resize(N,0);
  301.         vector<int> solution = possible[0];
  302.         for (int i=0;i<N;i++) {
  303.           (*perm)[solution[i]] = target[i];
  304.         }
  305.       }
  306.  
  307.       return true;
  308.  
  309.     }
  310.   private:
  311.     void buildTargetFromCardinality(  ){
  312.       target.resize(0);
  313.       for (int i=0;i<N;i++) {
  314. 	vector<int>* s = &r.cardinality[i]; 
  315. 	int n = s->size();
  316. 	for (int j=0;j<n;j++) target.push_back((*s)[j]);
  317.       }
  318.     }
  319.     void combineNewWithOld(const vector<int>& p) {
  320.       for (int i=0,posslen=possible.size();i<posslen;i++){
  321.         if (equalsCombination(possible[i],hlen,p,plen)) {
  322.           vector<int> slice(possible[i]);
  323.           slice.insert(slice.end(),p.begin(),p.end());
  324.           new_possible.push_back(slice);
  325.         }
  326.       }
  327.     }
  328.     bool equalOnSubgroup(const vector<int>& p) {
  329.       for (int i=0;i<plen;i++) {
  330.         for (int j=i+1;j<plen;j++) {
  331.           if (!equals(p[i],p[j],targetslice->at(i),targetslice->at(j)))
  332.             return false;
  333.         }
  334.       }
  335.       return true;
  336.     }
  337.     bool equalsCombination(const vector<int>& head, const int hlen,
  338.                            const vector<int>& tail, const int tlen) {
  339.       for (int j=0;j<hlen;j++) {
  340. 	for (int k=0;k<tlen;k++) {
  341. 	  if (!equals(head[j],tail[k],target[j],targetslice->at(k))) {
  342. 	    return false;
  343. 	  }
  344. 	}
  345.       }
  346.       return true; 
  347.     }
  348.     bool equals(const int i1,const int j1,const int i2,const int j2) {
  349.       return s.holds(i1,j1) == r.holds(i2,j2);
  350.     }
  351.   }; 
  352.   string getSignature(const vector<int>* v) {
  353.     string retval = "(";
  354.     bool doSpace = false;
  355.     int vsize = v->size();
  356.     for (int i=0; i<vsize;i++) {
  357.       int val = v->at(i);
  358.       if (val) {
  359.         if (doSpace) retval += ' ';
  360.         retval += (to_string(i) + ':' + to_string(val));
  361.         doSpace = true;
  362.       }
  363.     }
  364.     retval += ")";
  365.     return retval;
  366.   } 
  367.   void _getCycles(vector<int> a,int next,vector<vector<int>>* _cycles) {
  368.     if (find(a.begin(),a.end(),next)!=a.end()) return;
  369.     int size = a.size();
  370.  
  371.     if (size > 0 && holds(next,a[0])) {
  372.       vector<int> cycle = a;
  373.       cycle.push_back(next);
  374.       _cycles->push_back( cycle );
  375.     }
  376.     if (size >= N) return;
  377.  
  378.     vector<int> slice = a;
  379.     slice.push_back(next);
  380.     for (int i=0;i<N;i++) {
  381.       if (i==next) continue;
  382.       if (holds(next,i)) _getCycles(slice,i,_cycles);    
  383.     }
  384.   }
  385. };
  386.  
  387. ostream & operator<<(ostream & Str, DominanceRelation const & d) { 
  388.   d.toStream(Str);
  389.   return Str;
  390. }
  391.  
  392. // ----------------------------------------------------------------------------
  393. // Permute the elements of an array and call a callback for each permutation
  394. // idea from le_m http://stackoverflow.com/questions/9960908/permutations-in-javascript/37580979#37580979
  395. // ----------------------------------------------------------------------------
  396. void permute(vector<int> numbers,function<bool(const vector<int>&)> f) {
  397.  
  398.   int n = numbers.size();
  399.   vector<int> c(n);
  400.   int i = 1, j = 1;
  401.  
  402.   if (f(numbers)) return;
  403.   while (i < n) {
  404.     if (c[i] < i) {
  405.       int k = (i % 2) ? c[i] : 0,
  406.           p = numbers[i];
  407.       numbers[i] = numbers[k];
  408.       numbers[k] = p;
  409.       c[i]++;
  410.       i = 1;
  411.       if(f(vector<int>(numbers)))return;
  412.     } else {
  413.       c[i] = 0;
  414.       i++;
  415.     }
  416.   }  
  417. }
  418.  
  419.  
  420.  
  421. // ----------------------------------------------------------------------------
  422. // Main program
  423. // Loop over all 2^N*(N-1)/2 dominance relations
  424. // Collect mutually inequivalent dom. rel. in a hash of arrays,
  425. // with the cardinality as key
  426. // ----------------------------------------------------------------------------
  427. int main(int argc, char**argv) {
  428.  
  429.   typedef struct {
  430.     DominanceRelation* rel;
  431.     int size;
  432.   } t_entry;
  433.  
  434.  
  435.   string input;
  436.   getline(cin,input);
  437.   int N = atoi(input.c_str()); 
  438.   if (N==0) {
  439.     cerr << "stdin should contain the number N > 0" << endl;
  440.     return 8;
  441.   }
  442.  
  443.   Bitset bs(N*(N-1)/2);
  444.   DominanceRelation r(N,&bs);
  445.  
  446.   unordered_map<string,vector<t_entry>> classes;
  447.  
  448.   do {
  449.     string cs = r.getCardinalitySignature();
  450.     bool found = false;
  451.     vector<t_entry>& c = classes[cs];
  452.     for (t_entry& f : c) {
  453.       if (f.rel->isEquivalent(r)) {
  454.         found = true;
  455.         f.size+=2;
  456.       } 
  457.     }
  458.     if (!found) {
  459.       auto ne = (t_entry) {
  460.         .rel  = new DominanceRelation(N,new Bitset(*(r.bits))),
  461.         .size = 2,
  462.       };
  463.       c.push_back(ne);
  464.     }
  465.  
  466.  
  467. // Increment by 2, since the last bit can be determined arbitrarily
  468. // up to equivalence
  469.     bs++;
  470.     if (bs.isZero()) break;
  471.     bs++;
  472.     r.bitsChanged();
  473.   } while( !bs.isZero() );
  474.  
  475.   int total = 0;
  476.   for (auto pair : classes) {
  477.     cout << pair.first << endl;
  478.     for (auto e : pair.second) {
  479.       cout << "  " << *(e.rel) << endl;
  480.       total ++;
  481.     }
  482.   }
  483.   cout << total << endl;
  484.  
  485.  
  486.   return 0;
  487. }
Success #stdin #stdout 0.37s 3440KB
comments (?)
stdin
copy 
6
stdout
copy 
(0:1 3:5)
  {0:2 0:4 0:5 1:2 1:5}
(0:1 1:1 2:1 4:3)
  {0:2 0:3 0:4 0:5 1:2 1:3 1:4}
(0:1 2:1 3:3 4:1)
  {0:2 0:4 0:5 1:2}
  {0:2 0:5 1:2 1:4}
  {0:2 0:3 0:5 1:2 1:4}
(0:1 1:1 3:3 5:1)
  {0:2 0:3 0:4 1:2 1:3}
(0:1 2:2 3:2 5:1)
  {0:2 0:4 1:2}
(1:2 2:2 4:1 5:1)
  {0:3}
(2:3 3:3)
  {0:4 0:5 1:5}
  {0:4 0:5 1:3 1:5}
  {0:3 0:5 1:4 1:5}
  {0:4 0:5 1:4 1:5}
  {0:3 0:5 1:2 1:5 2:4}
(1:3 3:1 4:1 5:1)
  {0:2}
(1:2 2:1 3:2 5:1)
  {0:4}
  {0:2 0:4}
(1:1 2:3 4:2)
  {0:3 0:5}
  {0:5 1:3}
  {0:2 0:5 1:3}
(2:5 5:1)
  {0:3 0:4 1:4}
(1:1 2:2 3:2 4:1)
  {0:4 0:5}
  {0:3 0:4 0:5}
  {0:2 0:4 0:5 1:3}
  {0:5 1:4}
  {0:2 0:5 1:4}
  {0:2 0:3 0:5 1:4}
  {0:2 0:4 0:5 1:4}
  {0:3 0:5 1:2 1:4}
  {0:3 0:4 0:5 1:2 1:4}
  {0:2 0:5 1:3 1:4}
  {0:2 0:4 0:5 1:3 1:4}
  {0:3 0:4 1:2 1:5}
(1:3 4:3)
  {0:3 1:2 1:3 1:5}
(1:1 2:1 3:4)
  {0:2 0:4 0:5 1:5}
  {0:3 0:4 0:5 1:2 1:5}
  {0:2 0:4 0:5 1:3 1:5}
  {0:2 0:3 0:5 1:4 1:5}
(2:4 3:1 4:1)
  {0:4 0:5 1:3}
  {0:3 0:5 1:4}
  {0:4 0:5 1:4}
  {0:3 0:4 0:5 1:4}
(0:1 1:1 2:1 3:1 4:1 5:1)
  {}
(1:1 2:3 3:1 5:1)
  {0:3 0:4}
  {0:4 1:3}
  {0:2 0:4 1:3}
(0:1 2:2 3:1 4:2)
  {0:2 0:5 1:2}
  {0:2 0:3 0:5 1:2}
(1:2 2:1 3:1 4:2)
  {0:5}
  {0:2 0:5}
  {0:2 0:3 0:5}
  {0:2 0:3 0:5 1:3}
(0:1 1:1 3:2 4:2)
  {0:2 0:3 0:5 1:2 1:3}
(1:2 3:3 4:1)
  {0:2 0:4 0:5}
  {0:5 1:2 1:4}
  {0:2 0:3 0:5 1:3 1:4}
(0:1 2:3 4:1 5:1)
  {0:2 0:3 1:2}
56
https://ideone.com/KBihCA
Explore dominance relations
(cf. http colon slash slash ruediger-plantiko.blogspot.ch/2016/07/jenseits-von-schere-stein-und-papier.html )
language:
C++14 (gcc 8.3)
created:
8 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!