permlib  0.2.9
Library for permutation computations
bsgs.h
1 // ---------------------------------------------------------------------------
2 //
3 // This file is part of PermLib.
4 //
5 // Copyright (c) 2009-2011 Thomas Rehn <thomas@carmen76.de>
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions
10 // are met:
11 // 1. Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 // 2. Redistributions in binary form must reproduce the above copyright
14 // notice, this list of conditions and the following disclaimer in the
15 // documentation and/or other materials provided with the distribution.
16 // 3. The name of the author may not be used to endorse or promote products
17 // derived from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
21 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
22 // IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
23 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
24 // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
28 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // ---------------------------------------------------------------------------
31 
32 
33 #ifndef BSGS_H_
34 #define BSGS_H_
35 
36 #include <map>
37 #include <list>
38 #include <vector>
39 
40 #include <boost/cstdint.hpp>
41 #include <boost/foreach.hpp>
42 #include <boost/next_prior.hpp>
43 #include <boost/scoped_ptr.hpp>
44 #include <boost/shared_ptr.hpp>
45 #include <boost/utility.hpp>
46 
47 #include <permlib/bsgs_core.h>
48 
49 #include <permlib/transversal/orbit_set.h>
50 #include <permlib/transversal/transversal.h>
51 #include <permlib/predicate/pointwise_stabilizer_predicate.h>
52 #include <permlib/predicate/stabilizes_point_predicate.h>
53 
54 #include <permlib/redundant_base_point_insertion_strategy.h>
55 
56 namespace permlib {
57 
58 template <class PERM, class TRANS>
59 struct BSGS;
60 
61 template <class PERM, class TRANS>
62 std::ostream &operator<< (std::ostream &out, const BSGS<PERM,TRANS> &bsgs) {
63  out << "BASE[" << bsgs.B.size() << "]" << std::endl;
64  BOOST_FOREACH(unsigned long beta, bsgs.B) {
65  out << static_cast<unsigned int>(beta+1) << ",";
66  }
67  out << std::endl;
68  out << "SGS[" << bsgs.S.size() << "]" << std::endl;
69  BOOST_FOREACH(const typename PERM::ptr &g, bsgs.S) {
70  out << *g << ",";
71  }
72  out << std::endl;
73  out << "U" << std::endl;
74  BOOST_FOREACH(const TRANS &U, bsgs.U) {
75  for (unsigned int i=0; i<bsgs.n; ++i)
76  // trigger transversal depth reload
77  boost::scoped_ptr<PERM> dummy(U.at(i));
78  out << U.size() << "{" << U.m_statMaxDepth << "}" << ",";
79  }
80  out << " = " << bsgs.order() << std::endl;
81  BOOST_FOREACH(const TRANS &U, bsgs.U) {
82  out << U << std::endl;
83  }
84  return out;
85 }
86 
88 template <class PERM, class TRANS>
89 struct BSGS : public BSGSCore<PERM,TRANS> {
90  typedef typename BSGSCore<PERM,TRANS>::PERMlist PERMlist;
91 
93  explicit BSGS(dom_int n);
95 
99  BSGS(const BSGS<PERM,TRANS>& bsgs);
100 
102 
107 
109 
112  template<typename Integer>
113  Integer order() const;
114 
116 
119  boost::uint64_t order() const;
120 
122 
127  unsigned int sift(const PERM& g, PERM& siftee, unsigned int j = 0) const;
129 
135  unsigned int sift(const PERM& g, PERM& siftee, unsigned int j, unsigned int k) const;
137  bool sifts(const PERM& g) const;
138 
140 
146  bool chooseBaseElement(const PERM &h, dom_int &beta) const;
148 
153  unsigned int insertRedundantBasePoint(unsigned int beta, unsigned int minPos = 0);
155  void stripRedundantBasePoints(int minPos = 0);
156 
158 
166 
168 
171  void orbit(unsigned int j, const PERMlist &generators);
173 
176  void orbitUpdate(unsigned int j, const PERMlist &generators, const typename PERM::ptr &g);
177 
179 
184  int insertGenerator(const typename PERM::ptr& g, bool updateOrbit);
186 
190  void updateOrbits(int pos);
191 
193 
196  PERM random(const int i = 0) const;
197 
199 
203  void conjugate(const PERM& g);
204 
206  friend std::ostream &operator<< <> (std::ostream &out, const BSGS<PERM,TRANS> &bsgs);
207 private:
209  template <class BaseIterator, class TransversalIterator>
210  unsigned int sift(const PERM& g, PERM& siftee, BaseIterator begin, BaseIterator end, TransversalIterator beginT, TransversalIterator endT) const;
211 
213  static int ms_bsgsId;
214 
216  void copyTransversals(const BSGS<PERM,TRANS>& bsgs);
217 };
218 
219 //
220 // ---- IMPLEMENTATION
221 //
222 
223 template <class PERM, class TRANS>
225 
226 template <class PERM, class TRANS>
228  : BSGSCore<PERM,TRANS>(++ms_bsgsId, n_, 0)
229 {}
230 
231 template <class PERM, class TRANS>
233  : BSGSCore<PERM,TRANS>(bsgs.m_id, bsgs.B, bsgs.U, bsgs.n)
234 {
235  copyTransversals(bsgs);
236 }
237 
238 template <class PERM, class TRANS>
240  if (this == &bsgs)
241  return *this;
242 
243  this->B = bsgs.B;
244  this->n = bsgs.n;
245  this->m_id = bsgs.m_id;
246 
247  copyTransversals(bsgs);
248  return *this;
249 }
250 
251 template <class PERM, class TRANS>
252 template <class BaseIterator, class TransversalIterator>
253 unsigned int BSGS<PERM, TRANS>::sift(const PERM& g, PERM& siftee, BaseIterator begin, BaseIterator end, TransversalIterator beginT, TransversalIterator endT) const{
254  unsigned int k = 0;
255  siftee = g;
256  BaseIterator baseIt;
257  TransversalIterator transIt;
258  for (baseIt = begin, transIt = beginT; baseIt != end && transIt != endT; ++baseIt, ++transIt) {
259  unsigned long b = *baseIt;
260  const TRANS& U_i = *transIt;
261  //std::cout << " ~~~ sift " << siftee << " b" << b << std::endl;
262  boost::scoped_ptr<PERM> u_b(U_i.at(siftee / b));
263  if (u_b == 0)
264  return k;
265  u_b->invertInplace();
266  siftee *= *u_b;
267  ++k;
268  }
269  return k;
270 }
271 
272 template <class PERM, class TRANS>
273 unsigned int BSGS<PERM, TRANS>::sift(const PERM& g, PERM& siftee, unsigned int j) const {
274  return sift(g, siftee, this->B.begin() + j, this->B.end(), this->U.begin() + j, this->U.end());
275 }
276 
277 template <class PERM, class TRANS>
278 unsigned int BSGS<PERM, TRANS>::sift(const PERM& g, PERM& siftee, unsigned int j, unsigned int k) const {
279  return sift(g, siftee, this->B.begin() + j, this->B.begin() + k, this->U.begin() + j, this->U.begin() + k);
280 }
281 
282 template <class PERM, class TRANS>
283 bool BSGS<PERM, TRANS>::sifts(const PERM& g) const {
284  PERM siftee(this->n);
285  unsigned int m = sift(g, siftee);
286  return this->B.size() == m && siftee.isIdentity();
287 }
288 
289 template <class PERM, class TRANS>
290 bool BSGS<PERM, TRANS>::chooseBaseElement(const PERM &h, dom_int &beta) const {
291  for (beta = 0; beta < this->n; ++beta) {
292  if (std::find(this->B.begin(), this->B.end(), beta) != this->B.end())
293  continue;
294  if (h / beta != beta)
295  return true;
296  }
297  return false;
298 }
299 
300 template <class PERM, class TRANS>
301 void BSGS<PERM, TRANS>::orbit(unsigned int j, const PERMlist &generators) {
302  this->U[j].orbit(this->B[j], generators);
303 }
304 
305 template <class PERM, class TRANS>
306 void BSGS<PERM, TRANS>::orbitUpdate(unsigned int j, const PERMlist &generators, const typename PERM::ptr &g) {
307  this->U[j].orbitUpdate(this->B[j], generators, g);
308 }
309 
310 template <class PERM, class TRANS>
311 PERM BSGS<PERM, TRANS>::random(const int i) const {
312  BOOST_ASSERT( i >= 0 );
313  PERM g(this->n);
314  for (int l = this->U.size()-1; l>=i ; --l) {
315  //std::cout << l << " : " << U[l] << " : " << U[l].size() << std::endl;
316  unsigned long beta = *(boost::next(this->U[l].begin(), randomInt(this->U[l].size())));
317  boost::scoped_ptr<PERM> u_beta(this->U[l].at(beta));
318  g *= *u_beta;
319  }
320  return g;
321 }
322 
323 template <class PERM, class TRANS>
324 void BSGS<PERM, TRANS>::conjugate(const PERM& g) {
325  PERM gInv(g);
326  gInv.invertInplace();
327 
328  //
329  // to conjugate a BSGS, all three components (B,S,U) have to be adjusted
330  //
331 
332  // STEP 1: conjugate generating set S
333  BOOST_FOREACH(typename PERM::ptr& p, this->S) {
334  *p ^= gInv;
335  *p *= g;
336  }
337 
338  std::vector<dom_int> oldB(this->B);
339  for (unsigned int i = 0; i < this->U.size(); ++i) {
340  // STEP 2: adapt base B
341  this->B[i] = g / oldB[i];
342  // STEP 3: conjugate transversal U
343  this->U[i].permute(g, gInv);
344  }
345 }
346 
347 template <class PERM, class TRANS>
348 int BSGS<PERM, TRANS>::insertGenerator(const typename PERM::ptr& g, bool updateOrbit) {
349  int pos = 0;
350  for (; static_cast<unsigned int>(pos) < this->B.size(); ++pos) {
351  if (*g / this->B[pos] != this->B[pos])
352  break;
353  }
354 
355  if (static_cast<unsigned int>(pos) == this->B.size()) {
356  dom_int beta;
357  bool newBaseElement __attribute__((unused)) = chooseBaseElement(*g, beta);
358  BOOST_ASSERT( newBaseElement );
359  this->B.push_back(beta);
360  this->U.push_back(TRANS(this->n));
361  }
362 
363  const int insertionPosition = pos;
364  this->S.push_back(g);
365 
366  if (updateOrbit) {
367  bool groupOrderChanged = false;
368  for (; pos >= 0; --pos) {
369  PERMlist orbitGenerators;
370  const unsigned int oldTransversalSize = this->U[pos].size();
371  //std::cout << "INSERT orbit @ " << pos << " : " << g << std::endl;
372  std::copy_if(this->S.begin(), this->S.end(), std::back_inserter(orbitGenerators),
373  PointwiseStabilizerPredicate<PERM>(this->B.begin(), this->B.begin() + pos));
374  if (orbitGenerators.size() > 0) {
375  orbitUpdate(pos, orbitGenerators, g);
376 
377  // group order can only increase by adding generators
378  if (this->U[pos].size() > oldTransversalSize)
379  groupOrderChanged = true;
380  }
381  }
382 
383  if (!groupOrderChanged) {
384  this->S.pop_back();
385  return -1;
386  }
387  }
388 
389  return insertionPosition;
390 }
391 
392 template <class PERM, class TRANS>
394  if (pos < 0)
395  return;
396  for (; pos >= 0; --pos) {
397  PERMlist orbitGenerators;
398  std::copy_if(this->S.begin(), this->S.end(), std::back_inserter(orbitGenerators),
399  PointwiseStabilizerPredicate<PERM>(this->B.begin(), this->B.begin() + pos));
400  if (orbitGenerators.size() > 0)
401  orbit(pos, orbitGenerators);
402  }
403 }
404 
405 template <class PERM, class TRANS>
406 template <typename Integer>
407 Integer BSGS<PERM, TRANS>::order() const {
408  Integer orderValue(1);
409  BOOST_FOREACH(const TRANS &Ui, this->U) {
410  orderValue *= Ui.size();
411  }
412  return orderValue;
413 }
414 
415 template <class PERM, class TRANS>
416 boost::uint64_t BSGS<PERM, TRANS>::order() const {
417  return order<boost::uint64_t>();
418 }
419 
420 
421 template <class PERM, class TRANS>
422 unsigned int BSGS<PERM, TRANS>::insertRedundantBasePoint(unsigned int beta, unsigned int minPos) {
423  PERMlist S_i;
425  int pos = is.findInsertionPoint(beta, S_i);
426  if (pos < 0)
427  return -(pos+1);
428  pos = std::max(static_cast<unsigned int>(pos), minPos);
429 
430  this->B.insert(this->B.begin() + pos, beta);
431  this->U.insert(this->U.begin() + pos, TRANS(this->n));
432  this->U[pos].orbit(beta, S_i);
433  return pos;
434 }
435 
436 template <class PERM, class TRANS>
438  for (int i = this->B.size()-1; i >= minPos; --i) {
439  if (this->U[i].size() <= 1) {
440  if (i == static_cast<int>(this->B.size()-1)) {
441  this->B.pop_back();
442  this->U.pop_back();
443  } else {
444  this->B.erase(this->B.begin() + i);
445  this->U.erase(this->U.begin() + i);
446  }
447  }
448  }
449 }
450 
451 
453 
457 template <class PERM>
458 class StrongGeneratingSetSorter : public std::binary_function<typename PERM::ptr, typename PERM::ptr, bool> {
459 public:
464  template<class InputIterator>
465  StrongGeneratingSetSorter(InputIterator baseBegin, InputIterator baseEnd) : m_base(baseBegin, baseEnd) { }
466 
468  bool operator()(const typename PERM::ptr& p1, const typename PERM::ptr& p2) const {
469  BOOST_FOREACH(const dom_int b, m_base) {
470  if ( p1->at(b) == b && p2->at(b) != b )
471  return true;
472  if ( p1->at(b) != b )
473  return false;
474  }
475  return false;
476  }
477 private:
478  std::vector<dom_int> m_base;
479 };
480 
481 template <class PERM, class TRANS>
483  PERMlist sortedSGS(this->S);
484  sortedSGS.sort(StrongGeneratingSetSorter<PERM>(this->B.begin(), this->B.end()));
485 
486  PERMlist filteredSGS;
488  dom_int oldBaseElement = static_cast<dom_int>(-1);
489  BOOST_FOREACH(const typename PERM::ptr& gen, sortedSGS) {
490  if (gen->isIdentity())
491  continue;
492  filteredSGS.push_back(gen);
493 
494  // Compute to which base element this strong generator belongs.
495  // That is, gen stabilizes all base points up to baseElement.
496  // No generator can stabilize all base elements because we excluded
497  // identity permutations before.
498  dom_int baseElement = this->B.front();
499  BOOST_FOREACH(const dom_int b, this->B) {
500  baseElement = b;
501  if (*gen / b != b)
502  break;
503  }
504  PERMLIB_DEBUG(std::cout << "gen " << *gen << " @ " << baseElement << std::endl;)
505 
507  newOrbit->orbit(baseElement, filteredSGS, typename Transversal<PERM>::TrivialAction());
508  if (oldBaseElement == baseElement && newOrbit->size() == oldOrbit->size()) {
509  delete newOrbit;
510  PERMLIB_DEBUG(std::cout << " removed\n";)
511  filteredSGS.pop_back();
512  } else {
513  delete oldOrbit;
514  oldOrbit = newOrbit;
515  }
516  oldBaseElement = baseElement;
517  }
518  delete oldOrbit;
519 
520  this->S = filteredSGS;
521 }
522 
523 template <class PERM, class TRANS>
525  std::map<PERM*,typename PERM::ptr> genMap;
526  BOOST_FOREACH(const typename PERM::ptr& p, bsgs.S) {
527  typename PERM::ptr deepcopy(new PERM(*p));
528  //std::cout << "found " << p.get() << " = " << *p << std::endl;
529  genMap.insert(std::make_pair(p.get(), deepcopy));
530  this->S.push_back(deepcopy);
531  }
532 
533  BOOST_ASSERT(this->B.size() == bsgs.B.size());
534  BOOST_ASSERT(bsgs.B.size() == bsgs.U.size());
535  this->U.clear();
536  this->U.resize(bsgs.U.size(), TRANS(bsgs.n));
537  BOOST_ASSERT(this->U.size() == bsgs.U.size());
538 
539  for (unsigned int i=0; i<this->U.size(); ++i) {
540  this->U[i] = bsgs.U[i].clone(genMap);
541  }
542 }
543 
544 }
545 
546 #endif // -- BSGS_H_
stores an orbit in a set for fast contains() operation
Definition: orbit_set.h:42
size_t size() const
number of orbit elements
Definition: orbit_set.h:60
predicate matching a permutation if it stabilizes a given list of points pointwise
Definition: pointwise_stabilizer_predicate.h:42
Class that can be used to sort a strong generating set.
Definition: bsgs.h:458
StrongGeneratingSetSorter(InputIterator baseBegin, InputIterator baseEnd)
Definition: bsgs.h:465
bool operator()(const typename PERM::ptr &p1, const typename PERM::ptr &p2) const
true iff p1 stabilizes more base points (in increasing order) than p2
Definition: bsgs.h:468
insertion position after first non-trivial transversal
Definition: redundant_base_point_insertion_strategy.h:68
virtual int findInsertionPoint(dom_int beta, std::list< typename PERM::ptr > &S_i) const
finds possible insertion point for base point
Definition: redundant_base_point_insertion_strategy.h:73
core data of a base and strong generating set (BSGS)
Definition: bsgs_core.h:42
dom_int n
degree of group
Definition: bsgs_core.h:61
std::vector< TRANS > U
transversals along the stabilizer chain
Definition: bsgs_core.h:59
std::vector< dom_int > B
base
Definition: bsgs_core.h:55
PERMlist S
strong generating set
Definition: bsgs_core.h:57
int m_id
id of this BSGS instance
Definition: bsgs_core.h:81
Represents a base and strong generating set (BSGS)
Definition: bsgs.h:89
PERM random(const int i=0) const
generates a uniformly distributed random element of
Definition: bsgs.h:311
BSGS< PERM, TRANS > & operator=(const BSGS< PERM, TRANS > &)
assignment operator
Definition: bsgs.h:239
void orbitUpdate(unsigned int j, const PERMlist &generators, const typename PERM::ptr &g)
updates the j-th fundamental orbit with the given orbit generators and a new generator g
Definition: bsgs.h:306
unsigned int insertRedundantBasePoint(unsigned int beta, unsigned int minPos=0)
inserts a redundant base beta
Definition: bsgs.h:422
void conjugate(const PERM &g)
conjugate group with a permutation
Definition: bsgs.h:324
BSGS(dom_int n)
constructs an empty group of degree n
Definition: bsgs.h:227
unsigned int sift(const PERM &g, PERM &siftee, unsigned int j=0) const
sifts an element through the specified transversal range
Definition: bsgs.h:273
void stripRedundantBasePoints(int minPos=0)
strips redundant base points from the end to the minPos-th base element
Definition: bsgs.h:437
bool sifts(const PERM &g) const
true iff g sifts through transversal system
Definition: bsgs.h:283
int insertGenerator(const typename PERM::ptr &g, bool updateOrbit)
adds a new group generator
Definition: bsgs.h:348
void orbit(unsigned int j, const PERMlist &generators)
re-computes the j-th fundamental orbit with the given orbit generators
Definition: bsgs.h:301
bool chooseBaseElement(const PERM &h, dom_int &beta) const
tries to find a new base element
Definition: bsgs.h:290
void stripRedundantStrongGenerators()
removes redundant generators
Definition: bsgs.h:482
void updateOrbits(int pos)
updates transversals/orbits
Definition: bsgs.h:393
Integer order() const
order of the group
Definition: bsgs.h:407
action of a PERM on unsigned long element
Definition: transversal.h:112