//====================================================================== // PROJECT: Sally's problem // FILE: perms.cc // PURPOSE: permutation and inversion utilities // UPDATE HISTORY: version 1.0 3/15/10 //====================================================================== //----------------------------- Includes ------------------------------- #include <iostream> #include <vector> #include <list> #include <iterator> #include <algorithm> #include <utility> #include <cassert> using namespace std; //--------------------- Pair sort utilities ------------------------ bool operator== (const pair<int, int> & p1, const pair<int, int> & p2) { return (p1.first == p2.first) && (p1.second == p2.second) || (p1.first == p2.second) && (p1.second == p2.first); } bool operator< (const pair<int, int> & p1, const pair<int, int> & p2) { // First, sort the pairs pair<int, int> q1 = p1; pair<int, int> q2 = p2; if (q1.first > q1.second) { int temp = q1.first; q1.first = q1.second; q1.second = temp; } if (q2.first > q2.second) { int temp = q2.first; q2.first = q2.second; q2.second = temp; } // Now compare in lexicographic order. if (q1.first < q2.first) { return true; } if ((q1.first == q2.first) && (q1.second < q2.second)) { return true; } return false; } //================= Permutation CLASS DECLARATION ================= /****************************************************************** A permutation is represented internally as a vector of ints <p[0], p[1], ... , p[n - 1]>, where each p[i] is a unique integer in the range 0 .. n - 1. This is different from the usual representation starting from 1 and ending with n, but we all know that computer scientists start with zero. This class declares the Inversion class as a friend, for ease of implementation. Along with the usual constructors, we allow a permutation to be constructed from an inversion list. We supply ==, !=, and < as member functions for use in global sorting routines. In addition we supply a global output routine for printing permutations. ******************************************************************/ class Permutation { public: friend class Inversion; // Construct an n-element identity permutation Permutation(int n); // Construct a permutation from a vector of ints // NOTE: this doesn't check whether the vector is a permutation Permutation(vector<int> v) : place(v) { } // Copy ctor Permutation(const Permutation & p); // Construct a permutation from an inversion list Permutation(const Inversion & inv); int size() const { return place.size(); } int at(int n) const { return place.at(n); } // RETURN the order of this permutation int ord() const; // RETURN true iff self and p are the same permutation bool operator== (const Permutation & p) const; // RETURN false iff self and p are the same permutation bool operator!= (const Permutation & p) const { return !(*this == p); } // RETURN true if self is lexicographically less than p // PRE: the perms are over the same base set bool operator< (const Permutation & p) const; // RETURN a copy of the product permutation self * p Permutation operator* (const Permutation & p) const; // RETURN the result of self applied to n int apply_to(int n) const { assert(n >= 0 && n < place.size()); return place[n]; } // RETURN the inverse of this permutation Permutation inverse(); private: vector<int> place; }; //=================== Eq_class CLASS DECLARATION ================== /****************************************************************** An equivalence class object is a vector of Permutations. In this program (though it doesn't have to be), two permutations \pi and \sigma in S_n are equivalent if there is a permutation \tau such that \tau applied to the inversion list of \pi gives the inversion list of \sigma. For example, if \pi = (1 0 2 3) with inversion list (0, 1) then \tau = (1 2 0 3) applied to (0, 1) gives (1, 2), which is the inversion list of (0 2 1 3) so (1 0 2 3) \equiv (0 2 1 3). This condition also requires that the application of \tau yields a list that is ordered: the (i, j) in the transformed list must all have i < j or must all have i > j. As with Permutations, we supply ==, and < as member functions for use in global sorting routines and supply a global output routine for printing permutations. ******************************************************************/ class Eq_class { public: Eq_class() : data() {} // PRE: the vector is sorted Eq_class(const vector<Permutation> & v) : data(v) {} // Put a new permutation into this equivalence class void append(const Permutation & p) { data.push_back(p); } // Merge two equivalence classes into one void join(const Eq_class & c); // RETURN the number of permutations in this equivalence class int size() const { return data.size(); } // RETURN the i-th permutation in this equivalence class. Permutation at(int i) const { return data.at(i); } // PRE: both equivalence classes are sorted by the underlying // permutation orders. bool operator== (const Eq_class & ec) const; // PRE: both equivalence classes are sorted by the underlying // permutation orders. bool operator< (const Eq_class & ec) const; private: vector<Permutation> data; }; //=================== Inversion CLASS DECLARATION =================== /****************************************************************** An Inversion object is a vector of pairs, representing the inversions of a given permutation. This class declares Permutation as a friend, so that we can construct an Inversion from its Permutation. ******************************************************************/ class Inversion { public: friend class Permutation; // Construct an empty list of inversions Inversion() : num(0), pairs() { } // Copy ctor Inversion(const Inversion & inv) : num(inv.num), pairs(inv.pairs) { } // Construct an inversion list from a permutation Inversion(const Permutation & p); bool operator== (const Inversion & inv) const; bool operator!= (const Inversion & inv) const { return !(*this == inv); } // RETURN the result of applying a perm to to this inversion list Inversion apply(const Permutation & p) const; // RETURN the inversion vector corresponding to this inversion list vector<int> inversion_vector() const; // RETURN true iff Sally's criterion is met, i.e., every inversion // pair is ordered (first < second) or (first > second) bool pair_ordered() const; // RETURN true iff this inversion is that of some permutation bool is_valid() const; // RETURN the number of pairs in this inversion int size() const { return pairs.size(); } // RETURN the i-th pair in this inversion list pair<int, int> at(int i) const { return pairs.at(i); } // RETURN the size of the base set of this inversion list int num_elts() const { return num; } private: int num; vector<pair<int, int> > pairs; Inversion restore_order() const; }; //-------------- Permutation CLASS MEMBER FUNCTIONS ---------------- Permutation::Permutation(int n) : place(n) { for (int i = 0; i < n; i++) { place[i] = i; } } Permutation::Permutation(const Permutation & p) : place(p.size()) { for (int i = 0; i < p.size(); i++) { place[i] = p.place[i]; } } // NOTE: not every list of pairs can result from some permutation // (like (0, 2), (1, 3)). This should check that a list actually // can be the list of a permutation. Permutation::Permutation(const Inversion & inv) : place(inv.num_elts()) { // First, build the inversion vector vector<int> in = inv.inversion_vector(); // Now use the inversion vector to build the permutation list list<int> p; p.push_back(in.size() - 1); list<int>::iterator it; for (int n = in.size() - 2; n >= 0; n--) { int count = 0; it = p.begin(); while (count < in[n]) { if (*it > n) { count++; } it++; } p.insert(it, n); } // Finally, use the list to make the permutation vector int k = 0; for (it = p.begin(); it != p.end(); it++, k++) { place[k] = *it; } } int Permutation::ord() const { Permutation id(place.size()); // construct the identity permutation Permutation q = *this; int i = 1; while (q != id) { i++; q = q * (*this); } return i; } bool Permutation::operator== (const Permutation & p) const { if (place.size() != p.size()) { return false; } for (int i = 0; i < place.size(); i++) { if (place[i] != p.place[i]) { return false; } } return true; } bool Permutation::operator< (const Permutation & p) const { assert(place.size() == p.place.size()); for (int i = 0; i < place.size(); i++) { if (place[i] < p.place[i]) { return true; } if (place[i] > p.place[i]) { return false; } } return false; } Permutation Permutation::operator* (const Permutation & p) const { assert(place.size() == p.size()); vector<int> v(place.size()); for (int i = 0; i < place.size(); i++) { v[i] = place[p.place[i]]; } return Permutation(v); } Permutation Permutation::inverse() { vector<int> v(place.size()); for (int i = 0; i < place.size(); i++) { v[place[i]] = i; } return Permutation(v); } //----------------- Eq_class CLASS MEMBER FUNCTIONS ----------------- void Eq_class::join(const Eq_class & c) { for (int i = 0; i < c.data.size(); i++) { data.push_back(c.data[i]); } sort(data.begin(), data.end()); vector<Permutation>::iterator it = unique(data.begin(), data.end()); data.erase(it, data.end()); } bool Eq_class::operator== (const Eq_class & ec) const { if (data.size() != ec.data.size()) { return false; } for (int i = 0; i < data.size(); i++) { if (data[i] != ec.data[i]) { return false; } } return true; } bool Eq_class::operator< (const Eq_class & ec) const { if (data.size() < ec.data.size()) { return true; } if (data.size() > ec.data.size()) { return false; } for (int i = 0; i < data.size(); i++) { if (data[i] < ec.data[i]) { return true; } if (ec.data[i] < data[i]) { return false; } } return false; } //---------------- Inversion CLASS MEMBER FUNCTIONS ---------------- Inversion::Inversion(const Permutation & p) : num(p.size()), pairs() { for (int i = 0; i < num; i++) { for (int j = 0; j < i; j++) { if (p.at(j) > p.at(i)) { pairs.push_back(make_pair(p.at(i), p.at(j))); } } } } bool Inversion::operator== (const Inversion & inv) const { if (pairs.size() != inv.pairs.size()) { return false; } vector<pair<int, int> > v1 = pairs; vector<pair<int, int> > v2 = inv.pairs; sort(v1.begin(), v1.end()); sort(v2.begin(), v2.end()); for (int i = 0; i < pairs.size(); i++) { if (v1[i] != v2[i]) { return false; } } return true; } Inversion Inversion::apply(const Permutation & p) const { Inversion r = *this; for (int i = 0; i < r.size(); i++) { r.pairs[i].first = p.apply_to(pairs[i].first); r.pairs[i].second = p.apply_to(pairs[i].second); } return r; } vector<int> Inversion::inversion_vector() const { vector<int> v(num); for (int i = 0; i < num; i++) { v[i] = 0; for (int j = 0; j < pairs.size(); j++) { int min = pairs[j].first; if (pairs[j].second < min) { min = pairs[j].second; } if (min == i) { v[i]++; } } } return v; } bool Inversion::pair_ordered() const { if (pairs.size() == 0) { return true; } bool less = true; if (pairs[0].first > pairs[0].second) { less = false; } for (int i = 1; i < pairs.size(); i++) { if (((pairs[i].first < pairs[i].second) && !less) || ((pairs[i].first > pairs[i].second) && less)) { return false; } } return true; } Inversion Inversion::restore_order() const { Inversion ret = *this; for (int i = 0; i < pairs.size(); i++) { if (ret.pairs[i].first > ret.pairs[i].second) { swap(ret.pairs[i].first, ret.pairs[i].second); } } return ret; } // Uses Sally's result: a list of pairs is the inversion list of a permutation // iff for all i < j < k // (1) The list is transitive: (i, j) and (i, k) in the list ==> (i, k) is too // (2) The median property: (i, k) is in ==> either (i, j) or (j, k) is too. bool Inversion::is_valid() const { Inversion r = restore_order(); bool b1, b2; pair<int, int> p; // First, check for transitivity for (int i = 0; i < r.pairs.size(); i++) { b1 = b2 = true; for (int j = i + 1; j < r.pairs.size(); j++) { if (r.pairs[i].second == r.pairs[j].first) { // See if (r.pairs[i].first, r.pairs[j].second) in r.pairs p = make_pair(r.pairs[i].first, r.pairs[j].second); b1 = find(r.pairs.begin(), r.pairs.end(), p) != r.pairs.end(); } if (r.pairs[i].first == r.pairs[j].second) { // See if (r.pairs[i].second, r.pairs[j].first) in r.pairs make_pair(r.pairs[i].second, r.pairs[j].first); b2 = find(r.pairs.begin(), r.pairs.end(), p) != r.pairs.end(); } if (!(b1 || b2)) { return false; } } } // Then check each pair for the median property for (int i = 0; i < r.pairs.size(); i++) { b1 = b2 = true; for (int j = r.pairs[i].first + 1; j < r.pairs[i].second; j++) { // See if (r.pairs[i].first, j) in ir.pairs p = make_pair(r.pairs[i].first, j); b1 = find(r.pairs.begin(), r.pairs.end(), p) != r.pairs.end(); // See if (j, r.pairs[i].second) in ir.pairs p = make_pair(j, r.pairs[i].second); b2 = find(r.pairs.begin(), r.pairs.end(), p) != r.pairs.end(); if (!(b1 || b2)) { return false; } } } return true; } //----------------------- GLOBAL UTILITIES ------------------------ // RETURN a vector of all possible permutations on {0, 1, ... , n - 1} vector<Permutation> make_all_perms(int n) { vector<Permutation> pv; vector<int> v; for (int i = 0; i < n; i++) { v.push_back(i); } Permutation p(v); pv.push_back(p); while (next_permutation(v.begin(), v.end())) { Permutation p(v); pv.push_back(p); } return pv; } //----------------------------------------------------------------------- // This changes from the internal representation of a perm on {0, ..., n-1} // to the customary notation on {1, ..., n} ostream & operator<< (ostream & os, const Permutation & p) { os << "("; for (int i = 0; i < p.size(); i++) { os << " " << p.at(i) + 1; } os << " )"; return os; } ostream & operator<< (ostream & os, const Eq_class & ec) { os << "["; for (int i = 0; i < ec.size(); i++) { os << " " << ec.at(i); } os << " ]"; return os; } // Ditto. See above ostream & operator<< (ostream & os, const Inversion & inv) { for (int i = 0; i < inv.size(); i++) { os << "(" << inv.at(i).first+1 << ", " << inv.at(i).second+1 << ") "; } return os; } // Ditto. See above void print_pair_list(list<pair<int, int> > ls) { list<pair<int, int> >::iterator it; for (it = ls.begin(); it != ls.end(); it++) { cout << "(" << it->first + 1 << "," << it->second + 1 << ") "; } cout << endl; } //----------------------------------------------------------------------- list<pair<int, int> > make_node_list(const Inversion & inv) { list<pair<int, int> > ret; vector<int> count(inv.num_elts(), 0); for (int i = 0; i < inv.size(); i++) { count[inv.at(i).second]++; } for (int i = 0; i < inv.num_elts(); i++) { ret.push_back(make_pair(i, count[i])); } return ret; } // Utility function for make_class(). This takes a list of nodes and // their in-degrees and returns the list that represents the // digraph with node x removed. list<pair<int, int> > trim( int x, list<pair<int, int> > pl, const Inversion & inv) { // First, kill the (x, -) node from pl list<pair<int, int> >::iterator it; for (it = pl.begin(); it->first != x; it++) ; pl.erase(it); // Now decrement the in-degrees where needed for (int i = 0; i < inv.size(); i++) { if (inv.at(i).first == x) { for (it = pl.begin(); it != pl.end(); it++) { if (it->first == inv.at(i).second) { (it->second)--; } } } } return pl; } // RETURN the equivalence class of all permutations // equivalent to the permutation of the inversion list using Sally's // first criterion. This uses topological sorting on the DAG associated // with the inversion list. // It works, but does a lot of recursion that leads to duplicate results // and it gives a lot of invalid results that don't get removed until // the very end. In fact, this runs in time O(n!). Eq_class make_class(const Inversion & inv, vector<int> perm, list<pair<int, int> > node_list) { Eq_class ret; if (node_list.empty()) { Permutation p(perm); p = p.inverse(); Inversion in = inv.apply(p); if (in.is_valid()) { Permutation q(in); ret.append(q); } return ret; } else { for (list<pair<int, int> >::iterator it = node_list.begin(); it != node_list.end(); it++) { if (it->second == 0) { vector<int> vv = perm; vv.push_back(it->first); ret.join(make_class(inv, vv, trim(it->first, node_list, inv))); } } return ret; } } // RETURN all permutations equivalent to the argument, using both // the permutation and its inverse. Eq_class find_all(Permutation p) { Inversion in1(p); vector<int> w; Eq_class eq1 = make_class(in1, w, make_node_list(in1)); if (p != p.inverse()) { p = p.inverse(); Inversion in2(p); Eq_class eq2 = make_class(in2, w, make_node_list(in2)); eq1.join(eq2); } return eq1; } //----------------------------------------------------------------------- // RETURN a list of all permutations on {0, ... , n-1} with k inversions. // Uses the Effler/Ruskey algorithm list<Permutation> make_perm_list(int N, int n, int k, vector<int> v) { list<Permutation> ret; if ((n == 0) && (k == 0)) { Permutation p(v); ret.push_back(p); return ret; } int binom = (n < 3) ? 0 : ((n - 1) * (n - 2)) / 2; // Make a vector of permissible candidates vector<int> places; for (int i = 0; i < N; i++) { bool not_in_v = true; for (int j = n; j < N; j++) { if (v[j] == i) { not_in_v = false; } } if (not_in_v) { places.push_back(i); } } // For all legal x that haven't been used yet, make the rest of // the permutation recursively for (int i = 0; i < places.size(); i++) { if ((n - 1 - i <= k) && (k <= binom + n - 1 - i)) { v[n - 1] = places[i]; list<Permutation> l = make_perm_list(N, n - 1, k - n + 1 + i, v); ret.splice(ret.end(), l); } } return ret; } int factorial(int n) { int r = 1; for (int i = 2; i <= n; i++) { r *= i; } return r; } // Build a permutation from std input, checking that it's legal Permutation get_perm() { cout << "Enter perm, ending with 0: "; vector<int> v; int d; cin >> d; while (d > 0) { v.push_back(d - 1); cin >> d; } // Check that the input is valid. Inefficient, but it's unlikely that // the user would enter a thousand-element permutation. bool legal_input = true; for (int i = 0; i < v.size(); i++) { bool found_it = false; for (int j = 0; j < v.size(); j++) { if (i == v[j]) { found_it = true; break; } } legal_input = legal_input && found_it; } assert(legal_input); Permutation p(v); return p; } //============================= MAIN ============================== // This produces all equivalence classes, grouped by size int main() { int base; cout << "N = "; cin >> base; int possible = factorial(base); vector<Permutation> perms = make_all_perms(base); vector<Eq_class> equivs; // the vector of all equivalence classes for (int i = 0; i < possible; i++) { Permutation p = perms[i]; Inversion in(p); vector<Permutation> vp; // the perms in this equivalence class for (int k = 0; k < possible; k++) { Permutation q = perms[k]; Inversion ii = in.apply(q); if (ii.is_valid() && ii.pair_ordered()) { Permutation r(ii); vp.push_back(r); } } // Eliminate duplicates among the vector of permutations. sort(vp.begin(), vp.end()); vector<Permutation>::iterator it = unique(vp.begin(), vp.end()); vp.erase(it, vp.end()); Eq_class ec(vp); equivs.push_back(ec); } sort(equivs.begin(), equivs.end()); vector<Eq_class>::iterator ite = unique(equivs.begin(), equivs.end()); equivs.erase(ite, equivs.end()); // Look at all the equivalence classes. // This produces a lot of output. It can be commented out. for (int i = 0; i < equivs.size(); i++) { cout << Inversion(equivs[i].at(0)).size() << ": \t"; cout << equivs[i] << endl; } cout << endl << endl; // Now we'll count the number of classes of each size vector<int> count(100, 0); for (int i = 0; i < equivs.size(); i++) { count[equivs[i].size()]++; } int total = 0; cout << "size\tcount" << endl; for (int i = 1; i < 100; i++) { if (count[i] != 0) { cout << " " << i << "\t " << count[i] << endl; total += count[i]; } } cout << endl << "Total classes = " << total << endl; return EXIT_SUCCESS; }